Simu: 8 moved functions into lib

This commit is contained in:
Eric Teunis de Boone 2022-08-04 17:22:54 +02:00
parent d629dcc6eb
commit 657d1d1870

View file

@ -18,7 +18,13 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.signal as signal\n",
"import scipy.fft as ft"
"\n",
"import os\n",
"import sys\n",
"# Append parent directory to import path so lib can be found\n",
"sys.path.append(os.path.dirname(os.path.abspath(os.getcwd())))\n",
"from lib.util import *\n",
"from lib.plotting import *\n"
]
},
{
@ -120,89 +126,6 @@
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"### signal generation\n",
"def fft_bandpass(signal, band, samplerate):\n",
" \"\"\"\n",
" Simple bandpassing function employing a FFT.\n",
"\n",
" Parameters\n",
" ----------\n",
" signal : arraylike\n",
" band : tuple(low, high)\n",
" Frequencies for bandpassing\n",
" samplerate : float\n",
" \"\"\"\n",
" signal = np.asarray(signal)\n",
"\n",
" fft = ft.rfft(signal)\n",
" freqs = ft.rfftfreq(signal.size, 1/samplerate)\n",
" fft[(freqs < band[0]) | (freqs > band[1])] = 0\n",
" \n",
" return ft.irfft(fft, signal.size), (fft, freqs)\n",
"\n",
"def deltapeak(timelength=1e3, samplerate=1, offset=None, peaklength=1):\n",
" N_samples = int(timelength * samplerate)\n",
" if offset is None:\n",
" offset = (np.random.random(1)*N_samples).astype(int) % N_samples\n",
" elif isinstance(offset, (tuple, list)):\n",
" offset_min = offset[0]\n",
" offset_max = offset[-1]\n",
" \n",
" offset = (np.random.random(1)*(offset_max - offset_min)+offset_min).astype(int) % N_samples\n",
" \n",
" position = (offset + np.arange(0, peaklength)).astype(int) % N_samples\n",
" \n",
" signal = np.zeros(N_samples)\n",
" signal[position] = 1\n",
" \n",
" return signal, position\n",
"\n",
"def sin_delay(f, t, t_delay=0, phase=0):\n",
" return np.sin(2*np.pi*f * (t + t_delay) + phase)\n",
"\n",
"def annotate_width(ax, name, x1, x2, y, text_kw={}, arrow_kw={}):\n",
" default_arrow_kw = dict(\n",
" xy = (x1, y),\n",
" xytext = (x2,y),\n",
" arrowprops = dict(\n",
" arrowstyle=\"<->\",\n",
" shrinkA=False,\n",
" shrinkB=False\n",
" ),\n",
" )\n",
"\n",
" default_text_kw = dict(\n",
" va='bottom',\n",
" ha='center',\n",
" xy=((x1+x2)/2, y)\n",
" )\n",
"\n",
" an1 = ax.annotate(\"\", **{**default_arrow_kw, **arrow_kw})\n",
" an2 = ax.annotate(name, **{**default_text_kw, **text_kw})\n",
"\n",
" return [an1, an2]\n",
"\n",
"def time2phase(time, frequency=1):\n",
" return 2*np.pi*frequency*time\n",
"\n",
"def phase2time(phase, frequency=1):\n",
" return phase/(2*np.pi*frequency)\n",
"\n",
"def time_roll(a, samplerate, time_shift, *roll_args, **roll_kwargs):\n",
" \"\"\"\n",
" Like np.roll, but use samplerate and time_shift to approximate\n",
" the offset to roll.\n",
" \"\"\"\n",
" shift = np.rint(time_shift*samplerate).astype(int)\n",
" return np.roll(a, shift, *roll_args, **roll_kwargs)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
@ -278,12 +201,13 @@
"print(\"Beacon difference [ns]:\", phase2time(beacon_phase_offset, f_beacon)/ns)\n",
"print(\"Beacon difference [phase]:\", beacon_phase_offset)\n",
"print(\"Impulse offsets [ns]:\", impulses_offsets[:,0]/ns)\n",
"print(\"Time difference Impulses [ns]: {}\".format( (impulses_offsets[1,0]-impulses_offsets[0,0])/ns ))"
"print(\"Time difference Impulses [ns]: {}\".format( (impulses_offsets[1,0]-impulses_offsets[0,0])/ns ))\n",
"print(\"Time difference Impulses [T]: {}\".format( (impulses_offsets[1,0]-impulses_offsets[0,0])*f_beacon ))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
@ -443,18 +367,18 @@
"\\quad = \\Delta A + (kT) + t_\\phi\n",
"$\n",
"\n",
", where $\\Delta A < T$ and $k \\in \\mathbb{Z}$ and $t_\\phi$ is minimisable.\n",
", where $\\Delta A < T$ and $k \\in \\mathbb{Z}$ and $t_\\phi$ is minimisable by synchronising the beacons.\n",
"\n",
"Then $\\Delta t$ can be determined by iteratively summing the signals, changing $k$, and finding the $k$ belonging to the maximum of the sums."
]
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def find_best_integer_periods_sum(samplerate, f_beacon, ref_impulse, impulse, k_step=1):\n",
"def find_beacon_integer_period_sum(samplerate, f_beacon, ref_impulse, impulse, k_step=1):\n",
" max_k = int( len(ref_impulse)*f_beacon/samplerate )\n",
" ks = np.arange(-max_k/2, max_k/2, step=k_step)\n",
" \n",
@ -470,12 +394,15 @@
" if maxima[i] > maxima[best_i]:\n",
" best_i = i\n",
" \n",
" return ks[best_i], (ks, maxima)"
" return ks[best_i], (ks, maxima)\n",
"\n",
"def find_beacon_integer_period(samplerate, f_beacon, ref_impulse, impulse, k_step=1):\n",
" return find_beacon_integer_period_sum(samplerate, f_beacon, ref_impulse, impulse, k_step=k_step)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 8,
"metadata": {},
"outputs": [
{
@ -510,9 +437,9 @@
" my_impulse = time_roll(my_impulse, samplerate, -t_phi)\n",
"\n",
" # $\\Delta A$ offset\n",
" my_impulse = time_roll(my_impulse, samplerate, Delta_A)\n",
" my_impulse = time_roll(my_impulse, samplerate, +Delta_A)\n",
"\n",
"best_k, (ks, maxima) = find_best_integer_periods_sum(samplerate, f_beacon, ref_impulse, my_impulse)\n",
"best_k, (ks, maxima) = find_beacon_integer_period(samplerate, f_beacon, ref_impulse, my_impulse)\n",
"print(\"Best k: {:0g}\".format(best_k))\n",
"print(\"Maximum: {}\".format(maxima[np.where(ks == best_k)][0]))\n",
"\n",
@ -532,7 +459,7 @@
" axes[i].plot(time/ns, my_impulse, label='impulse')\n",
" axes[i].legend()\n",
"\n",
"axes[-1].set_ylabel(\"Sum\")\n",
"axes[-1].set_ylabel(\"Coherence Sum\")\n",
"\n",
"best_maximum = np.max(maxima)\n",
"axes[-1].axhline(best_maximum, alpha=0.7)\n",
@ -544,7 +471,7 @@
" summed_impulse = ref_impulse + augmented_impulses\n",
" if True or k%2 == 1:\n",
" axes[-1].plot(time/ns, summed_impulse, label='k={:.0f}'.format(k),\n",
" alpha=0.1 + 0.9*1/(1+4*abs(best_maximum-maxima[i]))\n",
" alpha=0.1 + 0.9*1/(1+2*abs(best_maximum-maxima[i]))\n",
" )\n",
" \n",
"axes[-1].legend()\n",
@ -558,12 +485,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Solve it"
"## 1. Solve it\n",
"\n",
" 1. Find $t_\\phi$\n",
" 2. Find $A_1$, $A_2$\n",
" 3. Find $B_1$, $B_2$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 9,
"metadata": {},
"outputs": [
{
@ -596,7 +527,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
@ -661,12 +592,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"##### 1.1 Beacon Phase Delays"
"##### 1.1 Beacon Phase Delay ($t_\\phi$)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 11,
"metadata": {},
"outputs": [
{
@ -679,7 +610,7 @@
],
"source": [
"beacon_phase_delays = np.array([\n",
" find_beacon_phase_delay(beacon_samplerate, f_beacon, beacons[0], beacon)\n",
" find_beacon_phase_delay(beacon_samplerate, f_beacon, ref_beacon, beacon)\n",
" for beacon in beacons\n",
"])\n",
"\n",
@ -690,7 +621,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 12,
"metadata": {},
"outputs": [
{
@ -751,7 +682,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"##### 1.2 Impulse vs beacon delays\n",
"##### 1.2 Impulse vs beacon delays ($A_1, A_2$)\n",
"\n",
"Find the delay within a single beacon period"
]