jaseg/master-thesis
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add some graphs, add frequency spectra comparison

Browse source 
compare between commercial measurements from Dr. Gobmaier GmbH and ours.
Turns out we agree!
This commit is contained in:
jaseg 2020-02-20 17:18:39 +00:00
parent 4eb0fab200
commit 809d6eeddc
3 changed files with 1214 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

225
lab-windows/dsss_experiments-ber.ipynb
View file

@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
@ -47,7 +47,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
@ -73,7 +73,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
@ -90,7 +90,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
@ -104,7 +104,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 40,
"metadata": {},
"outputs": [
{
@ -124,12 +124,12 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 77,
"metadata": {},
"outputs": [],
"source": [
"def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0):\n",
" test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration)\n",
"def generate_test_signal(duration, nbits=6, signal_amplitude=2.0e-3, decimation=10, seed=0, data=None):\n",
" test_data = np.random.RandomState(seed=seed).randint(0, 2 * (2**nbits), duration) if data is None else data\n",
" \n",
" signal = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n",
" noise = np.resize(mains_noise, len(signal))\n",
@ -139,7 +139,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
@ -153,13 +153,21 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-42-882fdfbdc9fa>:4: 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": "ac41374377284feaa0d8dc20bcb4a341",
"model_id": "cddadc4204f54789905875fcfb8e4126",
"version_major": 2,
"version_minor": 0
},
@ -177,7 +185,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
@ -186,7 +194,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
@ -366,13 +374,21 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 217,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-217-c761a61291a7>: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=(12, 9))\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b4952032e6ed4d62b23ddf7889a0efb9",
"model_id": "42da8a07c0e241f2b1d0f3358e68935b",
"version_major": 2,
"version_minor": 0
},
@ -382,19 +398,57 @@
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7efd98cd17f0>]"
]
},
"execution_count": 217,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(12, 9))\n",
"\n",
"params = dict(\n",
" nbits=8,\n",
" signal_amplitude=0.00015,\n",
" decimation=10,\n",
" threshold_factor=6.0,\n",
" power_avg_width=2.5,\n",
" max_lookahead=6.5)\n",
"decimation = 100\n",
"extra_dec = 10\n",
"nbits = 5\n",
"signal_amplitude=3e-3,\n",
"\n",
"seed = 42\n",
"\n",
"test_data, signal = generate_test_signal(duration, nbits, signal_amplitude, decimation, seed, data=np.array([4,5,6,7] * 8))\n",
"\n",
"#sosh = sig.butter(6, 1/(2**nbits), btype='highpass', output='sos', fs=decimation)\n",
"#sosl = sig.butter(6, 1.0, btype='lowpass', output='sos', fs=decimation)\n",
"#filtered = sig.sosfilt(sosh, sig.sosfilt(sosl, signal))\n",
"#filtered = sig.sosfilt(sosh, signal)\n",
"\n",
"cor_an1 = correlate(signal, nbits=nbits, decimation=decimation)\n",
"#cor_an2 = correlate(filtered, nbits=nbits, decimation=decimation)\n",
"#cor_an2 = correlate(sig.decimate(signal, 9), nbits=nbits, decimation=decimation//9)\n",
"cor_an2 = correlate(sig.decimate(signal, extra_dec), nbits=nbits, decimation=int(round(decimation/extra_dec)))\n",
"#ax.plot(cor_an[2])\n",
"#ax.matshow(sig.cwt(cor_an[2], sig.ricker, np.arange(1, 64)), aspect='auto')\n",
"cwt_ed1 = sig.cwt(cor_an1[2], sig.ricker, np.arange(1, 130))\n",
"cwt_ed2 = sig.cwt(cor_an2[2], sig.ricker, np.arange(1, 130))\n",
"#for f in [0.73, 1.0]:\n",
"# ax.plot(cwt_ed[int(round(f * decimation))], label=f'{f}')\n",
"\n",
"#ax.plot(signal)\n",
"#ax.twiny().plot(sig.decimate(signal, 9), color='orange')\n",
"#ax.twiny().plot(sig.decimate(signal, 10), color='orange')\n",
"\n",
"#ax.matshow(cwt_ed2, aspect='auto')\n",
"ax.plot(cwt_ed1[int(round(0.73 * decimation))], label=f'unfiltered')\n",
"ax.twinx().twiny().plot(cwt_ed2[int(round(0.73 * decimation / extra_dec))], label=f'filtered', color='orange')\n",
"\n",
"#ax.legend()\n",
"#ax.matshow(cor_an[2:4], aspect='auto')\n",
"#ser, br = run_ser_test(**params, sample_duration=100, print=print, seed=seed, ax=ax) #, debug_range=(16100, 16700))\n",
"#seed = 0xcbb3b8cf\n",
"#for seed in range(10):\n",
"# ser, br = run_ser_test(**params, sample_duration=32, print=print, seed=seed) #, debug_range=(16100, 16700))\n",
@ -612,21 +666,21 @@
},
{
"cell_type": "code",
"execution_count": 103,
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-103-3e101b35c97a>: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-36-8e813d331cd8>: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, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "84e1acbb31914f5f8b3a031654fd7733",
"model_id": "c05f590a7c0040628c5c5d4032bb7c8e",
"version_major": 2,
"version_minor": 0
},
@ -805,7 +859,7 @@
" std = data[:,1]\n",
" \n",
" ax.set_xlim([min(x), max(x)])\n",
" l = ax.plot(x, y, label='SER @ a=0.5', color='orange')\n",
" l = ax.plot(x, y, label='Amplitude at SER=0.5', color='orange')\n",
" \n",
" x, y, std = zip(*[ (le_x, le_y, le_std) for le_x, le_y, le_std in zip(x, y, std) if le_y is not None ])\n",
" y, std = np.array(y), np.array(std)\n",
@ -815,6 +869,7 @@
" ax.fill_between([-1, min(ser_valid)], 0, 1, facecolor='red', alpha=0.2, transform=trans, zorder=1)\n",
" ax.fill_between([max(ser_valid), max(ser_valid)*10], 0, 1, facecolor='red', alpha=0.2, transform=trans)\n",
" ax.set_ylim([min(y)*0.9, max(y)*1.1])\n",
" ax.grid()\n",
" return l\n",
"\n",
"l1 = plot_intercepts(intercept_ax)\n",
@ -835,6 +890,124 @@
"source": [
"#fig.savefig('dsss_prototype_symbol_error_rate_5-8_bit.svg')"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-35-cafaa6062c72>: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, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fa93fb31355840148f941dab8b60f51c",
"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": [
"<ipython-input-35-cafaa6062c72>:16: RuntimeWarning: divide by zero encountered in log10\n",
" cm_func = lambda x: cmap(np.log10(x - min(decimations)) / (np.log10(max(decimations)) - np.log10(min(decimations))))\n"
]
}
],
"source": [
"fig, ((ax, cbar_ax), (intercept_ax, empty)) = plt.subplots(2, 2, figsize=(12, 9), gridspec_kw={'width_ratios': [1, 0.05]})\n",
"empty.axis('off')\n",
"#fig.tight_layout()\n",
"\n",
"results = []\n",
"for fn in [\n",
" '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-14-10-13.json',\n",
" '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-13-21-57.json',\n",
" '/mnt/c/Users/jaseg/shared/dsss_experiments_res-2020-02-20-13-23-47.json',\n",
"]:\n",
" with open(fn, 'r') as f:\n",
" results += json.load(f)\n",
"\n",
"decimations = [decimation for (_nbits, thf, _reps, _points, _duration, decimation), series in results]\n",
"cmap = matplotlib.cm.viridis\n",
"cm_func = lambda x: cmap(np.log10(x - min(decimations)) / (np.log10(max(decimations)) - np.log10(min(decimations))))\n",
"\n",
"decimation_sers = {}\n",
"for (nbits, thf, reps, points, duration, decimation), series in results:\n",
" data = [ [ mean for mean, _std, _msg in reps if mean is not None ] for _amp, reps in series ]\n",
" amps = [ amp for amp, _reps in series ]\n",
" sers = np.array([ np.mean(values) for values in data ])\n",
" stds = np.array([ np.std(values) for values in data ])\n",
" decimation_sers[decimation] = list(zip(amps, sers, stds))\n",
" \n",
" l, = ax.plot(amps, np.clip(sers, 0, 1), label=f'decimation={decimation}', color=cm_func(decimation))\n",
" ax.fill_between(amps, np.clip(sers + stds, 0, 1), np.clip(sers - stds, 0, 1), facecolor=l.get_color(), alpha=0.2)\n",
" ax.axhline(0.5, color='gray', ls=(0, (3, 4)), lw=0.8)\n",
"ax.grid()\n",
"ax.set_xlabel('Amplitude in mHz')\n",
"ax.set_ylabel('Symbol error rate')\n",
"\n",
"norm = matplotlib.colors.Normalize(vmin=np.log10(min(decimations)), vmax=np.log10(max(decimations)))\n",
"cb1 = matplotlib.colorbar.ColorbarBase(cbar_ax, cmap=cmap, norm=norm, orientation='vertical', label=\"Decimation\", ticks=[np.log10(d) for d in decimations])\n",
"cb1.ax.set_yticklabels([f'{d:.1f}' for d in decimations])\n",
"\n",
"def plot_intercepts(ax, SER_TH = 0.5):\n",
" intercepts = {}\n",
" for dec, sers in decimation_sers.items():\n",
" last_ser, last_amp, last_std = 0, 0, 0\n",
" for amp, ser, std in sorted(sers):\n",
" if last_ser > SER_TH and ser < SER_TH:\n",
" icp = last_amp + (SER_TH - last_ser) / (ser - last_ser) * (amp - last_amp)\n",
" ic_std = abs(last_amp - amp) / 2# np.sqrt(np.mean(last_std**2 + std**2))\n",
" intercepts[dec] = (icp, ic_std)\n",
" break\n",
" last_amp, last_ser = amp, ser\n",
" else:\n",
" intercepts[dec] = None, None\n",
" \n",
" ser_valid = [dec for dec, (ser, _std) in intercepts.items() if ser is not None]\n",
" #ax.axvline(min(ser_valid), color='red')\n",
" #ax.axvline(max(ser_valid), color='red')\n",
" \n",
" x = sorted(intercepts.keys())\n",
" data = np.array([ intercepts[dec] for dec in x ])\n",
" y = data[:,0]\n",
" std = data[:,1]\n",
" \n",
" ax.set_xlim([min(x), max(x)])\n",
" l = ax.plot(x, y, label='Amplitude at SER=0.5', color='orange')\n",
" ax.legend(loc=3)\n",
" ax.grid()\n",
" \n",
" x, y, std = zip(*[ (le_x, le_y, le_std) for le_x, le_y, le_std in zip(x, y, std) if le_y is not None ])\n",
" y, std = np.array(y), np.array(std)\n",
" ax.fill_between(x, y-std, y+std, color=l[0].get_color(), alpha=0.3)\n",
" \n",
" trans = matplotlib.transforms.blended_transform_factory(ax.transData, ax.transAxes)\n",
" ax.fill_between([-1, min(ser_valid)], 0, 1, facecolor='red', alpha=0.2, transform=trans, zorder=1)\n",
" ax.fill_between([max(ser_valid), max(ser_valid)*10], 0, 1, facecolor='red', alpha=0.2, transform=trans)\n",
" ax.set_ylim([min(y)*0.9, max(y)*1.1])\n",
" ax.set_xscale('log')\n",
" ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, _: '{:g}'.format(x)))\n",
" ax.set_xticks([1, 2, 5, 10, 20, 50])\n",
" ax.set_xlim([1, 60])\n",
" return l\n",
"\n",
"l1 = plot_intercepts(intercept_ax)\n"
]
}
],
"metadata": {

340
lab-windows/grid_frequency_spectra.ipynb Normal file
View file

@ -0,0 +1,340 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import csv\n",
"\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"import scipy.fftpack"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data = np.genfromtxt('data/Netzfrequenz_Sekundenwerte_2012_KW37.csv', delimiter=',')[1:,1:]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "36f1f4d7970e41afa4737b6f63b7c449",
"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 0x7f952a6e4580>]"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(data[:3600*24, 0])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.02051102806199375"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.std(data[:,0])"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "d7fe0512f4254efeb15235a5617ef064",
"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": [
"(1e-06, 0.5)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Number of samplepoints\n",
"N = len(data[:,0])\n",
"# sample spacing\n",
"T = 1.0\n",
"x = np.linspace(0.0, N*T, N)\n",
"yf = scipy.fftpack.fft(data[:,0])\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",
"#yf = sum(yf[s::10] for s in range(10)) / 10\n",
"#xf = sum(xf[s::10] for s in range(10)) / 10\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/1000000, 0.5])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fd94bc97a276400db0539b703c4eeeac",
"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": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Number of samplepoints\n",
"N = len(data[:,0])\n",
"# sample spacing\n",
"T = 1.0\n",
"x = np.linspace(0.0, N*T, N)\n",
"yf = scipy.fftpack.fft(data[:,0])\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 = 20\n",
"average_start = 100\n",
"yf = average_from(yf, average_start, average_width)\n",
"xf = average_from(xf, average_start, average_width)\n",
"yf = average_from(yf, 300, average_width)\n",
"xf = average_from(xf, 300, 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])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "object of type <class 'float'> cannot be safely interpreted as an integer.",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m~/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/function_base.py\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(start, stop, num, endpoint, retstep, dtype, axis)\u001b[0m\n\u001b[1;32m 116\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 117\u001b[0;31m \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moperator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\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 118\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mTypeError\u001b[0m: 'float' object cannot be interpreted as an integer",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-7-75728c9461c4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mys\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconvolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'valid'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\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 8\u001b[0m \u001b[0;31m#xs = np.linspace(len(data)/2, 1, len(data)/2)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
"\u001b[0;32m~/safety-reset/lab-windows/env/lib/python3.8/site-packages/numpy/core/function_base.py\u001b[0m in \u001b[0;36mlinspace\u001b[0;34m(start, stop, num, endpoint, retstep, dtype, axis)\u001b[0m\n\u001b[1;32m 117\u001b[0m \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moperator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 118\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 119\u001b[0;31m raise TypeError(\n\u001b[0m\u001b[1;32m 120\u001b[0m \u001b[0;34m\"object of type {} cannot be safely interpreted as an integer.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 121\u001b[0m .format(type(num)))\n",
"\u001b[0;31mTypeError\u001b[0m: object of type <class 'float'> cannot be safely interpreted as an integer."
]
}
],
"source": [
"ys = scipy.fftpack.fft(data[:,0])\n",
"ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
"s = 60\n",
"\n",
"ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
"\n",
"xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
"#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.loglog(xs[s//2:-s//2+1], ys)\n",
"ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ys = scipy.fftpack.fft(data[:,0])\n",
"ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
"s = 1\n",
"\n",
"ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
"\n",
"xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
"#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.loglog(xs[s//2:-s//2+1 if s > 1 else None], ys)\n",
"ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ys = scipy.fftpack.fft(data[:,0])\n",
"ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
"s = 1\n",
"\n",
"ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
"\n",
"xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
"\n",
"ys *= 2*np.pi*xs\n",
"#xs = np.linspace(len(data)/2, 1, len(data)/2)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.loglog(xs[s//2:-s//2+1 if s > 1 else None], ys)\n",
"ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ys = scipy.fftpack.fft(data[:,0])\n",
"ys = 2.0/len(data) * np.abs(ys[:len(data)//2])\n",
"s = 30\n",
"\n",
"ys = np.convolve(ys, np.ones((s,))/s, mode='valid')\n",
"\n",
"xs = np.linspace(0, 1.0/2.0, len(data)/2)\n",
"\n",
"ys *= 2*np.pi*xs[s//2:-s//2+1]\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.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _pos: 1/x))\n",
"ax.grid()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"1/0.0628"
]
}
],
"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
}

675
lab-windows/grid_scope.ipynb Normal file
View file

@ -0,0 +1,675 @@
{
"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='<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.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": [
"<ipython-input-81-a751e13723ea>: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": [
"<ipython-input-156-ddde6af5dee1>: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
}