{ "cells": [ { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [], "source": [ "import math\n", "import sqlite3\n", "import struct\n", "import datetime\n", "import scipy.fftpack\n", "from scipy import signal as sig\n", "\n", "import matplotlib\n", "from matplotlib import pyplot as plt\n", "from matplotlib import patches\n", "import numpy as np\n", "from scipy import signal, optimize\n", "from tqdm.notebook import tnrange, tqdm" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib widget" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "db = sqlite3.connect('/mnt/c/Users/jaseg/shared/waveform-raspi.sqlite3')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Run 000: 2020-01-31 18:05:24 - 2020-02-01 00:13:45 ( 6:08:21.589, 22126080sp)\n" ] } ], "source": [ "for run_id, start, end, count in db.execute('SELECT run_id, MIN(rx_ts), MAX(rx_ts), COUNT(*) FROM measurements GROUP BY run_id'):\n", " foo = lambda x: datetime.datetime.fromtimestamp(x/1000)\n", " start, end = foo(start), foo(end)\n", " print(f'Run {run_id:03d}: {start:%Y-%m-%d %H:%M:%S} - {end:%Y-%m-%d %H:%M:%S} ({str(end-start)[:-3]:>13}, {count*32:>9d}sp)')\n", "last_run, n_records = run_id, count\n", "sampling_rate = 1000.0" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "par = lambda *rs: 1/sum(1/r for r in rs) # resistor parallel calculation\n", "\n", "# FIXME: These are for the first prototype only!\n", "vmeas_source_impedance = 330e3\n", "vmeas_source_scale = 0.5\n", "\n", "vcc = 15.0\n", "vmeas_div_high = 27e3\n", "vmeas_div_low = par(4.7e3, 10e3)\n", "vmeas_div_voltage = vcc * vmeas_div_low / (vmeas_div_high + vmeas_div_low)\n", "vmeas_div_impedance = par(vmeas_div_high, vmeas_div_low)\n", "\n", "#vmeas_overall_factor = vmeas_div_impedance / (vmeas_source_impedance + vmeas_div_impedance)\n", "v0 = 1.5746\n", "v100 = 2.004\n", "vn100 = 1.1452\n", "\n", "adc_vcc = 3.3 # V\n", "adc_fullscale = 4095\n", "\n", "adc_val_to_voltage_factor = 1/adc_fullscale * adc_vcc\n", "\n", "adc_count_to_vmeas = lambda x: (x*adc_val_to_voltage_factor - v0) / (v100-v0) * 100" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f960454ab93244db97e039ae7ebd02ee", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(FloatProgress(value=0.0, max=691440.0), HTML(value='')))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "limit = n_records\n", "record_size = 32\n", "skip_dropped_sections = False\n", "\n", "data = np.zeros(limit*record_size)\n", "data[:] = np.nan\n", "\n", "last_seq = None\n", "write_index = 0\n", "for i, (seq, chunk) in tqdm(enumerate(db.execute(\n", " 'SELECT seq, data FROM measurements WHERE run_id = ? ORDER BY rx_ts LIMIT ? OFFSET ?',\n", " (last_run, limit, n_records-limit))), total=n_records):\n", " \n", " if last_seq is None or seq == (last_seq + 1)%0xffff:\n", " last_seq = seq\n", " idx = write_index if skip_dropped_sections else i\n", " data[idx*record_size:(idx+1)*record_size] = np.frombuffer(chunk, dtype=' last_seq:\n", " last_seq = seq\n", " # nans = np.empty((record_size,))\n", " # nans[:] = np.nan\n", " # data = np.append(data, nans) FIXME\n", " \n", "data = (data * adc_val_to_voltage_factor - v0) / (v100-v0) * 100" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "227.68691180713367" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_not_nan = data[~np.isnan(data)]\n", "np.sqrt(np.mean(np.square(data_not_nan)))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "264a9f8478e449289c1592c9595dc6cc", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (top, bottom) = plt.subplots(2, figsize=(9,6))\n", "fig.tight_layout(pad=3, h_pad=0.1)\n", "\n", "range_start, range_len = -300, 60 # [s]\n", "\n", "data_slice = data[ int(range_start * sampling_rate) : int((range_start + range_len) * sampling_rate) ]\n", "\n", "top.grid()\n", "top.plot(np.linspace(0, range_len, int(range_len*sampling_rate)), data_slice, lw=1.0)\n", "top.set_xlim([range_len/2-0.25, range_len/2+0.25])\n", "mean = np.mean(data_not_nan)\n", "rms = np.sqrt(np.mean(np.square(data_not_nan - mean)))\n", "peak = np.max(np.abs(data_not_nan - mean))\n", "top.axhline(mean, color='red')\n", "bbox = {'facecolor': 'black', 'alpha': 0.8, 'pad': 2}\n", "top.text(0, mean, f'mean: {mean:.3f}', color='white', bbox=bbox)\n", "top.text(0.98, 0.2, f'V_RMS: {rms:.3f}', transform=top.transAxes, color='white', bbox=bbox, ha='right')\n", "top.text(0.98, 0.1, f'V_Pk: {peak:.3f}', transform=top.transAxes, color='white', bbox=bbox, ha='right')\n", "\n", "bottom.grid()\n", "bottom.specgram(data_slice, Fs=sampling_rate)\n", "top.set_ylabel('U [V]')\n", "bottom.set_ylabel('F [Hz]')\n", "bottom.set_xlabel('t [s]')\n", "None" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "fs = sampling_rate # Hz\n", "ff = 50 # Hz\n", "\n", "analysis_periods = 10\n", "window_len = fs * analysis_periods/ff\n", "nfft_factor = 4\n", "sigma = window_len/8 # samples\n", "\n", "f, t, Zxx = signal.stft(data,\n", " fs = fs,\n", " window=('gaussian', sigma),\n", " nperseg = window_len,\n", " nfft = window_len * nfft_factor)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "7c2382eb8e124ef9b29546ecaed3e155", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(9, 3))\n", "fig.tight_layout(pad=2, h_pad=0.1)\n", "\n", "ax.pcolormesh(t[-200:-100], f[:250], np.abs(Zxx[:250,-200:-100]))\n", "ax.set_title(f\"Run {last_run}\", pad=-20, color='white')\n", "ax.grid()\n", "ax.set_ylabel('f [Hz]')\n", "ax.set_ylim([30, 75]) # Hz\n", "ax.set_xlabel('simulation time t [s]')\n", "None" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d29e8ee5bf7f4e94a1749239aec24d6d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(FloatProgress(value=0.0, max=221260.0), HTML(value='')))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "f_t = t\n", "\n", "n_f, n_t = Zxx.shape\n", "# start, stop = 180, 220\n", "# start, stop = 90, 110\n", "# start, stop = 15, 35\n", "# bounds_f = slice(start // 4 * nfft_factor, stop // 4 * nfft_factor)\n", "f_min, f_max = 30, 70 # Hz\n", "bounds_f = slice(np.argmax(f > f_min), np.argmin(f < f_max))\n", "\n", "\n", "f_mean = np.zeros(Zxx.shape[1])\n", "for le_t in tnrange(1, Zxx.shape[1] - 1):\n", " frame_f = f[bounds_f]\n", " frame_step = frame_f[1] - frame_f[0]\n", " time_step = f_t[1] - f_t[0]\n", " #if t == 10:\n", " # axs[-1].plot(frame_f, frame_Z)\n", " frame_Z = np.abs(Zxx[bounds_f, le_t])\n", " # frame_f = f[180:220]\n", " # frame_Z = np.abs(Zxx[180:220, 40])\n", " # frame_f = f[15:35]\n", " # frame_Z = np.abs(Zxx[15:35, 40])\n", " # plt.plot(frame_f, frame_Z)\n", "\n", " # peak_f = frame_f[np.argmax(frame)]\n", " # plt.axvline(peak_f, color='red')\n", "\n", "# def gauss(x, *p):\n", "# A, mu, sigma, o = p\n", "# return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + o\n", "\n", " def gauss(x, *p):\n", " A, mu, sigma = p\n", " return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n", "\n", " f_start = frame_f[np.argmax(frame_Z)]\n", " A_start = np.max(frame_Z)\n", " p0 = [A_start, f_start, 1.]\n", " try:\n", " coeff, var = optimize.curve_fit(gauss, frame_f, frame_Z, p0=p0)\n", " # plt.plot(frame_f, gauss(frame_f, *coeff))\n", " #print(coeff)\n", " A, mu, sigma, *_ = coeff\n", " f_mean[le_t] = mu\n", " except Exception:\n", " f_mean[le_t] = np.nan" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "47c40f28b5a34a94b4a631fd62107521", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(9, 5), sharex=True)\n", "fig.tight_layout(pad=2.2, h_pad=0, w_pad=1)\n", "\n", "label = f'Run {last_run}'\n", "ax.plot(f_t[1:-1], f_mean[1:-1])\n", "\n", "# b, a = signal.butter(3,\n", "# 1/5, # Hz\n", "# btype='lowpass',\n", "# fs=1/time_step)\n", "# filtered = signal.lfilter(b, a, f_mean[1:-1], axis=0)\n", "# ax.plot(f_t[1:-1], filtered)\n", "\n", "ax.set_title(label, pad=-20)\n", "ax.set_ylabel('f [Hz]')\n", "ax.grid()\n", "if not label in ['off_frequency', 'sweep_phase_steps']:\n", " ax.set_ylim([49.90, 50.10])\n", " var = np.var(f_mean[~np.isnan(f_mean)][1:-1])\n", " ax.text(0.5, 0.08, f'σ²={var * 1e3:.3g} mHz²', transform=ax.transAxes, ha='center', color='white', bbox=bbox)\n", " ax.text(0.5, 0.15, f'σ={np.sqrt(var) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center', color='white', bbox=bbox)\n", "\n", "# ax.text(0.5, 0.2, f'filt. σ²={np.var(filtered) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n", "else:\n", " f_min, f_max = min(f_mean[1:-1]), max(f_mean[1:-1])\n", " delta = f_max - f_min\n", " ax.set_ylim(f_min - delta * 0.1, f_max + delta * 0.3)\n", "\n", "for i in np.where(np.isnan(f_mean))[0]:\n", " ax.axvspan(f_t[i], f_t[i+1], color='lightblue')\n", "\n", "ax.set_xlabel('recording time t [s]')\n", "None" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "4a4cb62296df496bad37d93547d3c26a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f_copy = np.copy(f_mean[1:-1])\n", "f_copy[np.isnan(f_copy)] = np.mean(f_copy[~np.isnan(f_copy)])\n", "b, a = signal.cheby2(7, 86, 100, 'low', output='ba', fs=1000)\n", "filtered = signal.lfilter(b, a, f_copy)\n", "\n", "b2, a2 = signal.cheby2(3, 30, 1, 'high', output='ba', fs=1000)\n", "filtered2 = signal.lfilter(b2, a2, filtered)\n", "\n", "fig, (ax2, ax1) = plt.subplots(2, figsize=(9,7))\n", "ax1.plot(f_t[1:-1], f_copy, color='lightgray')\n", "ax1.set_ylim([49.90, 50.10])\n", "ax1.grid()\n", "formatter = matplotlib.ticker.FuncFormatter(lambda s, x: str(datetime.timedelta(seconds=s)))\n", "ax1.xaxis.set_major_formatter(formatter)\n", "zoom_offx = 7000 # s\n", "zoom_len = 300 # s\n", "ax1.set_xlim([zoom_offx, zoom_offx + zoom_len])\n", "\n", "ax1.plot(f_t[1:-1], filtered, color='orange')\n", "ax1r = ax1.twinx()\n", "ax1r.plot(f_t[1:-1], filtered2, color='red')\n", "ax1r.set_ylim([-0.015, 0.015])\n", "ax1.set_title(f'Zoomed trace ({datetime.timedelta(seconds=zoom_len)})', pad=-20)\n", "\n", "\n", "ax2.set_title(f'Run {last_run}')\n", "ax2.plot(f_t[1:-1], f_copy, color='orange')\n", "\n", "ax2r = ax2.twinx()\n", "ax2r.set_ylim([-0.1, 0.1])\n", "ax2r.plot(f_t[1:-1], filtered2, color='red')\n", "#ax2.plot(f_t[1:-1], filtered, color='orange', zorder=1)\n", "ax2.set_ylim([49.90, 50.10])\n", "ax2.set_xlim([0, f_t[-2]])\n", "ax2.grid()\n", "formatter = matplotlib.ticker.FuncFormatter(lambda s, x: str(datetime.timedelta(seconds=s)))\n", "ax2.xaxis.set_major_formatter(formatter)\n", "\n", "ax2.legend(handles=[\n", " patches.Patch(color='lightgray', label='Raw frequency'),\n", " patches.Patch(color='orange', label='low-pass filtered'),\n", " patches.Patch(color='red', label='band-pass filtered')])\n", "\n", "#ax2r.spines['right'].set_color('red')\n", "ax2r.yaxis.label.set_color('red')\n", "#ax2r.tick_params(axis='y', colors='red')\n", "\n", "#ax1r.spines['right'].set_color('red')\n", "ax1r.yaxis.label.set_color('red')\n", "#ax1r.tick_params(axis='y', colors='red')\n", "\n", "ax1.set_ylabel('f [Hz]')\n", "ax1r.set_ylabel('band-pass Δf [Hz]')\n", "ax2.set_ylabel('f [Hz]')\n", "ax2r.set_ylabel('band-pass Δf [Hz]')\n", "\n", "# Cut out first 10min of filtered data to give filters time to settle\n", "rms_slice = filtered2[np.where(f_t[1:] > 10*60)[0][0]:]\n", "rms = np.sqrt(np.mean(np.square(rms_slice)))\n", "ax1.text(0.5, 0.1, f'RMS (band-pass): {rms*1e3:.3f}mHz', transform=ax1.transAxes, color='white', bbox=bbox, ha='center')\n", "None" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "chunk_size = 256\n", "#\n", "#with open('filtered_freq.bin', 'wb') as f:\n", "# for chunk in range(0, len(rms_slice), chunk_size):\n", "# out_data = rms_slice[chunk:chunk+chunk_size]\n", "# f.write(struct.pack(f'{len(out_data)}f', *out_data))\n", "# \n", "#with open('raw_freq.bin', 'wb') as f:\n", "# for chunk in range(0, len(f_copy), chunk_size):\n", "# out_data = f_copy[chunk:chunk+chunk_size]\n", "# f.write(struct.pack(f'{len(out_data)}f', *out_data))" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":17: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n", " fig, ax = plt.subplots(figsize=(9,5))\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3809a1a83b5844e3906f1d74bcd15b5d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "text/plain": [ "5.0" ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = f_copy\n", "ys = scipy.fftpack.fft(data)\n", "ys = scipy.fftpack.fftshift(ys)\n", "#ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n", "#s = 3\n", "\n", "#ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n", "\n", "#xs = np.linspace(0, 5, len(data)//2)\n", "xs = np.linspace(-5, 5, len(data))\n", "\n", "#ys *= 2*np.pi*xs[s//2:-s//2+1]\n", "#ys *= xs\n", "\n", "#xs = np.linspace(len(data)/2, 1, len(data)/2)\n", "\n", "fig, ax = plt.subplots(figsize=(9,5))\n", "#ax.loglog(xs[s//2:-s//2+1], ys)\n", "#ax.loglog(xs[s//2:-s//2+1], ys)\n", "#ax.loglog(xs, ys)\n", "#ys[len(xs)//2] = 0\n", "#ax.set_yscale('log')\n", "ax.plot(xs, ys)\n", "#ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n", "ax.grid()\n", "#plt.show()\n", "xs[-1]" ] }, { "cell_type": "code", "execution_count": 156, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":20: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n", " fig, ax = plt.subplots()\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d137ae59ed7947ce8e3f7295e102f2f0", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "(1.6666666666666667e-05, 0.5)" ] }, "execution_count": 156, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of samplepoints\n", "N = len(data)\n", "# sample spacing\n", "T = 1.0 / 10.0\n", "x = np.linspace(0.0, N*T, N)\n", "yf = scipy.fftpack.fft(data)\n", "xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n", "\n", "yf = 2.0/N * np.abs(yf[:N//2])\n", "\n", "average_from = lambda val, start, average_width: np.hstack([val[:start], [ np.mean(val[i:i+average_width]) for i in range(start, len(val), average_width) ]])\n", "\n", "average_width = 6\n", "average_start = 20\n", "yf = average_from(yf, average_start, average_width)\n", "xf = average_from(xf, average_start, average_width)\n", "yf = average_from(yf, 200, average_width)\n", "xf = average_from(xf, 200, average_width)\n", "\n", "fig, ax = plt.subplots()\n", "ax.loglog(xf, yf)\n", "ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n", "ax.set_xlabel('T in s')\n", "ax.set_ylabel('Amplitude Δf')\n", "ax.grid()\n", "ax.set_xlim([1/60000, 0.5])" ] } ], "metadata": { "kernelspec": { "display_name": "labenv", "language": "python", "name": "labenv" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.1" } }, "nbformat": 4, "nbformat_minor": 4 }