diff options
Diffstat (limited to 'lab-windows/grid_scope.ipynb')
| -rw-r--r-- | lab-windows/grid_scope.ipynb | 388 |
1 files changed, 222 insertions, 166 deletions
diff --git a/lab-windows/grid_scope.ipynb b/lab-windows/grid_scope.ipynb index 1ddde37..0cb2083 100644 --- a/lab-windows/grid_scope.ipynb +++ b/lab-windows/grid_scope.ipynb @@ -1,6 +1,13 @@ { "cells": [ { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Grid frequency estimation" + ] + }, + { "cell_type": "code", "execution_count": 1, "metadata": {}, @@ -31,12 +38,22 @@ ] }, { + "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", @@ -93,19 +122,28 @@ ] }, { + "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": "markdown", + "metadata": {}, + "source": [ + "## Step 2: Calculate Short-Time Fourier Transform of capture" + ] + }, + { "cell_type": "code", - "execution_count": 9, + "execution_count": 38, "metadata": {}, - "outputs": [], + "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 }, @@ -292,19 +353,33 @@ ] }, { + "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", "\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": { |