diff options
| -rw-r--r-- | lab-windows/dsss_experiments.ipynb | 320 |
1 files changed, 282 insertions, 38 deletions
diff --git a/lab-windows/dsss_experiments.ipynb b/lab-windows/dsss_experiments.ipynb index bb56362..9566edf 100644 --- a/lab-windows/dsss_experiments.ipynb +++ b/lab-windows/dsss_experiments.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 10, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -63,13 +63,13 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "7560730a2391425ab9dad7a1f22e5fb2", + "model_id": "2450440508db4069b3df05fa08346b39", "version_major": 2, "version_minor": 0 }, @@ -90,10 +90,10 @@ { "data": { "text/plain": [ - "<matplotlib.image.AxesImage at 0x7f0496c50b80>" + "<matplotlib.image.AxesImage at 0x7fe4abaf7490>" ] }, - "execution_count": 5, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -105,13 +105,14 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def modulate(data, nbits=5, code=29):\n", " # 0, 1 -> -1, 1\n", " mask = gold(nbits)[code]*2 - 1\n", + " \n", " # same here\n", " data_centered = (data*2 - 1)\n", " return (mask[:, np.newaxis] @ data_centered[np.newaxis, :] + 1).T.flatten() //2" @@ -119,13 +120,13 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ - "def correlate(sequence, nbits=5, code=29, decimation=1):\n", + "def correlate(sequence, nbits=5, code=29, decimation=1, mask_filter=lambda x: x):\n", " # 0, 1 -> -1, 1\n", - " mask = np.repeat(gold(nbits)[code]*2 -1, decimation)\n", + " mask = mask_filter(np.repeat(gold(nbits)[code]*2 -1, decimation))\n", " # center\n", " sequence -= np.mean(sequence)\n", " return np.correlate(sequence, mask, mode='full')" @@ -133,7 +134,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 19, "metadata": {}, "outputs": [ { @@ -146,7 +147,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "eb7e7e5d7dfe4e00b18c4e5038c11182", + "model_id": "750255b26cc74aa0974d288fda4c7142", "version_major": 2, "version_minor": 0 }, @@ -167,10 +168,10 @@ { "data": { "text/plain": [ - "[<matplotlib.lines.Line2D at 0x7f0494537dc0>]" + "[<matplotlib.lines.Line2D at 0x7fe4aa733850>]" ] }, - "execution_count": 8, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -183,7 +184,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 20, "metadata": {}, "outputs": [ { @@ -196,7 +197,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "c693349edfe843d6adab192d6a95c4dd", + "model_id": "996e7034b52d47409ad4548c354b72bd", "version_major": 2, "version_minor": 0 }, @@ -218,10 +219,10 @@ { "data": { "text/plain": [ - "(2.0, 1.0311014124075548)" + "(2.0, 1.0490216904842018)" ] }, - "execution_count": 9, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -265,7 +266,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 21, "metadata": {}, "outputs": [ { @@ -285,7 +286,7 @@ }, { "cell_type": "code", - "execution_count": 54, + "execution_count": 22, "metadata": {}, "outputs": [ { @@ -294,21 +295,14 @@ "text": [ "(63,) (63,)\n", "(63,) (63,)\n", + "(63,) (63,)\n", "(63,) (63,)\n" ] }, { - "name": "stderr", - "output_type": "stream", - "text": [ - "<ipython-input-54-34e6ee3f3fc5>:22: 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, ((ax1, ax3, ax5), (ax2, ax4, ax6)) = plt.subplots(2, 3, figsize=(16, 9))\n" - ] - }, - { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "f58125333c294cb1b426b735829c30c5", + "model_id": "abdce9c2d307402f8578eb83d9ce9b79", "version_major": 2, "version_minor": 0 }, @@ -325,7 +319,7 @@ "(0.002, 0.012591236)" ] }, - "execution_count": 54, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -338,19 +332,21 @@ "foo = np.repeat(modulate(np.array([0, 1, 0, 0, 1, 1, 1, 0]), nbits) * 2.0 - 1, decimation) * signal_amplitude\n", "noise = np.resize(mains_noise, len(foo))\n", "\n", - "sosh = sig.butter(12, 0.05, btype='highpass', output='sos', fs=decimation)\n", - "sosl = sig.butter(6, 1.0, btype='lowpass', output='sos', fs=decimation)\n", + "sosh = sig.butter(3, 0.01, btype='highpass', output='sos', fs=decimation)\n", + "sosl = sig.butter(3, 0.8, btype='lowpass', output='sos', fs=decimation)\n", "#filtered = sig.sosfilt(sosh, sig.sosfilt(sosl, foo + noise))\n", "filtered = sig.sosfilt(sosh, foo + noise)\n", "\n", "cor1 = correlate(foo + noise, nbits=nbits, decimation=decimation)\n", "cor2 = correlate(filtered, nbits=nbits, decimation=decimation)\n", "\n", + "cor2_pe = correlate(filtered, nbits=nbits, decimation=decimation, mask_filter=lambda mask: sig.sosfilt(sosh, sig.sosfiltfilt(sosl, mask)))\n", + "\n", "sosn = sig.butter(12, 4, btype='highpass', output='sos', fs=decimation)\n", "#cor1_flt = sig.sosfilt(sosn, cor1)\n", "#cor2_flt = sig.sosfilt(sosn, cor2)\n", - "cor1_flt = cor1[1:] - cor1[:-1]\n", - "cor2_flt = cor2[1:] - cor2[:-1]\n", + "#cor1_flt = cor1[1:] - cor1[:-1]\n", + "#cor2_flt = cor2[1:] - cor2[:-1]\n", "\n", "fig, ((ax1, ax3, ax5), (ax2, ax4, ax6)) = plt.subplots(2, 3, figsize=(16, 9))\n", "fig.tight_layout()\n", @@ -371,17 +367,265 @@ "ax4.set_title('corr filtered')\n", "ax4.grid()\n", "\n", - "ax5.plot(cor1_flt)\n", - "ax5.set_title('corr raw (highpass)')\n", - "ax5.grid()\n", + "#ax5.plot(cor1_flt)\n", + "#ax5.set_title('corr raw (highpass)')\n", + "#ax5.grid()\n", + "\n", + "#ax6.plot(cor2_flt)\n", + "#ax6.set_title('corr filtered (highpass)')\n", + "#ax6.grid()\n", "\n", - "ax6.plot(cor2_flt)\n", - "ax6.set_title('corr filtered (highpass)')\n", + "ax6.plot(cor2_pe)\n", + "ax6.set_title('corr filtered w/ mask preemphasis')\n", "ax6.grid()\n", "\n", "rms = lambda x: np.sqrt(np.mean(np.square(x)))\n", "rms(foo), rms(noise)" ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "8feea8e305004d33a39f2541a59a0ffa", + "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": "stdout", + "output_type": "stream", + "text": [ + "(63,) (63,)\n" + ] + }, + { + "data": { + "text/plain": [ + "[<matplotlib.lines.Line2D at 0x7fe4a9934c40>]" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "seq = np.repeat(gold(6)[29]*2 -1, decimation)\n", + "sosh = sig.butter(3, 0.01, btype='highpass', output='sos', fs=decimation)\n", + "sosl = sig.butter(3, 0.8, btype='lowpass', output='sos', fs=decimation)\n", + "seq_filtered = sig.sosfilt(sosh, sig.sosfiltfilt(sosl, seq))\n", + "#seq_filtered = sig.sosfilt(sosh, seq)\n", + "\n", + "ax.plot(seq)\n", + "ax.plot(seq_filtered)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "01df90b27d57470d9216ae9c549413c8", + "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": "stdout", + "output_type": "stream", + "text": [ + "(63,) (63,)\n" + ] + }, + { + "data": { + "text/plain": [ + "[<matplotlib.lines.Line2D at 0x7fe4a9596070>]" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig, axs = plt.subplots(3, 1, figsize=(9, 7), sharex=True)\n", + "fig.tight_layout()\n", + "axs = axs.flatten()\n", + "for ax in axs:\n", + " ax.grid()\n", + "\n", + "seq = np.repeat(gold(6)[29]*2 -1, decimation)\n", + "sosh = sig.butter(3, 0.1, btype='highpass', output='sos', fs=decimation)\n", + "sosl = sig.butter(3, 0.8, btype='lowpass', output='sos', fs=decimation)\n", + "cor2_pe_flt = sig.sosfilt(sosh, cor2_pe)\n", + "cor2_pe_flt2 = sig.sosfilt(sosh, sig.sosfiltfilt(sosl, cor2_pe))\n", + "\n", + "axs[0].plot(cor2_pe)\n", + "axs[1].plot(cor2_pe_flt)\n", + "axs[2].plot(cor2_pe_flt2)\n", + "\n", + "#for idx in np.where(np.abs(cor2_pe_flt2) > 0.5)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<ipython-input-49-9776d553457e>:5: 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, (ax1, ax2, ax3) = plt.subplots(3, figsize=(12, 12))\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "183f086844054252bc8ec77f7b99abfc", + "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": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "30db7a9c75834d64a31bf8e0fb5f5911", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=10.0, description='w', max=30.0, min=-10.0), Output()), _dom_classes=(…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "[<matplotlib.lines.Line2D at 0x7fe482e19190>]" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "threshold_factor = 5.0\n", + "\n", + "import ipywidgets\n", + "\n", + "cor_an = cor1\n", + "\n", + "fig, (ax1, ax2, ax3) = plt.subplots(3, figsize=(12, 12))\n", + "fig.tight_layout()\n", + "\n", + "ax1.matshow(sig.cwt(cor_an, sig.ricker, np.arange(1, 31)), aspect='auto')\n", + "\n", + "for i in np.linspace(1, 10, 19):\n", + " offx = 5*i\n", + " ax2.plot(sig.cwt(cor_an, sig.ricker, [i]).flatten() + offx, color='red')\n", + "\n", + " ax2.text(-50, offx, f'{i:.1f}',\n", + " horizontalalignment='right',\n", + " verticalalignment='center',\n", + " color='black')\n", + "ax2.grid()\n", + "\n", + "ax3.grid()\n", + "\n", + "cwt_res = sig.cwt(cor_an, sig.ricker, [7.3]).flatten()\n", + "line, = ax3.plot(cwt_res)\n", + "def update(w=10.0):\n", + " line.set_ydata(sig.cwt(cor_an, sig.ricker, [w]).flatten())\n", + " fig.canvas.draw_idle()\n", + "ipywidgets.interact(update)\n", + "\n", + "import itertools\n", + "th = np.convolve(np.abs(cwt_res), np.ones((500,))/500, mode='same')\n", + "peaks = [ list(group) for val, group in itertools.groupby(enumerate(zip(th, cwt_res)), lambda elem: abs(elem[1][1]) > elem[1][0]*threshold_factor) if val ]\n", + "for group in peaks:\n", + " pos = np.mean([idx for idx, _val in group])\n", + " pol = np.mean([val for _idx, (_th, val) in group])\n", + " ax3.axvline(pos, color='red', alpha=0.5)\n", + " ax3.text(pos-20, 2.0, f'{0 if pol < 0 else 1}', horizontalalignment='right', verticalalignment='center', color='black')\n", + " \n", + "ax3.plot(th)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "43b48925433b4464928805812cfebc24", + "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 0x7fe4a120ba60>]" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig, axs = plt.subplots(2, 1, figsize=(9, 7))\n", + "fig.tight_layout()\n", + "axs = axs.flatten()\n", + "for ax in axs:\n", + " ax.grid()\n", + " \n", + "axs[0].plot(cor2_pe_flt2[1::10] - cor2_pe_flt2[:-1:10])\n", + "a, b = cor2_pe_flt2[1::10] - cor2_pe_flt2[:-1:10], np.array([0.0, -0.5, 1.0, -0.5, 0.0])\n", + "axs[1].plot(np.convolve(a, b, mode='full'))" + ] } ], "metadata": { |