770 lines
24 KiB
Text
770 lines
24 KiB
Text
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"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('data/waveform-raspi-2-2.sqlite3')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Run 000: 2020-03-25 16:07:36 - 2020-03-26 00:15:13 ( 8:07:37.266, 29261120sp)\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": "582c4360e293466e9baed5bc66a47883",
|
||
"version_major": 2,
|
||
"version_minor": 0
|
||
},
|
||
"text/plain": [
|
||
"HBox(children=(FloatProgress(value=0.0, max=914410.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)%0x10000:\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='<H')\n",
|
||
" write_index += 1\n",
|
||
" \n",
|
||
" elif seq > 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.0922848236977"
|
||
]
|
||
},
|
||
"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": "ecf3e3e261c54d87b169c1eb391a67f9",
|
||
"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": {
|
||
"text/plain": [
|
||
"(200.0, 800.0)"
|
||
]
|
||
},
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"window_len, window_len * nfft_factor"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "cc101709475d440ea77e68bcb56ce3b7",
|
||
"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": 12,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "b0202d59643548cf83b6fa6fd7580d2d",
|
||
"version_major": 2,
|
||
"version_minor": 0
|
||
},
|
||
"text/plain": [
|
||
"HBox(children=(FloatProgress(value=0.0, max=292611.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 as e:\n",
|
||
" f_mean[le_t] = np.nan"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "59d50a2876634433bdfb751a7a66d9ca",
|
||
"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",
|
||
"formatter = matplotlib.ticker.FuncFormatter(lambda s, x: str(datetime.timedelta(seconds=s)))\n",
|
||
"ax.xaxis.set_major_formatter(formatter)\n",
|
||
"ax.set_xlabel('recording time t [s]')\n",
|
||
"None"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 14,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "95cfea97c3b946ed9cdf351b16c2c45e",
|
||
"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": 15,
|
||
"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": 53,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"<ipython-input-53-c54c3e4ac2be>: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(figsize=(5, 2))\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "ad5454f664734681adb8640f480ce69a",
|
||
"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": [
|
||
"Text(20, 1, '50 Hz')"
|
||
]
|
||
},
|
||
"execution_count": 53,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"# Number of samplepoints\n",
|
||
"N = len(data)\n",
|
||
"# sample spacing\n",
|
||
"T = 1.0 / sampling_rate\n",
|
||
"x = np.linspace(0.0, N*T, N)\n",
|
||
"yf = scipy.fftpack.fft(data * sig.blackman(N))\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)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 68,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"<ipython-input-68-21b49a5af249>:1: 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=(6, 3))\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "4a3edb0925fe47eb8751150aa7da8c22",
|
||
"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": [
|
||
"The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"fig, ax = plt.subplots(figsize=(6, 3))\n",
|
||
"fig.tight_layout()\n",
|
||
"ax.loglog(xf, yf)\n",
|
||
"#ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: f'{1/x:.1f}'))\n",
|
||
"ax.set_xlabel('f [Hz]')\n",
|
||
"ax.set_ylabel('Amplitude V [V]')\n",
|
||
"ax.grid()\n",
|
||
"ax.set_xlim([0.1, 500])\n",
|
||
"fig.subplots_adjust(bottom=0.2)\n",
|
||
"\n",
|
||
"for le_f in (50, 150, 250, 350, 450):\n",
|
||
" ax.axvline(le_f, color=(1, 0.5, 0.5), zorder=-2)\n",
|
||
"ax.annotate('50 Hz', xy=(20, 1), xycoords='data', bbox=dict(fc='white', alpha=0.8, ec='none'))\n",
|
||
"font = {'family' : 'normal',\n",
|
||
" 'weight' : 'normal',\n",
|
||
" 'size' : 10}\n",
|
||
"matplotlib.rc('font', **font)\n",
|
||
"fig.savefig('fig_out/mains_voltage_spectrum.eps')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 43,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"<ipython-input-43-2e31f0cb9460>:21: 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": "37e159955f114f36824dd42da060b9ea",
|
||
"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": [
|
||
"[<matplotlib.lines.Line2D at 0x7fafc108cdc0>]"
|
||
]
|
||
},
|
||
"execution_count": 43,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"# Number of samplepoints\n",
|
||
"newcopy = np.copy(f_mean[1:-2])\n",
|
||
"N = len(newcopy)\n",
|
||
"# sample spacing\n",
|
||
"T = 1.0 / 10\n",
|
||
"x = np.linspace(0.0, N*T, N)\n",
|
||
"yf = scipy.fftpack.fft(newcopy * sig.blackman(N))\n",
|
||
"xf = np.linspace(0.0, 10/2, 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_width1, average_start1 = 3, 40\n",
|
||
"average_width2, average_start2 = 4, 100\n",
|
||
"yf = average_from(yf, average_start1, average_width1)\n",
|
||
"xf = average_from(xf, average_start1, average_width1)\n",
|
||
"yf = average_from(yf, average_start2, average_width2)\n",
|
||
"xf = average_from(xf, average_start2, average_width2)\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 [s]')\n",
|
||
"ax.set_ylabel('Amplitude Δf [Hz]')\n",
|
||
"\n",
|
||
"for i, t in enumerate([60, 300, 450, 1200, 1800]):\n",
|
||
" ax.axvline(1/t, color='red', alpha=0.5)\n",
|
||
" ax.annotate(f'{t} s', xy=(1/t, 3e-5), xytext=(-15, 0), xycoords='data', textcoords='offset pixels', rotation=90)\n",
|
||
"#ax.text(1/60, 10,'60 s', ha='left')\n",
|
||
"ax.grid()\n",
|
||
"#ax.set_xlim([1/60000, 0.5])\n",
|
||
"ax.set_ylim([5e-7, 2e-2])\n",
|
||
"ax.plot(xf[1:], 2e-6/xf[1:])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 18,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"application/vnd.jupyter.widget-view+json": {
|
||
"model_id": "671ae919bf124e72b54144310ea1602d",
|
||
"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()\n",
|
||
"ax.plot(np.linspace(0, (len(f_mean)-3)/10, len(f_mean)-3) , f_mean[1:-2])\n",
|
||
"ax.grid()"
|
||
]
|
||
}
|
||
],
|
||
"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.2"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 4
|
||
}
|