diff options
-rw-r--r-- | lab-windows/dsss_experiments.ipynb | 408 |
1 files changed, 408 insertions, 0 deletions
diff --git a/lab-windows/dsss_experiments.ipynb b/lab-windows/dsss_experiments.ipynb new file mode 100644 index 0000000..bb56362 --- /dev/null +++ b/lab-windows/dsss_experiments.ipynb @@ -0,0 +1,408 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 10, + "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", + "\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": "7560730a2391425ab9dad7a1f22e5fb2", + "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": [ + "<matplotlib.image.AxesImage at 0x7f0496c50b80>" + ] + }, + "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, code=29):\n", + " # 0, 1 -> -1, 1\n", + " mask = gold(nbits)[code]*2 - 1\n", + " # same here\n", + " data_centered = (data*2 - 1)\n", + " return (mask[:, np.newaxis] @ data_centered[np.newaxis, :] + 1).T.flatten() //2" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def correlate(sequence, nbits=5, code=29, decimation=1):\n", + " # 0, 1 -> -1, 1\n", + " mask = np.repeat(gold(nbits)[code]*2 -1, decimation)\n", + " # center\n", + " sequence -= np.mean(sequence)\n", + " return np.correlate(sequence, mask, mode='full')" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(31,) (31,)\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "eb7e7e5d7dfe4e00b18c4e5038c11182", + "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": [ + "[<matplotlib.lines.Line2D at 0x7f0494537dc0>]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "foo = modulate(np.array([0, 1, 0, 0, 1, 1, 1, 0])).astype(float)\n", + "fig, ax = plt.subplots()\n", + "ax.plot(correlate(foo))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(31,) (31,)\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c693349edfe843d6adab192d6a95c4dd", + "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.0311014124075548)" + ] + }, + "execution_count": 9, + "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": 19, + "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": 54, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(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", + "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.012591236)" + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "decimation = 10\n", + "signal_amplitude = 2.0e-3\n", + "nbits = 6\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 = 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", + "#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", + "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, ax5), (ax2, ax4, ax6)) = plt.subplots(2, 3, 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)\n", + "ax3.set_title('corr raw')\n", + "ax3.grid()\n", + "\n", + "ax4.plot(cor2)\n", + "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", + "\n", + "ax6.plot(cor2_flt)\n", + "ax6.set_title('corr filtered (highpass)')\n", + "ax6.grid()\n", + "\n", + "rms = lambda x: np.sqrt(np.mean(np.square(x)))\n", + "rms(foo), rms(noise)" + ] + } + ], + "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 +} |