diff options
Diffstat (limited to 'notebooks/Modulation Experiment Calculations.ipynb')
| -rw-r--r-- | notebooks/Modulation Experiment Calculations.ipynb | 680 |
1 files changed, 680 insertions, 0 deletions
diff --git a/notebooks/Modulation Experiment Calculations.ipynb b/notebooks/Modulation Experiment Calculations.ipynb new file mode 100644 index 0000000..37e7ccb --- /dev/null +++ b/notebooks/Modulation Experiment Calculations.ipynb @@ -0,0 +1,680 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import struct\n", + "import random\n", + "import itertools\n", + "import datetime\n", + "import multiprocessing\n", + "from collections import defaultdict\n", + "import json\n", + "import traceback\n", + "import glob\n", + "\n", + "from matplotlib import pyplot as plt\n", + "import matplotlib\n", + "from matplotlib import ticker\n", + "import numpy as np\n", + "from scipy import signal as sig\n", + "from scipy import fftpack as fftpack\n", + "import ipywidgets\n", + "from IPython.display import set_matplotlib_formats\n", + "import scipy\n", + "\n", + "from tqdm.notebook import tqdm\n", + "import colorednoise\n", + "import accupy\n", + "\n", + "np.set_printoptions(linewidth=240)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "#%matplotlib inline\n", + "#set_matplotlib_formats('png', 'pdf')\n", + "#font = {'family' : 'normal',\n", + "# 'weight' : 'normal',\n", + "# 'size' : 6}\n", + "#matplotlib.rc('font', **font)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "sampling_rate = 10 # sp/s" + ] + }, + { + "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", + " 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", + " \n", + " min_key, max_key = min(preferred_pairs.keys()), max(preferred_pairs.keys())\n", + " \n", + " if n < min_key:\n", + " raise KeyError('preferred pairs for %s bits unknown' % str(n))\n", + " \n", + " if n > max_key:\n", + " t0, t1 = preferred_pairs[max_key]\n", + " (seq0, _st0), (seq1, _st1) = sig.max_len_seq(max_key, taps=t0), sig.max_len_seq(max_key, taps=t1)\n", + " raw_codes = gen_gold(seq0, seq1)\n", + " # Hackety hack: This is super-experimental test code and guaranteed sub-optimal!\n", + " \n", + " raw_codes = np.array(raw_codes)\n", + " out = raw_codes[:-1]\n", + " for i in range(n - max_key):\n", + " out = np.hstack([out[:out.shape[0]//2], out[out.shape[0]//2:]])\n", + " return out\n", + " \n", + " else:\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 np.array(gen_gold(seq0, seq1))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def modulate(data, nbits=5, pad=True):\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", + " signal = (np.multiply(sel, np.tile(data_lsb_centered, (sel.shape[1], 1)).T).flatten() + 1) // 2\n", + " if pad:\n", + " print('mod: padding by', mask.shape[1])\n", + " return np.hstack([ np.zeros(mask.shape[1]), signal, np.zeros(mask.shape[1]) ])\n", + " else:\n", + " return signal" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def correlate(sequence, nbits=5, decimation=1, mask_filter=lambda x: x, mask_max=None):\n", + " codes = gold(nbits)\n", + " mask = np.tile(np.array(codes)[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((codes.shape[0], codes.shape[1] * decimation))\n", + " mask = mask.astype(np.float64)\n", + "\n", + " # Our input signal has large DC bias. Remove DC bias to reduce numerical errors during correlation.\n", + " sequence -= np.mean(sequence)\n", + " \n", + " #return np.array([np.correlate(sequence, row, mode='valid') for row in mask[:mask_max]])\n", + " return np.array([sig.correlate(sequence, row, mode='valid') for row in mask[:mask_max]]) / len(sequence)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mean: 50.00341 len: 1946174\n" + ] + } + ], + "source": [ + "with open('data/fmeas_export_ocxo_2day.bin', 'rb') as f:\n", + " meas_data = np.copy(np.frombuffer(f.read(), dtype='float32'))\n", + " print('mean:', np.mean(meas_data), 'len:', len(meas_data))\n", + " meas_data -= np.mean(meas_data)\n", + " meas_data = np.hstack([meas_data, meas_data[::-1]] * 50)\n", + " def mains_noise_meas(n):\n", + " last_valid = len(meas_data) - n\n", + " start = random.randint(0, last_valid)\n", + " return meas_data[start:start+n]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def load_noise_meas_params(capture_file):\n", + " with open(capture_file, 'rb') as f:\n", + " meas_data = np.copy(np.frombuffer(f.read(), dtype='float32'))\n", + " print('mean:', np.mean(meas_data), 'len:', len(meas_data))\n", + " meas_data -= np.mean(meas_data)\n", + " return {'meas_data': meas_data}\n", + "\n", + "def mains_noise_measured(n, meas_data):\n", + " last_valid = len(meas_data) - n\n", + " start = np.random.randint(last_valid)\n", + " return meas_data[start:start+n]\n", + "\n", + "def load_noise_synth_params(specfile):\n", + " with open(specfile) as f:\n", + " d = json.load(f)\n", + " return {'spl_x': np.linspace(*d['x_spec']),\n", + " 'spl_N': d['x_spec'][2],\n", + " 'psd_spl': (d['t'], d['c'], d['k']) }\n", + " \n", + "def mains_noise_synthetic(n, psd_spl, spl_N, spl_x):\n", + " noise = np.random.normal(size=spl_N) * 2\n", + " spec = scipy.fftpack.fft(noise) **2\n", + " spec *= np.exp(scipy.interpolate.splev(spl_x, psd_spl))\n", + "\n", + " spec **= 1/2\n", + " \n", + " renoise = scipy.fftpack.ifft(spec)\n", + " return np.hstack([renoise[10000:], renoise[10000:][::-1]])[:n]\n", + "\n", + "noise_params = load_noise_synth_params('data/grid_freq_psd_spl_108pt.json')\n", + "mains_noise = lambda n: mains_noise_synthetic(n, **noise_params)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "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, pad=True) * 2.0 - 1, decimation) * signal_amplitude\n", + " noise = mains_noise_meas(len(signal))\n", + " \n", + " return test_data, signal.astype(np.longdouble) + noise.astype(np.longdouble)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.50806058, 0.49242794])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gold(12)[:2].mean(axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "amplitude @P_mod=10000.0kW: Δf=400.0µHz\n", + "mod: padding by 1023\n", + "[0 1 2 3] (49104,)\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "70c8a0e9c03441fab580b5b6b219e065", + "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": [ + "signal strength (off/on peak): +4.900e-05 | -6.406e-04\n", + "signal strength (off/on peak): +7.699e-05 | +5.627e-04\n", + "signal strength (off/on peak): +8.441e-05 | +4.206e-05\n", + "signal strength (off/on peak): +8.178e-05 | -8.413e-05\n" + ] + } + ], + "source": [ + "nbits = 10\n", + "decimation = 8\n", + "# working params: decimation=15 nbits=21 power=200kW\n", + "# working params: decimation=8 nbits=10 power=10MW\n", + "grid_frequency_characteristic = 25e9 # watt per hertz\n", + "modulation_power = 10000e3 # watt\n", + "data = np.array([0, 1, 2, 3])\n", + "amplitude = modulation_power / grid_frequency_characteristic\n", + "print(f'amplitude @P_mod={modulation_power/1e3}kW: Δf={amplitude*1e6}µHz')\n", + "\n", + "# Generate test data for this simulation run\n", + "test_data, signal = generate_test_signal(duration=2, nbits=nbits, signal_amplitude=amplitude, decimation=decimation, seed=0, data=data)\n", + "print(test_data, signal.shape)\n", + "\n", + "def correlate_test(sequence, nbits, decimation, mask_max):\n", + " codes = gold(nbits)\n", + " mask = np.tile(np.array(codes)[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((codes.shape[0], codes.shape[1] * decimation))\n", + " mask = mask.astype(np.float64)\n", + "\n", + " # Our input signal has large DC bias. Remove DC bias to reduce numerical errors during correlation.\n", + " sequence -= np.mean(sequence)\n", + " \n", + " #return np.array([np.correlate(sequence, row, mode='valid') for row in mask[:mask_max]])\n", + " return np.array([accupy.fsum(np.multiply(sequence, row)) for row in mask[:mask_max]]) / len(sequence)\n", + "\n", + "signal = signal.astype(np.float64)\n", + "fig, axs = plt.subplots(len(data), figsize=(12, 5))\n", + "for i, ax in enumerate(axs.flatten()):\n", + " gen_per = len(data) + 2 # one period for each symbol, plus one period of padding at the beginning and end\n", + " symbol_length = signal.shape[0]//gen_per\n", + " l, r = int(max(0, (i + 1) * symbol_length - 1000)), int((i + 2) * symbol_length + 1000)\n", + " signal_slice = signal[l:r]\n", + " cor_an = correlate(signal_slice, nbits=nbits, decimation=decimation, mask_max=2)\n", + " \n", + " l, r = int((i + 1) * symbol_length), int((i + 2) * symbol_length)\n", + " signal_slice = signal[l:r]\n", + " #print('on peak [numpy]:', correlate(signal_slice, nbits=nbits, decimation=decimation, mask_max=2).flatten() / amplitude)\n", + " #print('on peak [accupy]:', correlate_test(signal_slice, nbits=nbits, decimation=decimation, mask_max=2).flatten() / amplitude)\n", + " \n", + " l, r = int((i + 1 - 0.2) * symbol_length), int((i + 2 - 0.2) * symbol_length)\n", + " signal_slice = signal[l:r]\n", + " #print('off peak [numpy]:', correlate(signal_slice, nbits=nbits, decimation=decimation, mask_max=2).flatten() / amplitude)\n", + " #print('off peak [accupy]:', correlate_test(signal_slice, nbits=nbits, decimation=decimation, mask_max=2).flatten() / amplitude)\n", + " #ax.plot(cor_an.T)\n", + "\n", + " cwt_res = np.array([ sig.cwt(row, sig.ricker, [0.73 * decimation]).flatten() for row in cor_an ])\n", + " print(f'signal strength (off/on peak): {np.sqrt(np.mean(cwt_res[0,cwt_res.shape[1]//6:cwt_res.shape[1]//3]**2)):+.3e} | {cwt_res[0,cwt_res.shape[1]//2]:+.3e}')\n", + " ax.plot(cwt_res.T / amplitude)" + ] + }, + { + "cell_type": "code", + "execution_count": 127, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "amplitude @P_mod=10.0kW: Δf=0.39999999999999997µHz\n", + "mod: padding by 31\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<ipython-input-95-7000bdc82d1b>:7: ComplexWarning: Casting complex values to real discards the imaginary part\n", + " return test_data, signal.astype(np.longdouble) + noise.astype(np.longdouble)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "mod: padding by 31\n", + "[0.00031921745429844274216, 0.0003466253347083378983, 0.00044946868526334726654, 0.00043601785453840703518, 0.00027364766834844621423, 0.00034429714350871021226, 0.00029257480354291570764, 0.00043761373963434835483, 0.00034653346037204395042, 0.00041369145870726819648, 0.00031663491690404936546, 0.0003630671579533060333, 0.0002934121262870590337, 0.00036264881740159871053, 0.00033973548163424448722, 0.00045000162460528714934, 0.00036345615545028226934, 0.00034203880459996953548, 0.0004340491996267533479, 0.00027402160564997715637]\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "mod: padding by 63\n", + "[0.00043183462240684953605, 0.00049995218014044821913, 0.00034320159159449236953, 0.00047574377773946273916, 0.00042495372363872566055, 0.0005340590070278771213, 0.00043128485159190824515, 0.00034867807874549918169, 0.00039072607874113340072, 0.00034884834981040312958, 0.00042730792396368059213, 0.0005043675584057222778, 0.0003499736642140912821, 0.00021263777836169804578, 0.00023831087214244111956, 0.00030008253740133703718, 0.00021268569995328232606, 0.00054148296194613603103, 0.00040150545723533495933, 0.00028661636535881167316]\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "mod: padding by 127\n", + "[0.00026772134709063096758, 0.0004023169589146378664, 0.00052784935496986912216, 0.00033900150125753280712, 0.00022029107067617116654, 0.00041990310449284588383, 0.0004415131846159276476, 0.00040143078518334457215, 0.00031109331861264841798, 0.0002811855751776715795, 0.00044570966535767569383, 0.0003301271164337369956, 0.00033158808534416807198, 0.00026578644644711210834, 0.00034742527826464402753, 0.0005426992707322571906, 0.00044335106857497668567, 0.00032225907153724712904, 0.0002813426161001733043, 0.00041113085027627075947]\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "mod: padding by 255\n", + "[0.00019884404254742679324, 0.000242880777230771204, 0.00019589281771018505908, 0.00022349702779183853734, 0.00020954727853424042772, 0.00020473570231859194066, 0.00020169541990839266527, 0.00025024395700750437305, 0.0002751261591544219094, 0.00019445654879589896521, 0.00020805938128876005823, 0.00020078107591647889947, 0.00019010551607677558, 0.00023584019352316058905, 0.0002678021299787742911, 0.00017845149465711736266, 0.00016695541802082155936, 0.00016698843247300128494, 0.00014251122213695801984, 0.00028090032636289179365]\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "mod: padding by 511\n", + "[0.00019975970170629950027, 0.00013186823465058419483, 0.00023953255381256545738, 0.00017255565602586637237, 0.00018396445334660804652, 0.00014540194644455677518, 0.00019876358084066945335, 0.0001738424376447688047, 0.0001765109371402779291, 0.00018196424414407873752, 0.00023044940092346266965, 0.00019731842839493123688, 0.00023705358802857269644, 0.00015328320915330276564, 0.00017778635677371888602, 0.00021212946232256081042, 0.00026136191723243009158, 0.00016766188832612281848, 0.000192498527508538207, 0.00024108248581102643531]\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "mod: padding by 1023\n", + "[0.00013530662801083781384, 0.0001497806916396048579, 0.00013664959261312633858, 0.00016589078984916737442, 0.00017011538645283205862, 0.00015336800379842749231, 0.0001694272184760911441, 0.00013875282881629243428, 0.000112969749633285926445, 0.00014931796422008114573, 0.0001372035855968592536, 0.00012312327359857258514, 0.00013905194521669023908, 0.00013901001422022015822, 0.00016339317685934381651, 0.00015711672887077335837, 0.00012437507736796685456, 0.00012540715390246491552, 0.00010657701913677854194, 0.00014676391008527084711]\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "mod: padding by 2047\n", + "[9.397038360574742737e-05, 9.414892638810825941e-05, 8.643228267549556685e-05, 8.4662709794848760505e-05, 7.018825855548282147e-05, 9.3962524208973514975e-05, 8.850854906632382727e-05, 9.490061324129898401e-05, 9.9421482075907846084e-05, 9.87805754535881666e-05, 8.77129165764244642e-05, 8.627947590758716296e-05, 9.310894721386338135e-05, 0.00010087838856131649663, 7.6747263039331168876e-05, 8.2710693121826614583e-05, 0.000106444803775372911015, 9.682856374335765839e-05, 8.172981177790056443e-05, 9.233645079286106117e-05]\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "mod: padding by 4094\n", + "[7.291179023159128505e-05, 6.921653169518604612e-05, 7.0212257862804909805e-05, 7.739735101427627361e-05, 7.05170185815625493e-05, 7.667140441838123699e-05, 6.684245257693558025e-05, 6.7274435477192927474e-05, 7.3140978179672748096e-05, 6.462540379490416508e-05, 7.8575870371873889974e-05, 7.0962264879147295656e-05, 6.887905626735870121e-05, 6.6741332451462875345e-05, 6.774053852476256235e-05, 8.4915545476866071356e-05, 6.755406828145308233e-05, 7.4150959965664274745e-05, 7.289475345040078811e-05, 7.417472792999099321e-05]\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "mod: padding by 8188\n", + "[4.809424705445009225e-05, 5.3333793693925671324e-05, 5.6174142194133624293e-05, 5.070037019380478062e-05, 4.2797138004346089653e-05, 5.6184084783894627412e-05, 4.8948695230438478363e-05, 5.231607985651584901e-05, 4.8563366184299055167e-05, 5.0229443138648695e-05, 5.140676133340082623e-05, 5.1571059935388156693e-05, 4.799175351462962519e-05, 5.0143994275182861474e-05, 5.176115169674431516e-05, 5.4781771632211411004e-05, 5.1420164246213953733e-05, 5.584101664179133697e-05, 5.2387951249217037658e-05, 4.8454927642515831705e-05]\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "mod: padding by 16376\n", + "[3.441001692841763238e-05, 3.9226152369629212293e-05, 3.5174358297006928404e-05, 3.858339666154567386e-05, 3.729609282512620058e-05, 3.8017607168678973304e-05, 3.6432437134237636557e-05, 3.5396590548610446987e-05, 4.2872778372966043538e-05, 3.881954270914244008e-05, 3.6638457035791679994e-05, 4.3987686908835178038e-05, 3.5810535332467513793e-05, 3.669588547824606978e-05, 3.9034162252655188724e-05, 3.686462573963708382e-05, 3.64018808798140893e-05, 4.140879713977473717e-05, 3.601891068378439981e-05, 4.2239933297790315508e-05]\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "mod: padding by 32752\n", + "[3.0164568033183993325e-05, 3.340002097279428693e-05, 3.1692825762129764792e-05, 3.2341558069057274648e-05, 3.5104645491986034208e-05, 3.3804378516131733095e-05, 3.0847430711483434965e-05, 3.163598515247249634e-05, 3.1877253973163090688e-05, 3.145134169534852806e-05, 2.9469157225219549716e-05, 3.5022547350051354115e-05, 3.185664970143889292e-05, 3.2250533412830754613e-05, 3.2968988013916950676e-05, 3.1167159682343826453e-05, 3.0519324313846599056e-05, 3.6983652848933391153e-05, 3.4137279541673856477e-05, 3.4656056379915872357e-05]\n" + ] + } + ], + "source": [ + "decimation = 10\n", + "grid_frequency_characteristic = 25e9 # watt per hertz\n", + "modulation_power = 10e3 # watt\n", + "amplitude = modulation_power / grid_frequency_characteristic\n", + "print(f'amplitude @P_mod={modulation_power/1e3}kW: Δf={amplitude*1e6}µHz')\n", + "\n", + "ps_off = []\n", + "xs = np.array(list(range(5, 16)))\n", + "for nbits in xs:\n", + " data = np.array([1])\n", + "\n", + " def correlate_test(sequence, nbits, decimation, mask_max):\n", + " codes = gold(nbits)\n", + " mask = np.tile(np.array(codes)[:,:,np.newaxis]*2 - 1, (1, 1, decimation)).reshape((codes.shape[0], codes.shape[1] * decimation))\n", + " mask = mask.astype(np.float64)\n", + "\n", + " # Our input signal has large DC bias. Remove DC bias to reduce numerical errors during correlation.\n", + " sequence -= np.mean(sequence)\n", + "\n", + " return np.array([sig.correlate(sequence, row, mode='valid') for row in mask[:mask_max]])\n", + " #return np.array([accupy.fsum(np.multiply(sequence, row)) for row in mask[:mask_max]]) / len(sequence)\n", + "\n", + " vals = []\n", + " for i in range(20):\n", + " test_data, signal = generate_test_signal(duration=2, nbits=nbits, signal_amplitude=amplitude, decimation=decimation, seed=random.randint(0, 0xffffffff), data=data)\n", + " cor_an = correlate(signal, nbits=nbits, decimation=decimation, mask_max=2)\n", + " vals.append(np.sqrt(np.mean(cor_an**2)))\n", + " ps_off.append((np.mean(vals), np.std(vals)))\n", + " print(vals)\n", + "\n", + "ps_off = np.array(ps_off)" + ] + }, + { + "cell_type": "code", + "execution_count": 129, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3dade90110af4097b123f3de508e843f", + "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": [ + "<ErrorbarContainer object of 3 artists>" + ] + }, + "execution_count": 129, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "ax.errorbar(x=xs, y=ps_off[:,0], yerr=ps_off[:,1])" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ma-thesis-env", + "language": "python", + "name": "ma-thesis-env" + }, + "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.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} |