{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "import numpy as np\n", "from scipy import signal as sig\n", "import struct\n", "import random\n", "import ipywidgets\n", "import itertools\n", "\n", "import colorednoise\n", "\n", "np.set_printoptions(linewidth=240)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib widget" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#colorednoise.powerlaw_psd_gaussian(1, int(1e4))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# From https://github.com/mubeta06/python/blob/master/signal_processing/sp/gold.py\n", "preferred_pairs = {5:[[2],[1,2,3]], 6:[[5],[1,4,5]], 7:[[4],[4,5,6]],\n", " 8:[[1,2,3,6,7],[1,2,7]], 9:[[5],[3,5,6]], \n", " 10:[[2,5,9],[3,4,6,8,9]], 11:[[9],[3,6,9]]}\n", "\n", "def gen_gold(seq1, seq2):\n", " print(seq1.shape, seq2.shape)\n", " gold = [seq1, seq2]\n", " for shift in range(len(seq1)):\n", " gold.append(seq1 ^ np.roll(seq2, -shift))\n", " return gold\n", "\n", "def gold(n):\n", " n = int(n)\n", " if not n in preferred_pairs:\n", " raise KeyError('preferred pairs for %s bits unknown' % str(n))\n", " t0, t1 = preferred_pairs[n]\n", " (seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)\n", " return gen_gold(seq0, seq1)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "01394154cc52483e9d4483f1178d94c3", "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": [ "(31,) (31,)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.matshow(gold(5))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def modulate(data, nbits=5):\n", " # 0, 1 -> -1, 1\n", " mask = np.array(gold(nbits))*2 - 1\n", " \n", " sel = mask[data>>1]\n", " data_lsb_centered = ((data&1)*2 - 1)\n", "\n", " return (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(31,) (31,)\n" ] }, { "data": { "text/plain": [ "array([-1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1,\n", " 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1,\n", " 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, 1, 1, 1,\n", " 1, -1, 1, 1, 1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1,\n", " -1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1,\n", " 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1,\n", " 1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1,\n", " -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1,\n", " -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.array(list(range(16)))\n", "\n", "mask = np.array(gold(5))*2 - 1\n", " \n", "sel = mask[data>>1]\n", "data_lsb_centered = ((data&1)*2 - 1)\n", "mask.shape, data.shape, sel.shape\n", "\n", "#fig, ax = plt.subplots()\n", "#ax.plot(\n", "np.multiply(sel, np.tile(data_lsb_centered, (2**5-1, 1)).T).flatten()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def correlate(sequence, nbits=5, decimation=1, mask_filter=lambda x: x):\n", " mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))\n", "\n", " sequence -= np.mean(sequence)\n", " \n", " return np.array([np.correlate(sequence, row, mode='full') for row in mask])" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(31,) (31,)\n", "(31,) (31,)\n", "shapes (1240,) (1240,)\n", "(31,) (31,)\n", "mask (33, 310)\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "182fd5ac86e74ad299a67e5f1d0b2b2b", "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": [ "(31,) (31,)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nbits = 5\n", "decimation = 10\n", "\n", "foo = np.repeat(modulate(np.array(list(range(4))), nbits).astype(float), decimation)\n", "bar = np.repeat(modulate(np.array(list(range(4))), nbits) * 2.0 - 1, decimation) * 1e-3\n", "print('shapes', foo.shape, bar.shape)\n", "\n", "mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))\n", "print('mask', mask.shape)\n", "\n", "fig, (ax1, ax2) = plt.subplots(2, figsize=(16, 5))\n", "fig.tight_layout()\n", "corr_m = np.array([np.correlate(foo, row, mode='full') for row in mask])\n", "#corr_m = np.array([row for row in mask])\n", "ax1.matshow(corr_m, aspect='auto')\n", "#ax.matshow(foo.reshape(32, 31)[::2,:])\n", "ax2.matshow(correlate(bar, decimation=decimation), aspect='auto')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "4cb2661eebb84478b06d285166ec13bc", "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": [ "(31,) (31,)\n", "(31,) (31,)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "decimation = 10\n", "\n", "fig, (ax1, ax2) = plt.subplots(2, figsize=(12, 5))\n", "fig.tight_layout()\n", "\n", "#mask = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((2**nbits + 1, (2**nbits-1) * decimation))\n", "#mask_stretched = np.tile(np.array(gold(nbits))[:,:,np.newaxis]*2 - 1, (1, 1, 1)).reshape((2**nbits + 1, (2**nbits-1) * 1))\n", "\n", "#ax1.matshow(mask)\n", "#ax2.matshow(mask_stretched, aspect='auto')\n", "\n", "foo = np.repeat(modulate(np.array(list(range(4)))).astype(float), 1).reshape((4, 31))\n", "foo_stretched = np.repeat(modulate(np.array(list(range(4)))).astype(float), 10).reshape(4, 310)\n", "\n", "ax1.matshow(foo)\n", "ax2.matshow(foo_stretched, aspect='auto')" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(31,) (31,)\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a2e2f747193b478bbfa792a0995ad4ed", "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": [ "(31,) (31,)\n", "(31,) (31,)\n" ] }, { "data": { "text/plain": [ "(2.0, 1.0234353995297893)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "decimation = 10\n", "signal_amplitude = 2.0\n", "nbits = 5\n", "\n", "foo = np.repeat(modulate(np.array([0, 1, 0, 0, 1, 1, 1, 0]), nbits) * 2.0 - 1, decimation) * signal_amplitude\n", "noise = colorednoise.powerlaw_psd_gaussian(1, len(foo))\n", "\n", "sosh = sig.butter(4, 0.01, 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, foo + noise))\n", "#filtered = sig.sosfilt(sosh, foo + noise)\n", "\n", "fig, ((ax1, ax3), (ax2, ax4)) = plt.subplots(2, 2, figsize=(16, 9))\n", "fig.tight_layout()\n", "\n", "ax1.plot(foo + noise)\n", "ax1.plot(foo)\n", "ax1.set_title('raw')\n", "\n", "ax2.plot(filtered)\n", "ax2.plot(foo)\n", "ax2.set_title('filtered')\n", "\n", "ax3.plot(correlate(foo + noise, nbits=nbits, decimation=decimation))\n", "ax3.set_title('corr raw')\n", " \n", "ax3.grid()\n", "\n", "ax4.plot(correlate(filtered, nbits=nbits, decimation=decimation))\n", "ax4.set_title('corr filtered')\n", "ax4.grid()\n", "\n", "rms = lambda x: np.sqrt(np.mean(np.square(x)))\n", "rms(foo), rms(noise)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean: 49.98625\n" ] } ], "source": [ "with open('/mnt/c/Users/jaseg/shared/raw_freq.bin', 'rb') as f:\n", " mains_noise = np.copy(np.frombuffer(f.read(), dtype='float32'))\n", " print('mean:', np.mean(mains_noise))\n", " mains_noise -= np.mean(mains_noise)" ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(31,) (31,)\n", "(31,) (31,)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":33: 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), (ax2, ax4)) = plt.subplots(2, 2, figsize=(16, 9))\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "246c19a3c1424e7eb0b675ff060ea5b3", "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": [ "(0.002, 0.013899708)" ] }, "execution_count": 96, "metadata": {}, "output_type": "execute_result" } ], "source": [ "decimation = 10\n", "signal_amplitude = 2.0e-3\n", "nbits = 5\n", "\n", "#test_data = np.random.randint(0, 2, 100)\n", "#test_data = np.array([0, 1, 0, 0, 1, 1, 1, 0])\n", "#test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**nbits), 64)\n", "#test_data = np.random.RandomState(seed=0).randint(0, 8, 64)\n", "#test_data = np.array(list(range(8)) * 8)\n", "#test_data = np.array([0, 1] * 32)\n", "test_data = np.array(list(range(64)))\n", "\n", "foo = np.repeat(modulate(test_data, nbits) * 2.0 - 1, decimation) * signal_amplitude\n", "noise = np.resize(mains_noise, len(foo))\n", "#noise = 0\n", "\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", "\n", "fig, ((ax1, ax3), (ax2, ax4)) = plt.subplots(2, 2, figsize=(16, 9))\n", "fig.tight_layout()\n", "\n", "ax1.plot(foo + noise)\n", "ax1.plot(foo)\n", "ax1.set_title('raw')\n", "\n", "ax2.plot(filtered)\n", "ax2.plot(foo)\n", "ax2.set_title('filtered')\n", "\n", "ax3.plot(cor1.T)\n", "ax3.set_title('corr raw')\n", "ax3.grid()\n", "\n", "#ax4.plot(cor2[:4].T)\n", "#ax4.set_title('corr filtered')\n", "#ax4.grid()\n", "ax4.matshow(cor1, aspect='auto')\n", "\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_pe[:4].T)\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": 14, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c08b2a1dbdef429eb22b598bd3dc0146", "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": [ "[]" ] }, "execution_count": 14, "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": 15, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "20117316e02548a99386c39e45e71ef1", "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" ] }, { "ename": "NameError", "evalue": "name 'cor2_pe' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0msosh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbutter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'highpass'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'sos'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdecimation\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[0msosl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbutter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.8\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'lowpass'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'sos'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdecimation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mcor2_pe_flt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msosfilt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msosh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcor2_pe\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 11\u001b[0m \u001b[0mcor2_pe_flt2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msosfilt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msosh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msosfiltfilt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msosl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcor2_pe\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 12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'cor2_pe' is not defined" ] } ], "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": 57, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":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()\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "25fe274dea83415491fd7f86d38188d7", "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": [ "[]" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, ax = plt.subplots()\n", "nonlinear_distance = lambda x: 100**(2*np.abs(0.5-x%1)) / (np.abs(x)+3)**2\n", "x = np.linspace(-1.5, 5.5, 10000)\n", "ax.plot(x, nonlinear_distance(x))" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":11: 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) = plt.subplots(2, figsize=(12, 5))\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "b083c661b5b441d6b7fc45201faa0576", "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": [ "cor_an (33, 20149)\n", "cwt_res (33, 20149)\n", "th (33, 20149)\n", "[((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,)), ((33,), (33,))]\n", "peaks: 180\n", "avg_peak 1.058897833824206\n", "skipped 3 symbols at 5889.0\n", "skipped 2 symbols at 8369.0\n", "skipped 2 symbols at 14568.5\n", "skipped 2 symbols at 16739.0\n", "decoding [ref|dec]:\n", " -1| -1 ✔ 1| 1 ✔ 2| 2 ✔ 3| 3 ✔ 4| 4 ✔ 5| 5 ✔ 6| 6 ✔ 7| 7 ✔ \n", " 8| 8 ✔ 9| 9 ✔ 10| 10 ✔ 11| 11 ✔ 12| 12 ✔ 13| 13 ✔ 14| 14 ✔ 15| 15 ✔ \n", " 16| -1 ✘ 17| -1 ✘ 18| 18 ✔ 19| 19 ✔ 20| 20 ✔ 21| 21 ✔ 22| 22 ✔ 23| 23 ✔ \n", " 24| 24 ✔ 25| -1 ✘ 26| 26 ✔ 27| 27 ✔ 28| 28 ✔ 29| 29 ✔ 30| 30 ✔ 31| 31 ✔ \n", " 32| 32 ✔ 33| 33 ✔ 34| 34 ✔ 35| 35 ✔ 36| 36 ✔ 37| 37 ✔ 38| 38 ✔ 39| 39 ✔ \n", " 40| 40 ✔ 41| 41 ✔ 42| 42 ✔ 43| 43 ✔ 44| 44 ✔ 45| -1 ✘ 46| 46 ✔ 47| 47 ✔ \n", " 48| 48 ✔ 49| 49 ✔ 50| 50 ✔ 51| 51 ✔ 52| -1 ✘ 53| 53 ✔ 54| 54 ✔ 55| 55 ✔ \n", " 56| 56 ✔ 57| 57 ✔ 58| 58 ✔ 59| 59 ✔ 60| 60 ✔ 61| 61 ✔ 62| 62 ✔ 63| 56 ✘ \n", "Symbol error rate r=0.09375\n" ] } ], "source": [ "threshold_factor = 4.0\n", "power_avg_width = 1024\n", "max_lookahead = 6.5\n", "\n", "bit_period = (2**nbits) * decimation\n", "peak_group_threshold = 0.1 * bit_period\n", "\n", "cor_an = cor1\n", "\n", "#fig, (ax1, ax2, ax3) = plt.subplots(3, figsize=(12, 12))\n", "fig, (ax1, ax3) = plt.subplots(2, figsize=(12, 5))\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", "print('cor_an', cor_an.shape)\n", "\n", "cwt_res = np.array([ sig.cwt(row, sig.ricker, [0.73 * decimation]).flatten() for row in cor_an ])\n", "ax3.plot(cwt_res.T)\n", "#def update(w = 1.0 * decimation):\n", "# line.set_ydata(sig.cwt(cor_an, sig.ricker, [w]).flatten())\n", "# fig.canvas.draw_idle()\n", "#ipywidgets.interact(update)\n", "\n", "print('cwt_res', cwt_res.shape)\n", "th = np.array([ np.convolve(np.abs(row), np.ones((power_avg_width,))/power_avg_width, mode='same') for row in cwt_res ])\n", "ax1.plot(th.T)\n", "print('th', th.shape)\n", "\n", "def compare_th(elem):\n", " idx, (th, val) = elem\n", " #print('compare_th:', th.shape, val.shape)\n", " return np.any(np.abs(val) > th*threshold_factor)\n", "\n", "print([ (a.shape, b.shape) for a, b in zip(th.T, cwt_res.T) ][:5])\n", "\n", "peaks = [ list(group) for val, group in itertools.groupby(enumerate(zip(th.T, cwt_res.T)), compare_th) if val ]\n", "print('peaks:', len(peaks))\n", "peak_group = []\n", "for group in peaks:\n", " pos = np.mean([idx for idx, _val in group])\n", " pol = np.mean([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group])\n", " pol_idx = np.argmax(np.bincount([ np.argmax(np.abs(val)) for _idx, (_th, val) in group ]))\n", " #print(f'group', pos, pol, pol_idx)\n", " #for pol, (_idx, (_th, val)) in zip([max(val.min(), val.max(), key=abs) for _idx, (_th, val) in group], group):\n", " # print(' ', pol, val)\n", " ax3.axvline(pos, color='cyan', alpha=0.3)\n", " \n", " if not peak_group or pos - peak_group[-1][1] > peak_group_threshold:\n", " if peak_group:\n", " peak_pos = peak_group[-1][3]\n", " ax3.axvline(peak_pos, color='red', alpha=0.6)\n", " #ax3.text(peak_pos-20, 2.0, f'{0 if pol < 0 else 1}', horizontalalignment='right', verticalalignment='center', color='black')\n", " \n", " peak_group.append((pos, pos, pol, pos, pol_idx))\n", " #ax3.axvline(pos, color='cyan', alpha=0.5)\n", " \n", " else:\n", " group_start, last_pos, last_pol, peak_pos, last_pol_idx = peak_group[-1]\n", " \n", " if abs(pol) > abs(last_pol):\n", " #ax3.axvline(pos, color='magenta', alpha=0.5)\n", " peak_group[-1] = (group_start, pos, pol, pos, pol_idx)\n", " else:\n", " #ax3.axvline(pos, color='blue', alpha=0.5)\n", " peak_group[-1] = (group_start, pos, last_pol, peak_pos, last_pol_idx)\n", "\n", "avg_peak = np.mean(np.abs(np.array([last_pol for _1, _2, last_pol, _3, _4 in peak_group])))\n", "print('avg_peak', avg_peak)\n", "\n", "noprint = lambda *args, **kwargs: None\n", "def mle_decode(peak_groups, print=print):\n", " peak_groups = [ (pos, pol, idx) for _1, _2, pol, pos, idx in peak_groups ]\n", " candidates = [ (0, [(pos, pol, idx)]) for pos, pol, idx in peak_groups if pos < bit_period*2.5 ]\n", " \n", " while candidates:\n", " chain_candidates = []\n", " for chain_score, chain in candidates:\n", " pos, ampl, _idx = chain[-1]\n", " score_fun = lambda pos, npos, npol: abs(npol)/avg_peak + nonlinear_distance((npos-pos)/bit_period)\n", " next_candidates = sorted([ (score_fun(pos, npos, npol), npos, npol, nidx) for npos, npol, nidx in peak_groups if pos < npos < pos + bit_period*max_lookahead ], reverse=True)\n", " \n", " print(f' candidates for {pos}, {ampl}:')\n", " for score, npos, npol, nidx in next_candidates:\n", " print(f' {score:.4f} {npos:.2f} {npol:.2f} {nidx:.2f}')\n", " \n", " nch, cor_len = cor_an.shape\n", " if cor_len - pos < 1.5*bit_period or not next_candidates:\n", " score = sum(score_fun(opos, npos, npol) for (opos, _opol, _oidx), (npos, npol, _nidx) in zip(chain[:-1], chain[1:])) / len(chain)\n", " yield score, chain\n", " \n", " else:\n", " print('extending')\n", " for score, npos, npol, nidx in next_candidates[:3]:\n", " if score > 0.5:\n", " new_chain_score = chain_score * 0.9 + score * 0.1\n", " chain_candidates.append((new_chain_score, chain + [(npos, npol, nidx)]))\n", " print('chain candidates:')\n", " for score, chain in sorted(chain_candidates, reverse=True):\n", " print(' ', [(score, [(f'{pos:.2f}', f'{pol:.2f}') for pos, pol, _idx in chain])])\n", " candidates = [ (chain_score, chain) for chain_score, chain in sorted(chain_candidates, reverse=True)[:10] ]\n", "\n", "res = sorted(mle_decode(peak_group, print=noprint), reverse=True)\n", "#for i, (score, chain) in enumerate(res):\n", "# print(f'Chain {i}@{score:.4f}: {chain}')\n", "(_score, chain), *_ = res\n", "\n", "def viz(chain):\n", " last_pos = None\n", " for pos, pol, nidx in chain:\n", " if last_pos:\n", " delta = int(round((pos - last_pos) / bit_period))\n", " if delta > 1:\n", " print(f'skipped {delta} symbols at {pos}')\n", " for i in range(delta-1):\n", " yield None\n", " ax3.axvline(pos, color='blue', alpha=0.5)\n", " decoded = nidx*2 + (0 if pol < 0 else 1)\n", " yield decoded\n", " ax3.text(pos-20, 0.0, f'{decoded}', horizontalalignment='right', verticalalignment='center', color='black')\n", "\n", " last_pos = pos\n", "\n", "decoded = list(viz(chain))\n", "print('decoding [ref|dec]:')\n", "failures = 0\n", "for i, (ref, found) in enumerate(itertools.zip_longest(test_data, decoded)):\n", " print(f'{ref or -1:>3d}|{found or -1:>3d} {\"✔\" if ref==found else \"✘\"}', end=' ')\n", " if ref != found:\n", " failures += 1\n", " if i%8 == 7:\n", " print()\n", "print(f'Symbol error rate r={failures/len(test_data)}')\n", "#ax3.plot(th)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "987be038c1b34e6e9509f7f224bbb620", "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": [ "[]" ] }, "execution_count": 15, "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": { "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 }