demod prototype

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": {