thesis: Add freq measurement notebook

This commit is contained in:
jaseg 2020-04-02 15:11:17 +02:00
parent 7f0041bf10
commit e104e5d2dc
9 changed files with 4650 additions and 9140 deletions

1
.gitignore vendored
View file

@ -13,3 +13,4 @@ gm_platform/platform/gerber
*.run.xml
*.toc
*.bib~
*-eps-converted-to.pdf

View file

@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
@ -18,7 +18,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
@ -27,7 +27,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
@ -63,7 +63,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
@ -76,7 +76,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
@ -99,7 +99,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 17,
"metadata": {},
"outputs": [
{
@ -108,7 +108,7 @@
"(129, 470)"
]
},
"execution_count": 6,
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
@ -119,7 +119,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 18,
"metadata": {},
"outputs": [
{
@ -128,7 +128,7 @@
"3.90625"
]
},
"execution_count": 7,
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
@ -139,13 +139,13 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "c7a38939cb3a42cfbfe99fe7d603e11a",
"model_id": "96e3c06718434807813d909d3a1f4d0e",
"version_major": 2,
"version_minor": 0
},
@ -174,13 +174,13 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "e6b769f15870407ba59cff627632d447",
"model_id": "af9cb4a349eb42c6b89d7951dabdd71e",
"version_major": 2,
"version_minor": 0
},
@ -209,7 +209,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 21,
"metadata": {},
"outputs": [
{
@ -239,7 +239,7 @@
" 492.1875 , 496.09375, 500. ])"
]
},
"execution_count": 10,
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
@ -250,13 +250,13 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "db4021b95410428cb07f76edd932907c",
"model_id": "930f346b46d64489bdc790ec11ef9a69",
"version_major": 2,
"version_minor": 0
},
@ -266,6 +266,36 @@
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/user/safety-reset/lab-windows/env/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
" warnings.warn('Covariance of the parameters could not be estimated',\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n",
"slice(12, 15, None)\n"
]
}
],
"source": [

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,153 @@
{
"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",
"from scipy import optimize as opt\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_1pps_debug.sqlite3')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Run 000: 2020-03-31 16:58:00 - 2020-03-31 16:58:36 ( 0:00:36.029, 36512sp)\n",
"Run 001: 2020-03-31 16:58:51 - 2020-03-31 17:05:19 ( 0:06:27.729, 392608sp)\n",
"Run 002: 2020-03-31 17:07:02 - 2020-03-31 17:41:34 ( 0:34:32.105, 37024sp)\n",
"Run 003: 2020-03-31 18:50:05 - 2020-03-31 18:50:43 ( 0:00:37.576, 38048sp)\n",
"Run 004: 2020-03-31 18:54:08 - 2020-03-31 19:14:32 ( 0:20:24.104, 1239424sp)\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": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "2712791d61b549ecb531f886d88e1d54",
"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": [
"histogram = np.array(db.execute('SELECT gps_1pps, COUNT(*) FROM measurements WHERE gps_1pps != -1 AND run_id = ? GROUP BY gps_1pps', (last_run,)).fetchall())\n",
"hist_plot = histogram.astype(float)[1:-1]\n",
"hist_plot[:, 0] *= 2 / 5 * 2\n",
"hist_plot[:, 1] /= (1000 / 32)\n",
"\n",
"f_nom = 19.440e6\n",
"\n",
"font = {'family' : 'normal',\n",
" 'weight' : 'normal',\n",
" 'size' : 10}\n",
"matplotlib.rc('font', **font)\n",
"fig, ax = plt.subplots(figsize=(5, 4))\n",
"ax.grid()\n",
"# We have a bug that causes our measurements to occassionally be out by +/- 65534 counts. For now, fix this by simply throwing away these (very obviously invalid) bins.\n",
"ax.bar(hist_plot[:,0] - f_nom , hist_plot[:, 1])\n",
"\n",
"def gauss(x, *p):\n",
" A, mu, sigma = p\n",
" return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n",
"\n",
"gauss_x = np.linspace(np.min(hist_plot[:,0]), np.max(hist_plot[:,0]), 10000)\n",
"coeff, var_matrix = opt.curve_fit(gauss, hist_plot[:,0], hist_plot[:,1], p0=[np.max(hist_plot[:,1]), np.mean(hist_plot[:,0]), 1])\n",
"hist_fit = gauss(gauss_x, *coeff)\n",
"ax.plot(gauss_x - f_nom, hist_fit, color='orange')\n",
"_A, mu, sigma = coeff\n",
"bbox_props = dict(fc='white', alpha=0.8, ec='none')\n",
"ax.annotate(f'σ² = {sigma**2 * 1e3:.1f} mHz ({sigma**2 / f_nom * 1e9:.2f} ppb)\\nμ = {mu-f_nom:+.1f} Hz ({(mu-f_nom)/f_nom * 1e6:+.2f} ppm)', xy=[0.6, 0.5], xycoords='figure fraction', bbox=bbox_props)\n",
"ax.set_xlabel('$f - f_{nom}$ [Hz]')\n",
"ax.set_ylabel('# observations')\n",
"\n",
"#ax.set_title('OCXO frequency derivation relative to GPS 1pps')\n",
"fig.savefig('fig_out/ocxo_freq_stability.eps', format='eps')"
]
}
],
"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
}

File diff suppressed because one or more lines are too long

View file

@ -1,5 +1,12 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Grid frequency estimation"
]
},
{
"cell_type": "code",
"execution_count": 1,
@ -30,13 +37,23 @@
"%matplotlib widget"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 1: Setup\n",
"### Load data series information from sqlite capture file\n",
"\n",
"One capture file may contain multiple runs/data series. Display a list of runs and their start/end time and sample count, then select the newest one in `last_run` variable."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"db = sqlite3.connect('data/waveform-raspi-2-2.sqlite3')"
"db = sqlite3.connect('data/waveform-raspi-ocxo-2.sqlite3')"
]
},
{
@ -48,7 +65,9 @@
"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"
"Run 000: 2020-04-01 14:00:25 - 2020-04-01 15:09:31 ( 1:09:05.846, 4197664sp)\n",
"Run 001: 2020-04-02 11:56:41 - 2020-04-02 11:57:59 ( 0:01:18.544, 79552sp)\n",
"Run 002: 2020-04-02 12:03:51 - 2020-04-02 12:12:54 ( 0:09:03.033, 549792sp)\n"
]
}
],
@ -57,8 +76,16 @@
" 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"
"last_run, n_records = run_id, count"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setup analog parameters\n",
"\n",
"Setup parameters of analog capture hardware here. This is used to scale samples from ADC counts to analog voltages. Also setup sampling rate here. Nominal sampling rate is 1 ksps."
]
},
{
@ -67,9 +94,11 @@
"metadata": {},
"outputs": [],
"source": [
"sampling_rate = 1000.0 * 48.6 / 48\n",
"\n",
"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",
"# Note: These are for the first prototype only!\n",
"vmeas_source_impedance = 330e3\n",
"vmeas_source_scale = 0.5\n",
"\n",
@ -92,20 +121,29 @@
"adc_count_to_vmeas = lambda x: (x*adc_val_to_voltage_factor - v0) / (v100-v0) * 100"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load run data from sqlite3 capture file\n",
"\n",
"Load measurement data for the selected run and assemble a numpy array containing one continuous trace. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "582c4360e293466e9baed5bc66a47883",
"model_id": "38970a137bed494d8c28c270471c73df",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(FloatProgress(value=0.0, max=914410.0), HTML(value='')))"
"HBox(children=(FloatProgress(value=0.0, max=17181.0), HTML(value='')))"
]
},
"metadata": {},
@ -115,7 +153,8 @@
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
"\n",
"RMS voltage: 228.5564548966498\n"
]
}
],
@ -145,39 +184,44 @@
" # nans[:] = np.nan\n",
" # data = np.append(data, nans) FIXME\n",
" \n",
"data = (data * adc_val_to_voltage_factor - v0) / (v100-v0) * 100"
"data = (data * adc_val_to_voltage_factor - v0) / (v100-v0) * 100\n",
"\n",
"# https://stackoverflow.com/questions/6518811/interpolate-nan-values-in-a-numpy-array\n",
"nan_helper = lambda y: (np.isnan(y), lambda z: z.nonzero()[0])\n",
"\n",
"# data rarely may contain NaNs where the capture script failed to read and acknowledge capture buffers from the sensor board fast enough.\n",
"# For RMS calculation and overall FFT fill these NaNs with interpolated values from their neighbors.\n",
"data_not_nan = np.copy(data)\n",
"nans, x = nan_helper(data)\n",
"data_not_nan[nans]= np.interp(x(nans), x(~nans), data[~nans])\n",
"\n",
"print('RMS voltage:', np.sqrt(np.mean(np.square(data_not_nan))))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"cell_type": "markdown",
"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)))"
"### Show a preview of loaded data"
]
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-33-4419e570bd12>: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, (top, bottom) = plt.subplots(2, figsize=(9,6))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "ecf3e3e261c54d87b169c1eb391a67f9",
"model_id": "6faa0da0b6b64d0094b9b683e5f2c434",
"version_major": 2,
"version_minor": 0
},
@ -191,7 +235,7 @@
],
"source": [
"fig, (top, bottom) = plt.subplots(2, figsize=(9,6))\n",
"fig.tight_layout(pad=3, h_pad=0.1)\n",
"fig.tight_layout(pad=3, h_pad=1.8)\n",
"\n",
"range_start, range_len = -300, 60 # [s]\n",
"\n",
@ -205,68 +249,85 @@
"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.02, 0.5, f'mean: {mean:.3f}', transform=top.transAxes, color='white', bbox=bbox, ha='left', va='center')\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",
"top.text(0.5, 0.9, f'Run {run_id}', transform=top.transAxes, color='white', bbox=bbox, ha='center', fontweight='bold')\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",
"\n",
"top.set_title('Voltage waveform')\n",
"bottom.set_title('Voltage frequency spectrum')\n",
"None"
]
},
{
"cell_type": "code",
"execution_count": 9,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"## Step 2: Calculate Short-Time Fourier Transform of capture"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Window length: 405 sp, zero-padded to 405 sp\n"
]
}
],
"source": [
"fs = sampling_rate # Hz\n",
"ff = 50 # Hz\n",
"\n",
"analysis_periods = 10\n",
"analysis_periods = 20\n",
"window_len = fs * analysis_periods/ff\n",
"nfft_factor = 4\n",
"nfft_factor = 1\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)"
" nfft = window_len * nfft_factor)\n",
"print(f'Window length: {window_len:.0f} sp, zero-padded to {window_len * nfft_factor:.0f} sp')"
]
},
{
"cell_type": "code",
"execution_count": 10,
"cell_type": "markdown",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(200.0, 800.0)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"window_len, window_len * nfft_factor"
"### Show a preview of STFT results\n",
"\n",
"Cut out our approximate frequency range of interest"
]
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 64,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-64-467ca72791b1>: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=(9, 3))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "cc101709475d440ea77e68bcb56ce3b7",
"model_id": "95c94ad504b248febb55f81b0e919464",
"version_major": 2,
"version_minor": 0
},
@ -291,20 +352,34 @@
"None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 3: Run Gasior and Gonzalez for precise frequency estimation\n",
"\n",
"Limit analysis to frequency range of interest. If automatic adaption to totally different frequency ranges\n",
"(e.g. 400Hz) would be necessary, we could switch here based on configuration or a lookup of the STFT bin\n",
"containing highest overall energy.\n",
"\n",
"As elaborated in the Gasior and Gonzalez Paper [1] the shape of the template function should match the expected peak shape.\n",
"Peak shape is determined by the STFT window function. As Gasior and Gonzalez note, a gaussian is a very good fit for a steep gaussian window."
]
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b0202d59643548cf83b6fa6fd7580d2d",
"model_id": "5443a7d157cc4a828cbffc91b2d645e3",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(FloatProgress(value=0.0, max=292611.0), HTML(value='')))"
"HBox(children=(FloatProgress(value=0.0, max=2708.0), HTML(value='')))"
]
},
"metadata": {},
@ -322,61 +397,65 @@
"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",
"# Frequency ROI\n",
"f_min, f_max = 30, 70 # Hz\n",
"# Indices of bins within ROI\n",
"bounds_f = slice(np.argmax(f > f_min), np.argmin(f < f_max))\n",
"\n",
"\n",
"# Initialize output array\n",
"f_mean = np.zeros(Zxx.shape[1])\n",
"\n",
"# Iterate over STFT time slices\n",
"for le_t in tnrange(1, Zxx.shape[1] - 1):\n",
" # Cut out ROI and compute magnitude of complex fourier coefficients\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",
" # Template function. We use a gaussian here. This function needs to fit the window above.\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",
" # Calculate initial values for curve fitting\n",
" f_start = frame_f[np.argmax(frame_Z)] # index of strongest bin index\n",
" A_start = np.max(frame_Z) # strongest bin value\n",
" p0 = [A_start, f_start, 1.]\n",
" try:\n",
" # Fit template to measurement data STFT ROI \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",
" _A, f_mean[le_t], _sigma, *_ = coeff # The measured frequency is the mean of the fitted gaussian\n",
" \n",
" except Exception as e:\n",
" # Handle fit errors\n",
" f_mean[le_t] = np.nan"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Produce a plot of measurement results\n",
"\n",
"Include measurements of mean, standard deviation and variance of measurement data"
]
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 56,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-56-4e12d8913585>: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=(9, 5), sharex=True)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "59d50a2876634433bdfb751a7a66d9ca",
"model_id": "c1c391d0577f4dbeb59db8d1fa9261de",
"version_major": 2,
"version_minor": 0
},
@ -392,49 +471,45 @@
"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",
"# Cut off invalid values at fringes\n",
"ax.plot(f_t[1:-2], f_mean[1:-2])\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",
"var = np.var(f_mean[~np.isnan(f_mean)][1:-1])\n",
"mean = np.mean(f_mean[~np.isnan(f_mean)][1:-1])\n",
"ax.text(0.5, 0.95, f'Run {run_id}', transform=ax.transAxes, ha='center', fontweight='bold', color='white', bbox=bbox)\n",
"ax.text(0.05, 0.95, f'μ={mean:.3g} Hz', transform=ax.transAxes, ha='left', color='white', bbox=bbox)\n",
"ax.text(0.05, 0.89, f'σ={np.sqrt(var) * 1e3:.3g} mHz', transform=ax.transAxes, ha='left', color='white', bbox=bbox)\n",
"ax.text(0.05, 0.83, f'σ²={var * 1e3:.3g} mHz²', transform=ax.transAxes, ha='left', color='white', bbox=bbox)\n",
"\n",
"# Indicated missing values\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",
"ax.set_xlabel('recording time t [hh:mm:ss]')\n",
"None"
]
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 57,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-57-8b77e38496af>:9: 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, (ax2, ax1) = plt.subplots(2, figsize=(9,7))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "95cfea97c3b946ed9cdf351b16c2c45e",
"model_id": "c1209b4895814c92a9b0fa01ad666667",
"version_major": 2,
"version_minor": 0
},
@ -444,6 +519,17 @@
},
"metadata": {},
"output_type": "display_data"
},
{
"ename": "IndexError",
"evalue": "index 0 is out of bounds for axis 0 with size 0",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-57-8b77e38496af>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 56\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 57\u001b[0m \u001b[0;31m# Cut out first 10min of filtered data to give filters time to settle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 58\u001b[0;31m \u001b[0mrms_slice\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfiltered2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwhere\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf_t\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m60\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 59\u001b[0m \u001b[0mrms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msquare\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrms_slice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 60\u001b[0m \u001b[0max1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mf'RMS (band-pass): {rms*1e3:.3f}mHz'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtransform\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransAxes\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'white'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbbox\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbbox\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mha\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'center'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mIndexError\u001b[0m: index 0 is out of bounds for axis 0 with size 0"
]
}
],
"source": [
@ -512,7 +598,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
@ -531,49 +617,16 @@
},
{
"cell_type": "code",
"execution_count": 53,
"execution_count": 73,
"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"
}
],
"outputs": [],
"source": [
"# Number of samplepoints\n",
"N = len(data)\n",
"N = len(data_not_nan)\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",
"yf = scipy.fftpack.fft(data_not_nan * 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",
@ -590,21 +643,21 @@
},
{
"cell_type": "code",
"execution_count": 68,
"execution_count": 77,
"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",
"<ipython-input-77-2f4bcf6b2d33>: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",
"model_id": "8b689c4f96fa40ffb5012764afb57564",
"version_major": 2,
"version_minor": 0
},
@ -631,7 +684,7 @@
"ax.set_xlabel('f [Hz]')\n",
"ax.set_ylabel('Amplitude V [V]')\n",
"ax.grid()\n",
"ax.set_xlim([0.1, 500])\n",
"ax.set_xlim([0.001, 500])\n",
"fig.subplots_adjust(bottom=0.2)\n",
"\n",
"for le_f in (50, 150, 250, 350, 450):\n",
@ -646,21 +699,21 @@
},
{
"cell_type": "code",
"execution_count": 43,
"execution_count": 84,
"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",
"<ipython-input-84-936ca777d145>:26: 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",
"model_id": "b58b8858dea1485fae236c9fbb6954d5",
"version_major": 2,
"version_minor": 0
},
@ -670,21 +723,15 @@
},
"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",
"nans, x = nan_helper(newcopy)\n",
"newcopy[nans]= np.interp(x(nans), x(~nans), newcopy[~nans])\n",
"\n",
"N = len(newcopy)\n",
"# sample spacing\n",
"T = 1.0 / 10\n",
@ -715,8 +762,8 @@
"#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:])"
"#ax.set_ylim([5e-7, 2e-2])\n",
"#ax.plot(xf[1:], 2e-6/xf[1:])"
]
},
{
@ -744,6 +791,15 @@
"ax.plot(np.linspace(0, (len(f_mean)-3)/10, len(f_mean)-3) , f_mean[1:-2])\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"1. **Gasior, M. & Gonzalez, J.** Improving FFT frequency measurement resolution by parabolic and gaussian interpolation *CERN-AB-Note-2004-021, CERN-AB-Note-2004-021, 2004*"
]
}
],
"metadata": {

View file

@ -1,4 +1,6 @@
LAB_PATH ?= ../lab-windows
SHELL := bash
.ONESHELL:
.SHELLFLAGS := -eu -o pipefail -c
@ -8,11 +10,16 @@ MAKEFLAGS += --no-builtin-rules
all: safety_reset.pdf
safety_reset.pdf: resources/grid_freq_estimation.pdf
%.pdf: %.tex %.bib
pdflatex -shell-escape $<
biber $*
pdflatex -shell-escape $<
resources/%.pdf: $(LAB_PATH)/%.ipynb
jupyter-nbconvert --to=pdf --output-dir=resources --output=$* --LatexExporter.template_file=resources/nbexport.tplx $^
.PHONY: clean
clean:
rm -f safety_reset.aux safety_reset.bbl safety_reset.bcf safety_reset.log safety_reset.blg

View file

@ -0,0 +1,13 @@
((*- extends 'article.tplx' -*))
((*- block header -*))
((( super() )))
\pagenumbering{gobble}
((*- endblock header -*))
((* block maketitle *))((* endblock maketitle *))

View file

@ -47,6 +47,7 @@
\usetikzlibrary{calc}
%\usepackage[pdftex]{graphicx,color}
\usepackage{epstopdf}
\usepackage{pdfpages}
% Needed for murks.tex
\usepackage{setspace}
\usepackage[draft=false,babel,tracking=true,kerning=true,spacing=true]{microtype} % optischer Randausgleich etc.
@ -56,6 +57,10 @@
\newcommand{\degree}{\ensuremath{^\circ}}
\newcolumntype{P}[1]{>{\centering\arraybackslash}p{#1}}
\usepackage{fancyhdr}
\fancyhf{}
\fancyfoot[C]{\thepage}
\begin{document}
% Beispielhafte Nutzung der Vorlage für die Titelseite (bitte anpassen):
@ -778,7 +783,7 @@ an \emph{arbitrary} signal this would highly limit our practical measurement acc
instead just amounts to an interpolation between output bins. Depending on the downstream analysis algorithm it may
still be sensible to use this property of the DFT for interpolation, but in general it will be computationally
expensive compared to other interpolation methods and in any case it will not yield any better frequency resolution
aside from a hypothetical numerical advantage\cite{gasior01}.
aside from a hypothetical numerical advantage\cite{gasior02}.
}.
For this reason all approaches to mains frequency estimation are based on a model of the mains voltage waveform.
Nominally, this waveform would be a perfect sine at $f = 50 \text{Hz}$. In practice it is a sine at $f \approx 50
@ -989,7 +994,8 @@ interface and its good tolerance of system resets due to unexpected power loss.
\subsection{Frequency sensor measurement results}
Captured raw waveform data is processed in the Jupyter Lab environment\cite{kluyver01} and grid frequency estimates are
extracted as described in sec. \ref{frequency_estimation} using the \textcite{gasior01} technique.
extracted as described in sec. \ref{frequency_estimation} using the \textcite{gasior01} technique. Appendix
\ref{grid_freq_est_notebook} contains the Jupyter notebook we used for frequency measurement.
% FIXME comparison against reference measurements?
@ -1088,6 +1094,13 @@ correctly configure than it is to simply use separate hardware and secure the in
\printbibliography
\newpage
\chapter{Transcripts of Jupyter notebooks used in this thesis}
\section{Grid frequency estimation}
\label{grid_freq_est_notebook}
\fancyhead[C]{Included Jupyter notebook: Grid frequency estimation}
\includepdf[pages=-, pagecommand=\thispagestyle{fancy}]{resources/grid_freq_estimation.pdf}
\chapter{Demonstrator schematics and code}
\chapter{Economic viability of countermeasures}