2021-12-07 19:08:54 +01:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import scipy.fft as ft\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.gridspec as gridspec\n",
"import matplotlib.ticker as tck\n",
"rng = np.random.default_rng()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# copied from 01_fourier 4988cf4f6e81b6b9510bf55a264011c37dc71872\n",
"def ft_spectrum( signal, sample_rate, fft=None, freq=None, mask_bias=False):\n",
" \"\"\"Return a FT of $signal$, with corresponding frequencies\"\"\"\n",
" n_samples = len(signal)\n",
" real_signal = np.isrealobj(signal)\n",
" \n",
" if fft is None:\n",
" if real_signal:\n",
" fft = ft.rfft\n",
" freq = ft.rfftfreq\n",
" else:\n",
" fft = ft.fft\n",
" freq = ft.fftfreq\n",
"\n",
" if freq is None:\n",
" freq = ft.fftfreq\n",
" \n",
" spectrum = fft(signal) / sample_rate\n",
" freqs = freq(n_samples, 1/sample_rate)\n",
" \n",
" if not mask_bias:\n",
" return spectrum, freqs\n",
" else:\n",
" return spectrum[1:], freqs[1:]\n",
"\n",
" \n",
"def plot_spectrum( ax, spectrum, freqs, plot_complex=False, plot_power=False, plot_amplitude=None):\n",
" \"\"\" Plot a signal's spectrum on an Axis object\"\"\"\n",
" plot_amplitude = plot_amplitude or (not plot_power and not plot_complex)\n",
" alpha = 1\n",
" \n",
" ax.set_title(\"Spectrum\")\n",
" ax.set_xlabel(\"f (Hz)\")\n",
" ylabel = \"\"\n",
" if plot_amplitude or plot_complex:\n",
" ylabel = \"Amplitude\"\n",
" if plot_power:\n",
" if ylabel:\n",
" ylabel += \"|\"\n",
" ylabel += \"Power\"\n",
" ax.set_ylabel(ylabel)\n",
"\n",
" if plot_complex:\n",
" alpha = 0.5\n",
" ax.plot(freqs, np.real(spectrum), '.-', label='Real', alpha=alpha)\n",
" ax.plot(freqs, np.imag(spectrum), '.-', label='Imag', alpha=alpha)\n",
"\n",
" if plot_power:\n",
" ax.plot(freqs, np.abs(spectrum)**2, '.-', label='Power', alpha=alpha)\n",
" \n",
" if plot_amplitude:\n",
" ax.plot(freqs, np.abs(spectrum), '.-', label='Abs', alpha=alpha)\n",
"\n",
" ax.legend()\n",
"\n",
" return ax\n",
"\n",
"\n",
"def plot_phase( ax, spectrum, freqs, ylim_epsilon=0.5):\n",
" ax.set_ylabel(\"Phase\")\n",
" ax.set_xlabel(\"f (Hz)\")\n",
"\n",
" ax.plot(freqs, np.angle(spectrum), '.-')\n",
" ax.set_ylim(-1*np.pi - ylim_epsilon, np.pi + ylim_epsilon)\n",
" \n",
" return ax\n",
"\n",
"\n",
"def plot_combined_spectrum(spectrum, freqs, \n",
" spectrum_kwargs={}, fig=None, gs=None):\n",
" \"\"\"Plot both the frequencies and phase in one figure.\"\"\"\n",
" \n",
" # configure plotting layout\n",
" if fig is None:\n",
" fig = plt.figure(figsize=(8, 16))\n",
"\n",
" if gs is None:\n",
" gs = gridspec.GridSpec(2, 1, figure=fig, height_ratios=[3,1], hspace=0)\n",
"\n",
" ax1 = fig.add_subplot(gs[:-1, -1])\n",
" ax2 = fig.add_subplot(gs[-1, -1], sharex=ax1)\n",
"\n",
" axes = np.array([ax1, ax2])\n",
" \n",
" # plot the spectrum \n",
" plot_spectrum(ax1, spectrum, freqs, **spectrum_kwargs)\n",
"\n",
" # plot the phase\n",
" plot_phase(ax2, spectrum, freqs)\n",
"\n",
" ax1.xaxis.tick_top()\n",
" [label.set_visible(False) for label in ax1.get_xticklabels()]\n",
" \n",
" return fig, axes"
]
},
2021-12-10 00:06:07 +01:00
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 3,
2021-12-10 00:06:07 +01:00
"metadata": {},
"outputs": [],
"source": [
"def phase_modulo(phase):\n",
" \"\"\"\n",
" Modulo phase such that it falls within the interval [\\pi, \\pi)\n",
" \"\"\"\n",
"\n",
" return (phase + np.pi) % (2*np.pi) - np.pi"
]
},
2021-12-07 19:08:54 +01:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Phase information in the Fourier Transform\n",
"\n",
"$$\n",
"u(t) = sin(2\\pi f t + \\varphi_t)\n",
"$$\n",
"\n",
"Define $f_\\mathrm{max}$ as the frequency with the highest power in the FT (it should be close to $f$).\n",
"Then $\\varphi_f$ is its associated phase."
]
},
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 4,
2021-12-07 19:08:54 +01:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Required signal length is: 0.005s\n",
"Required number of samples: 50.0\n"
]
}
],
"source": [
"sample_rate = 1/1e-4 # Hz\n",
"f = 200 # Hz\n",
"required_N_samples = sample_rate/f\n",
"\n",
"signal_func = lambda phase: np.sin(phase)\n",
"\n",
2021-12-09 20:24:59 +01:00
"# set signal_func to exp(i*phi)\n",
"if False:\n",
" signal_func = lambda phase: np.exp(1j*phase)\n",
"\n",
2021-12-07 19:08:54 +01:00
"print(\"Required signal length is: {}s\".format(1/f))\n",
"print(\"Required number of samples: {}\".format(required_N_samples))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\varphi_f$ vs $f_\\mathrm{max}$ for differing $\\Delta f$"
]
},
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 5,
2021-12-07 19:08:54 +01:00
"metadata": {},
"outputs": [
{
"data": {
2021-12-10 00:06:07 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAIeCAYAAACV/nMhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5wkdZ3/8denqsOEzYkNLCw5KUFJEswKmODOU9EzYPyZwxkwnJ7nGc8znPE8xRNREZQgKCiggKKgLriEXdK6hF02p9lJnao+vz++NTM9szOzs2l6Zvb9fDz6Md3V1dWfqg7Tn/p+v5+vuTsiIiIiIiIiE0XU6ABERERERERE9iQluiIiIiIiIjKhKNEVERERERGRCUWJroiIiIiIiEwoSnRFRERERERkQlGiKyIiIiIiIhOKEl0RmdDMrMPMDm50HIMxs+eb2dV7ads/MLNP7+FtXmBmt+3Jbe4JZvZRM/veGIijaGYPmNmcRseyO8zsdDN7OPvsnNfoePYEM/sfM/t4o+MYipk908xW7a31Bzz2TDN7cFceO8T2xuT3goiIEl0RaRgze9TMKmY2a8DyJWbmZrZod5/D3Se5+4rd3c5e8lng840OYrxz98+6+5vGQBxl4PvAhUOtY2afNLNqlkRuNbM/mdnTRi/KEfkU8I3ss7NXTsSMNnd/q7v/R6Pj2FXZ9+Ghe2Jb7v4Hdz9iT2xLRGQsU6IrIo32CPDKnhtm9mSguXHhjA4zOwmY6u53NDoW2aN+ArzOzIrDrHOZu08CZgO3AVeamY1KdHXMLDfEXQcCS4d4jJmZfjtMUMO8J0RExh39sxKRRrsEeG3d7dcBP6xfwcxeaGZ/M7NtZrbSzD5Zd98rzGyFmU3Jbp9jZmvNbHZ2u7clJOvO+y0zuz5rUfujmc01s6+a2Zas2+kJddvu14pS3x24p+ugmX3IzNab2RozO8/MXmBmD5nZZjP76DD7fQ5w64D9/O9s/7aZ2Z1mdmbdfZ80s8vN7Idm1m5mS83sxLr7TzCzu7L7LgOahnpiMzvUzG41szYz25itj5ktyvY5V7fuLWb2pv4Pt69nj33AzJ4zzPM8amYfNLN7zKzTzC4ys/2y499uZjeZ2fS69X+WvXZtZvZ7MzsmW17IWvnfld2Os9fuE3XH5kcD9uH12bHcYmZvNbOTsji2mtk3BhzXH9Xd7ncMsv3/dNby2mFm15rZTDP7cfY6/dXqeh64+ypgC3DqUMelbt0qcDEwF5hpZpGZ/auZPZa9p35oZlOzOC42s/dn1xdkMb697vXcbBaSZTN7UXa8elqMjx3wmlxoZvcAnTYgsTGzvwMHA9dm+1vMjsFnzOyPQBdwsJlNzV7PNWb2RHaM4rrX57+y99YKM3vHgGP6qJk9d5jX4NQs7q1mdreZPbPuvlvM7D+y17/dzG6wuh4hZnZG3WNXmtkF2fJ+Xfl3cIwuzPap3cweHOo9bsN/L/W8j15nZo9nx+Jjdfc3ZzFtMbNlwEmDv0vAzH6fXb07e01eUXff+63v++f1dcuL2WvwuJmts9B1uzm7r1+35x29J7J13Mzenb2eG83sizbghEf2fFvM7BEzO6du+evN7P7seK4ws/9Xd98sM/tl9jpsNrM/9GzXzOab2RVmtiHb5ruHOkYiIoNRoisijXYHMMXMjsp+KL8C+NGAdToJyfA04IXA2ywbO+julwG3A18zs5nARcCb3H3DEM/3cuBfgVlAOXvsXdntnwNf3onY5xISygXAJ4DvAq8GngqcCXzChh4f/GRg4Di5vwLHAzMILYM/M7P6hPUlwE8Jx+Ea4BsQEkHgasJJgxnAz4CXDhP3fwA3ANOB/YGv73hXe50CrCAcr38jtEbOGGb9lwLPAw4HXgxcD3w0e3wE1P94vR44DJhDeE1+DODuFcJx/ZSZHQV8GIiBz+wgzsMI76evAh8DngscA7zczJ4x4j2G84HXEF7nQwjvmf8jHOv7Cceh3v3AcTvaqIVW3wuAVe6+Mbt+AfAsQrI5iew1JpwUeWZ2/RmE16BnH54O/MHd3cyeQug+/f+AmcB3gGusfwvzKwmfo2nuXquPyd0PAR4HXpx1XS5nd70GeAswGXiMkKDXgEOBE4DnAz0nRN4MvChbfiLwTzs6FnXHZAHwK+DThOP7AeAKy05cZV4FvJ7wPilk62BmBxDeQ18ntJYfDywZ5DmGPEZmdgTwTuAkd58MnAU8OkS4Q34v1TkDOAJ4DuH74Khs+b8R3kuHZM/xuqGOibs/Pbt6XPaaXJbdngtMJbwv3wh80/pOHH2B8Jk7nvAa9XxHDWXI90SdfyC8nk8BzgXeUHffKYTvs1nAfwIXmfX2UlhPeD9MIbxuX8leA4D3A6sIr9d+hO8Gz5Lda4G7s9ifA7zXzM4aZh9ERPpRoisiY0FPq+7zgAeAJ+rvdPdb3P1ed0/d/R7gUvp+5AO8A3g2cAtwrbv/cpjnusrd73T3EnAVUHL3H7p7AlxG+HE+UlXgM1nL3E8JP/L+293b3X0pofvnsUM8dhrQPmA/f+Tum9y95u5fAoqEH8k9bnP367JYL6EvmToVyANfdfequ/+ckDQPF/eBwHx3L7n7zhSSWV/3PJcRfty+cJj1v+7u69z9CeAPwJ/d/W9ZAnUVdcfb3b+fHbsy8EngOMtaNN39PkLycxUhsXlNdhyG8h/Zvt1ASEgudff1dXHszOv8f+7+d3dvIyRSf3f3m7KE4GeDbKud8PoO5eVmthVYSTgp0pMc/TPwZXdf4e4dwEeA87MWtluBM7ME4OmEZOL07HHPoK93wJuB77j7n909cfeLCSd06luYv+buK929eyeOwQ/cfWm2zzMIPRLe6+6d7r4e+ArhhACEk0lfzZ5jM/C5nXieVwPXZe/z1N1vBBYDL6hb5//c/aEs/ssJyRyE43eTu1+avT83uft2iS7DH6OE8Lk72szy7v6ou/99sEBH8L0E8O/u3u3udxOStp7P7MsJ3x2b3X0l8LWdOEY9qsCnsn29DugAjsgSzDcD78u2306oB3D+MNsayXviC9n2HiecPHpl3X2Puft3s8/kxcA8QuKKu/8q+/y4u99KOMnW01ulmq17YLYff3B3J7Rwz3b3T7l7xUOdhe/uYB9ERPpRoisiY8ElhFaaCxjQbRnAzE4xs5uzLmxtwFsJSSUA7r6VkHA8CfjSDp5rXd317kFuT9qJuDfVJVs9PxBHur0thNaxXlk3xPstdN3dSmitqS/UtbbuehfQlCVB84Ensh+IPR4bJu4PAQb8xUIX6DcMs+5Agz3P/GHWH9HxttDd9fNm9ncz20ZfK1r9/l8MLCIkQg/vIM49+Trv7LYmA1uH2d7l7j7N3ee4+7Pd/c5s+Xz6v26PATlgvyzZ6iAkdWcCvwRWZy2Q9YnugcD7s66gW7P30UL6v0Yrh9/dQdU/5kDCiZU1dc/xHUILa89+1K8/3HtxoAOBlw2I/wxCMtRj4Oeg5/gvBAZNSgd5jkGPkbsvB95LONGy3sx+amaDvr939L20g1h35xj12DSg9bVn+7OBFuDOuv37dbZ8KCN5TwyMt/649O6nu3dlV3s+2+eY2R1Z1+SthJMWPcfpi8By4IasW/OHs+UHAvMHvEYfJUueRURGQomuiDScuz9GKEr1AuDKQVb5CaGr7kJ3nwr8DyFRA8DMjid0o7uUXWsZGUoX4Qdjj7l7cNv3ELoWAmHKD0K13pcD0919GtBG3X4OYw2woK6rIMABQ63s7mvd/c3uPp/QffNbFsYid2arDLfPgz3P6hHEuCOvInSHfC4hwV+ULa9/rm8REryzzOyMPfCcEPZ5T7/GRxFa73bWasIP/B4HELoH9yTWtxK6AReylulbCT0hptPXRXcloaVwWt2lxd0vrdtu/YmKkap/zEpCC+isuueY4u7HZPevISSO9ftRb7hjvhK4ZED8re4+kurkKwldgUey3pDHyN1/4u5nEF4LJ3QDHsyw30s7sKNjtDs2Ek7AHFO3f1M9FEAbykjeEwPj3eHnPusyfwXwX4QTNtOA68iOU9aD4/3ufjB
2021-12-07 19:08:54 +01:00
"text/plain": [
2021-12-09 20:24:59 +01:00
"<Figure size 1152x576 with 4 Axes>"
2021-12-07 19:08:54 +01:00
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# \\phi_f vs f_max for differing \\Delta f\n",
"N_delta_f = 200\n",
"\n",
2021-12-09 20:24:59 +01:00
"plot_submax = True\n",
"\n",
2021-12-07 19:08:54 +01:00
"\n",
"Ns_samples = required_N_samples//1 + np.arange(0, N_delta_f)\n",
"phi_f = np.empty(N_delta_f)\n",
"f_max = np.empty(N_delta_f)\n",
"\n",
2021-12-09 20:24:59 +01:00
"if plot_submax:\n",
" phi_f_sub = np.empty(N_delta_f)\n",
" f_submax = np.empty(N_delta_f)\n",
2021-12-07 19:08:54 +01:00
"\n",
"for i, N_sample in enumerate(Ns_samples):\n",
" time = np.arange(N_sample) / sample_rate\n",
" \n",
" fft, freqs = ft_spectrum(signal_func(2*np.pi*f*time), sample_rate)\n",
" \n",
" fft_power = np.abs(fft)**2\n",
" id_max = np.argmax(fft_power)\n",
" \n",
" phi_f[i] = np.angle(fft[id_max])\n",
" f_max[i] = freqs[id_max]\n",
" \n",
2021-12-09 20:24:59 +01:00
" if plot_submax:\n",
" fft_power[id_max] = 0\n",
" id_submax = np.argmax( fft_power )\n",
" \n",
" phi_f_sub[i] = np.angle(fft[id_submax])\n",
" f_submax[i] = freqs[id_submax]\n",
"\n",
"if plot_submax:\n",
" fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2, figsize=(16,8), sharex=\"col\")\n",
" fig.suptitle(\"Maximum (and sub maximum) Power frequencies and their phase\")\n",
"else:\n",
" fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,4))\n",
" fig.suptitle(\"Maximum Power frequencies and their phase\")\n",
" \n",
2021-12-07 19:08:54 +01:00
" \n",
"# Maximum values\n",
"ax1.set_xlabel('$f_\\\\mathrm{max}$')\n",
"ax1.set_ylabel('$\\\\varphi_f$')\n",
"ax1.plot(f_max, phi_f, '--', alpha=0.1)\n",
"sc = ax1.scatter(f_max, phi_f, c=Ns_samples, cmap='Spectral')\n",
"ax1.axvline(f, color='r', alpha=0.5, label=\"Signal frequency\")\n",
2021-12-09 20:24:59 +01:00
"for hline in [0, -np.pi/2]:\n",
" ax1.axhline(hline, color='k', alpha=0.5)\n",
2021-12-07 19:08:54 +01:00
"\n",
"ax2.set_xlabel('$N_\\\\mathrm{samples}$')\n",
"ax2.set_ylabel('$f_\\\\mathrm{max}$')\n",
"ax2.scatter(Ns_samples, f_max, c=Ns_samples, cmap='Spectral')\n",
"ax2.axhline(f, color='r', alpha=0.5, label=\"Signal frequency\")\n",
"\n",
2021-12-09 20:24:59 +01:00
"# SubMaximum values\n",
"if plot_submax:\n",
" \n",
" # filter submax frequencies above twice the frequency\n",
" if True:\n",
" idx_submax = np.argwhere(np.abs(f_submax) < 2*f)\n",
" \n",
" f_submax = f_submax[idx_submax]\n",
" phi_f_sub = phi_f_sub[idx_submax]\n",
" \n",
" Ns_samples_submax = Ns_samples[idx_submax]\n",
" else:\n",
" Ns_samples_submax = Ns_samples\n",
" \n",
" \n",
" ax3.set_xlabel('$f_\\\\mathrm{submax}$')\n",
" ax3.set_ylabel('$\\\\varphi_{f\\_sub}$')\n",
" ax3.plot(f_submax, phi_f_sub, '--', alpha=0.1)\n",
" sc = ax3.scatter(f_submax, phi_f_sub, c=Ns_samples_submax, cmap='Spectral')\n",
" ax3.axvline(f, color='r', alpha=0.5, label=\"Signal frequency\")\n",
" for hline in [0, -np.pi/2]:\n",
" ax3.axhline(hline, color='k', alpha=0.5)\n",
"\n",
" ax4.set_xlabel('$N_\\\\mathrm{samples}$')\n",
" ax4.set_ylabel('$f_\\\\mathrm{submax}$')\n",
" ax4.scatter(Ns_samples_submax, f_submax, c=Ns_samples_submax, cmap='Spectral')\n",
" ax4.axhline(f, color='r', alpha=0.5, label=\"Signal frequency\")\n",
"\n",
"\n",
"\n",
2021-12-10 00:06:07 +01:00
"if False:\n",
2021-12-09 20:24:59 +01:00
" res = 50\n",
2021-12-10 00:06:07 +01:00
" ax1.set_xlim(f-res, f+res)\n",
" ax2.set_ylim(f-res, f+res)\n",
" if plot_submax:\n",
" ax4.set_ylim(f-res, f+res)\n",
2021-12-07 19:08:54 +01:00
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2021-12-09 20:24:59 +01:00
"## $\\varphi_f$ vs $\\varphi_t$ and the effect of $f/f_\\mathrm{sample}$"
2021-12-07 19:08:54 +01:00
]
},
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 6,
2021-12-07 19:08:54 +01:00
"metadata": {},
"outputs": [
{
"data": {
2021-12-10 17:39:00 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAEuCAYAAAC3Tv7YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xc5ZXw8d+5M6NRL65ykwsu2GCwTY9NN8EmhJoECAnZBEKym92Q8O6m74Zkw27YN3mXBJJNCCyQQiCUUIITmrFNMTbuTa6ybFm2bMnqberz/nHvyCNppBnVmbHPN5/5eGZuO/femTBH5ylijEEppZRSSimllEp1VrIDUEoppZRSSimlEqEJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKKaWUUkqptKAJrFJKqQERkSdE5EdJOO52EblsCPZbLiKLB3u/g0VEZonIRhFpEpGvJjuedKLXTiml0p872QEopVQ6EpFyYCwQinp7pjHmcHIiOvUYY85IdgxJ8g1ghTFmfrIDSUN67ZRSKs1pBVYppfrv48aY3KhHp+RVRPSPhGooTAa2JzuIVNbLd6/Ha6ffV6WUSg+awCql1CBymp9+U0S2AC0i4haR8SLyvIhUi8j+rk0XRWS+iGxwmjU+IyJPR5rkiogRkelR63Zqrtvbvp1Y/llEtohIg7PvzKjlk0TkBWfb4yLysPP+v4jI811ifEhEHowVLxC9z9kiskJE6p0mvtfFuD7/4sTUIiKPichYEfmrs783RaTIWfdbIrLPeX+HiNwYY1+LEznXGNt929lnnYg8HmPdeb1ctx7jcu59pbNsl4hcGe8+xYivx2soIsuBy4GHRaRZRGb2cH4DvsYicpqI1IrIgqhzqJEemm3Hu649nZeIfF5EXolab6+I/CnqdYWIzIt3HSXGd69LfN2uXaxt4t2rrp9/SfHva4z79Bknhjud+/SaiOTGWlcppVKSMUYf+tCHPvTRxwdQDizu4f1NwCQgC/sPheuBfwMygGlAGXC1s34GcAD4OuABPgEEgB85yw0wPWr/T0Qti7fvcmAtMB4YAZQCX3aWuYDNwH8DOdhJ6CJn2TigBSh0XruBY8A5vcXrvN4LfMdZ7wqgCZjV5fp8gN38eoKz3w3AfMALLAe+76z7SSd2C7jFiWlcrHvQ27n2cI+2OfdoBPBe5Jomsq+e4gJmARXAeGe9KcBp8e5Tl9gSuYYrgLvifDYH6xp/0Tn/bOA14CdxjhvzuvZ2Xs71qHdiGIf9+ap0tpsG1DnLEvm8d3z3eoix07Xruk0Cx0ir72uM878WqATmOtd1LrALuDvZ/5+qD33oQx+JPrQCq5RS/feiU02qF5EXo97/uTGmwhjTBpwHjDbG/NAY4zfGlAG/AW511r0Q+4fwg8aYgDHmOeDDBI8fb9+RWA4bY2qBV4B5zvvnY/9Q/hdjTIsxpt0Y8y6AMeYIsAo7uQFYAtQYY9bHifdCIBf4sRPPcuAvwG1d4n7IGHPUGFMJvAOsMcZsNMb4gD9jJ1oYY551Yg8bY54B9jhx96Snc43lYece1QL3x4ixx331ElcIO0GcIyIeY0y5MWYfid2niESvYTyDco2NMb9x3luDnSh9N85xe7quPZ6Xcz2asK/xpdiJcqWInO68fscYEybxz3vku5eok/n72tUPsPsBW0CLMWYr8D52Ao9TqX9XRJYleE5KKTXstL+HUkr13w3GmDdjvF8R9XwyMF5E6qPec2EnFWD/KK00xpio5QcSPH68fQNURT1vdY4H9g/WA8aYYA/7fhL4e+wf2J8BfpdAvOOBCifZiF42ocu+j0Y9b4vxOhdARO4A7sWuZOK8P6qHeKHnc40l+h4diLFuj/vqKS5jzF4R+RpwH3CGiLzmrJfIfYpI9BrGM5jX+DfAy9hVOp+z3e3Ar53l7xhjljrPe7qu8c5rJXAZMN15Xo+dvF7kvIbErmP08RN1Mn9fO4jIGGABdmJ8PXblGWA0doUe7D84POgk5koplZK0AquUUoMv+sdtBbDfGFMY9cgzxlzjLD8CTBARidqmJOp5K3bzzYjiPuy7NxVASdd+glFeBM4SkTOxmx3+IYF4DwOTRMTqsqwygXg6EZHJ2D/G/xEYaYwpxG6eKr1umLhJUc9LsGMfcFzGmKeMMYuwkxUDPEDf7tOgXcOBnouzTi7wIPAYcJ+IjHDO8w/mxOBlS6N229N1jXdekQT2Yuf5SuwE9lJOJLCJXMfo716iTubva7RRQBC72n02sFlEcoBLgGUicjF2E++vi8i1CcSklFJJoQmsUkoNrbVAozNQTJaIuETkTBE5z1m+GvtH5VedAWRuonMz2U3Ap53tlmD/oE903/HiOgL8WERyRCRTRBZGFhpj2oHngKeAtcaYgwnEuwa7L943RMQj9mA/HweeTuxSdZKDnVhUgz3QD3BmP/bTk6+IyEQnIfsO8MxA4xJ7jtErRMQLtGNXOkP07T4N5jXs97lE+Rmw3hhzF/Aq8Ks4++zpusY7r5XYAyxlGWMOYVcllwAjgY3OOgP5vCfqZPu+RisHfNjV17Ox+9/+L/B7Y8w+Y8w7wFbgYmPMXxKISSmlkkITWKWUGkLGmBD2D/V5wH6gBngUKHCW+4GbgL/DHlTlFuCFqF3c42xfD9yOXWlJaN8JxjUdOAgcco4d7UnsQV5+F7Vdj/E6y64Dljqx/BK4wxizM148MeLbAfwUO2E46sTxXl/304ungNexf8SXYQ9CNdC4vMCPsc+9ChgDfKcv92kwr+EAzwURuR47ifyy89a9wAKn+XBPYl7XeOdljNkNNOM0pzXGNDrbv+dcvwF93hN1sn1fu+yjFbgLeAi4Evve7Ae+BhCpjndp5q2UUilHOnfjUEoplWwi8gRwyBjzvSTHUQLsBIqdhOKkICLl2CPRxuq/rPrpVL2u6fZ9FZFxwD4gJ7ovr4jMwB4k6u4hD1YppQZAK7BKKaW6caox9wJPn0zJq1Inoz5+X08H9pjuFYyzsafqUUqplKajECullOrEGdjlKPboqkuSHI5Sqhf9+L7OAnZ3fVNHHlZKpQttQqyUUkoppZRSKi1oE2KllFJKKaWUUmlBE1illFJKKaWUUmlBE1illFJKKaWUUmlBE1illFJKKaWUUmlBE1illFJKKaWUUmlBE1illDoFich2EblsGI7zhIj8aKiPE+O4Q3Z+IlIuIouHYt+pIt71S9VrMNSfaxEZLSJviEidiDzW1+VKKaUGTueBVUqpk4yINEe9zAZ8QMh5/SVjzB+MMWcMf2TD52Q/v6HW9fqJSDlwlzHmzeRElJhhuO/fBvYYY64CEJHxwFpjzMRYy5VSSg0+rcAqpdRJxhiTG3kAB4GPR733h2THp1R/iEgq/NF9MfBs1OtrgL/1slwppdQg0wRWKaVSiIh8S0Se6/Lez0Tk587zb4pIpYg0icguEbmyn8fpaALqPP8XEdkiIi0i8piIjBWRvzrHeVNEiqK2HS8iz4tItYjsF5GvRi2bLyIbnO2eATK7HHe2iKwQkXqnued1XWLqSxzfEpF9zrIdInJjrPOLev3Pzr4bROQZEekUW5d1v+3ss05EHo+x7rxY++otJmd5zPvX2zXtIcbPOOvf6cT4mojkxtnm8yLyStTrvSLyp6jXFSIyr+v1E5HfASXAKyLSLCLf6O0axDhuvM90vPv4TRHZArQ4n4/nu+zrIRF5MEbcvd5zEVkgIhud4z7rLI/Z3F1EMkSkAZjrXIetzqJrgGW9LFdKKTXYjDH60Ic+9KGPFHkAk4FWIN957QKOABcCs4AKYLyzbApwWpz9lQOLe3vfef4BMBaYABwDNgDzAS+wHPi+s64FrAf+DcgApgFlwNXO6wPA1wEP8AkgAPzI2dYD7AW+46x7BdAEzOprHM76nwTGOzHdArQA42Kdt/N6rbP+CKAU+HIv12wbMMlZ973IOcTbV5yYYt6/3q5pD/FdC1RiJ0t1zr+7gLvjfBa
2021-12-07 19:08:54 +01:00
"text/plain": [
"<Figure size 1152x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
2021-12-09 20:24:59 +01:00
"sample_rate = 1/1e-4 # Hz\n",
"\n",
"phase_offsets = np.linspace(-np.pi, np.pi, 500, endpoint=True)# rad\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-09 20:24:59 +01:00
"frequencies = sample_rate * np.array([0.5, 0.49, 0.45, 0.3, 0.28, 0.25, 0.05, 0.001])\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-09 20:24:59 +01:00
"signal_func = lambda phase: np.sin(phase)\n",
2021-12-07 19:08:54 +01:00
"\n",
"\n",
2021-12-09 20:24:59 +01:00
"# Precreate figure and axis\n",
2021-12-07 19:08:54 +01:00
"fig, (ax1) = plt.subplots(1,1, figsize=(16,4))\n",
"\n",
2021-12-09 20:24:59 +01:00
"for f in frequencies[::-1]:\n",
" \n",
" required_N_samples = sample_rate/f\n",
"\n",
" phi_f = np.empty_like(phase_offsets)\n",
"\n",
" time = np.arange(required_N_samples) / sample_rate\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-09 20:24:59 +01:00
" for i, offset in enumerate(phase_offsets):\n",
"\n",
" fft, freqs = ft_spectrum(signal_func(2*np.pi*f*time + offset), sample_rate) \n",
" id_max = np.argmax(np.abs(fft)**2)\n",
"\n",
" phi_f[i] = np.angle(fft[id_max])\n",
2021-12-07 19:08:54 +01:00
" \n",
" \n",
2021-12-09 20:24:59 +01:00
" ax1.plot(phase_offsets, phi_f, '.--', label=\"$f/f_\\\\mathrm{{sample}} = {}$\".format(f/sample_rate))\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-09 20:24:59 +01:00
" if True:\n",
" id_phi_f_min = np.argmin(phi_f)\n",
" ylocation = (np.max(phi_f) + np.min(phi_f)) /2\n",
" ax1.text(phase_offsets[id_phi_f_min], ylocation, \"${:.2g}\\\\pi$\".format(phase_offsets[id_phi_f_min]/np.pi), horizontalalignment='center')\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-09 20:24:59 +01:00
" \n",
2021-12-10 17:39:00 +01:00
"ax1.set_title(\"Frequencydomain phase of max-power frequency $\\\\varphi_f$ \\n vs Timedomain phase $\\\\varphi_t$ with varying $f / f_\\\\mathrm{{sample}}$\")\n",
2021-12-07 19:08:54 +01:00
"ax1.set_xlabel('$\\\\varphi_t$')\n",
"ax1.set_ylabel('$\\\\varphi_f$')\n",
2021-12-09 20:24:59 +01:00
"ax1.legend(loc='lower right')\n",
2021-12-07 19:08:54 +01:00
"\n",
"# grid lines\n",
2021-12-09 20:24:59 +01:00
"## vertical lines\n",
2021-12-07 19:08:54 +01:00
"vlines = [\n",
" (-np.pi, r'$-\\pi$'),\n",
" (-np.pi/np.sqrt(2), r'$\\frac{-\\pi}{\\sqrt{2}}$'),\n",
" (-np.pi/np.sqrt(3), r'$\\frac{-\\pi}{\\sqrt{3}}$'),\n",
" (-np.pi/2, r'$\\frac{-\\pi}{2}$'),\n",
" (np.pi, r'$\\pi$'),\n",
"]\n",
"\n",
"xtrans = ax1.get_xaxis_transform()\n",
"ax1.axhline(0, alpha=0.1, color='k')\n",
"for location, label in vlines:\n",
" ax1.axvline(location, alpha=0.1, color='k')\n",
" ax1.text(location, -0.06, label, transform=xtrans, horizontalalignment='center')\n",
"\n",
2021-12-09 20:24:59 +01:00
"## horizontal lines\n",
2021-12-07 19:08:54 +01:00
"hlines = [\n",
2021-12-09 20:24:59 +01:00
" (1, ''),\n",
" (-2, ''),\n",
2021-12-07 19:08:54 +01:00
" (-np.pi/2, r'$\\frac{-\\pi}{2}$'),\n",
" (np.pi/2, r'$\\frac{\\pi}{2}$'),\n",
"]\n",
"\n",
"ytrans = ax1.get_yaxis_transform()\n",
"ax1.axvline(0, alpha=0.1, color='k')\n",
"for location, label in hlines:\n",
" ax1.axhline(location, alpha=0.1, color='k')\n",
" ax1.text(-0.03, location, label, transform=ytrans, verticalalignment='center')\n",
" \n",
"plt.show()"
]
},
2021-12-09 20:24:59 +01:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $f_\\mathrm{sample} \\geq 3f$ the relationship between $\\varphi_t$ and $\\varphi_f$ is (almost) linear.\n",
"\n",
"From $f_\\mathrm{sample} \\geq 4f$ onwards, this relationship is stable with\n",
"\n",
"$$\n",
"\\varphi_f = \\varphi_t - \\frac{\\pi}{2} \\delta_\\mathrm{sin}\n",
",\n",
"$$\n",
"where $\\delta_\\mathrm{sin}$ is 1 if the signal was a sine, and 0 for a cosine."
]
},
2021-12-07 19:08:54 +01:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
2021-12-10 17:39:00 +01:00
"# reconstruct phase from off-frequency ft\n",
"\n",
"Require atleast $f_\\mathrm{sample} \\geq 4 f$"
2021-12-07 19:08:54 +01:00
]
},
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 7,
2021-12-07 19:08:54 +01:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2021-12-14 16:17:42 +01:00
"Required signal length is: 0.01s\n",
"Required number of samples: 100.0\n",
2021-12-10 00:06:07 +01:00
"Phase to be retrieved: -1.5707963267948966\n"
2021-12-07 19:08:54 +01:00
]
}
],
"source": [
"sample_rate = 1/1e-4 # Hz\n",
2021-12-10 17:39:00 +01:00
"f = 100 # Hz\n",
2021-12-07 19:08:54 +01:00
"required_N_samples = sample_rate/f\n",
"\n",
2021-12-10 00:06:07 +01:00
"phase_to_retrieve = phase_modulo(-np.pi/2)\n",
"signal_func = lambda phase: np.cos(phase + phase_to_retrieve)\n",
2021-12-07 19:08:54 +01:00
"\n",
"print(\"Required signal length is: {}s\".format(1/f))\n",
2021-12-10 00:06:07 +01:00
"print(\"Required number of samples: {}\".format(required_N_samples))\n",
"print(\"Phase to be retrieved: {}\".format(phase_to_retrieve))"
2021-12-07 19:08:54 +01:00
]
},
{
"cell_type": "code",
2021-12-10 17:39:00 +01:00
"execution_count": 8,
2021-12-07 19:08:54 +01:00
"metadata": {},
"outputs": [
{
"data": {
2021-12-14 16:17:42 +01:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAEPCAYAAABsnLiCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZBd133g9+85d31b9+sd3Y2lsRPcAO7UZkmUJUsyJcqiLdtJPHEyHse2PNZMFkUqR65Y4xkVkyrnDyWlklhMxS473u2xpSjW2LFoi9pIcANAgMTeABrovfvtdz/54zVfgxTNxtJEN5q/T9WtOu+90+f+7q3XwP312ZQxBiGEEEIIIYQQYqPQax2AEEIIIYQQQgixmiTRFUIIIYQQQgixoUiiK4QQQgghhBBiQ5FEVwghhBBCCCHEhiKJrhBCCCGEEEKIDUUSXSGEEEIIIYQQG8q6SHSVUh9WSr2ilDqplPrcG3z+i0qpGaXUC0vHL61FnEIIIVafUmqLUurbSqljSqmXlFKfWeuYhBBCCHFzU2u9j65SygKOAx8ELgDPAD9vjDl6WZ1fBO41xvz6mgQphBDiLaOUGgaGjTHPKaVKwLPAJy7/f0AIIYQQ4mqshx7d+4GTxpjTxpgI+GPgkTWOSQghxA1ijLlkjHluqVwDjgGjaxuVEEIIIW5m6yHRHQXOX/b6Am/8gPOoUuqQUurPlVJbbkxoQgghbiSl1BhwF/DDtY1ECCGEEDcze60DANQbvPf68dRfB/7IGBMqpX4F+D3goTdsTKlfBn4ZoFAo3HPLLbesZqzideafP0mWZQBkRR9ViskyRb4VkXl95AZ6cHOr8zVbOHyaNEoAaAx1U0gDkrwiSRuU/UGK5b7XxvbCSbI0AwWqv0BASsk4JIstTJpijMF0WUS+RzGzSBebZMYQDzmQaqxqDEUXD4uW08QYGwcL4gTjZiSpBQoMGjKDURpIUWlKmk+wjIU90b43uXKR0vbhVbkPQlyLZ599dtYYM7CabY4O3WnCsLZivbnK2ZeA4LK3vmaM+drr6ymlisBfAP/GGFNdtUDFutLf32/GxsbWOgwhhBAbwJs936yHRPcCcHkP7Wbg4uUVjDFzl718HHjsn2ts6eHpawD33nuvOXjw4OpFKn7Etz7yP3LxW+17HG3aDD+R40j/IGrqZR681OIT/+Ex8rfeuirn+vbPfpGzf/aPAJz/6Y9Q5CXCW7up157hvujdfOILr/1afOM9n2Hmu0cwQOHRPXw77/KJ2iD175ygduICZIbGlgLTH9vKtkYP2Q/G0QcvEv/cbny9k+wfjzP5wQS9oNnkzzCvNjPYH5I7ljC+p0o97sdKmoRWARNGqCxPlITYUUBz/wVKYZnBz8UAfPJv/xeG77ttVe6DENdCKTW+2m0GUZ2PPvTvVqz3B3/1C4Ex5t43q6OUcmgnuX9ojPnLVQpRrENjY2PI/81CCCFWw5s936yHocvPALuVUtuVUi7wc8DfXF5haaGSV32c9vwtsQ7s+a8/0imXpqcYnN6JScc4WLqX6XCOqcNPr9q5Rj98X6c8duQQ9dEEa2GeeypF3BfO/Uj9nn1bgfaQgZcU9Pc1mNQLVL0UZbW/+rnzUOwJuOjWqPZrNAZttbCw0GUP7SiMgsS28JIGSWaR2KC1htQmUS5pZuPrAq6b4tkZRqXccjyj/2LaieXI7/zRqt0HIdYTo9WKx0qUUgp4AjhmjPndtzxoIYQQQmx4a57oGmMS4NeBb9FOYP/UGPOSUuqLSqmPL1X7jaUtJ14EfgP4xbWJVrze8Afu7gw+j+oRamKaPY0mOa+fS5s8jo3PcuHY3Js3coV679yJsi0A3OMzWHmwI0PFyROm07x+BfHhh+7qlJuVEj3FjItOlXoxQlntdkytgUpBGUPc47br2mBlBtWbIxfa7GEXl6w+asZQ1XAyl1Je2IvXikjiEfzAoxBZEKcoKySNA3ZPJIzMR53z2wV/Ve6BEOuKgsxSKx5X4F3ALwAPXbaN3Eff2uCFEEIIsZGth6HLGGO+CXzzde/91mXlzwOfv9FxiZX5vV14A2XC6UWMUVTnjrGtrrl7V45+FXI0nGFgvMLInh60dX1/V+ndv7OTVIdzAfmFjOH57Xxzs0/l1he448Rhdu25s1O/sHUQZWlMmtEzF6IDH6UT3P5uoD39zwLiyLA9GSIupxhO0Ui2kK9XqGEopAkjDHIkrRDGE9jKULe78RKF49TIMUTBqqKUwkpsWpmPtlOMAnXZVPPhD95zXdcuxHpkgOwKemxXbMeYp3jj9RqEEEIIIa7Jmvfoipvfto+/q1OeNfMEM+fYFjYwuc3UGi8yO3OJykzrus+jbYv8SD8AaawwNYtiK6VpxpgqlDl55MXX1O/evRmTtZPNzfOzNIImBxq93FG+B3VZ568iw88MyvEwQDluEZoYJ+cT6hgMeCToqIFlUrRtk9l1osxgewlRFlNPI6K0RZopdJLwgl/gJa/QOcfLX/6r675+IdYd1U50VzqEEEIIIW40SXTFddv36eVtjwvnXeYal+ibiWnSx47JKrmXjjJ97NKqnGvTe/cDoFA4Zy3mrXk2ZS1KvuHYke+/pq4/UMbynfaL4/PYClKrhXIc7Fy+U88LNZEJ2r2/nk0Y15nKhUx2R1RNF4vZJJkOmL97msGpLRR7Qypeg9m4RMsz1COPRuKTRC46TYizkIrqopH1dM4x8ODqLMglxHpilCJ19IqHEEIIIcSNJk8g4rr17t/ZWdxJxy4zw0XOLdQppAsU/ALO3Almj02QxOkKLa1s03uXhyb3ndD8Q7nI7eok7/LOM1898iP1c8PtLYeiekC55nLajWnpOna5CEs9TV0zIWdMg4FGjOu4LEQ5KoWAuhPhJj5BFmIiH+1oFClZpklSF4WDjjW5rIBtwMXGtxQ6l8PYBqWXrzeYrVz3tQuxHkmPrhBCCCHWI0l0xarov28vAHaasuX4FIkxFLOE8a4hJuqvkC0eJ4lWIdF9z53Li1+ZsyS5QfLzdbTVQzPfpNpcfE39vrt2t+PKNLlmyAuFbgJTwel1YWnxqrRRpelYgMEueOSOVFAqQRuL3HSVUMXYkWL0XJGZQoDTykh0Rr6/iRfXMG4T5YWYfIKDJrEtmgdqWHvmO3HMP3f8uq9diPXGKMi0XvEQQgghhLjR5AlErIo9v7S8QOqlaC+u06K4cBf/UPgQf9ndRdw6iJ6bvO7zdO0aRbvt4cjJuZjdbpVqEqNjn9HEo3n2tdsM9e7f2SkHqoF2uwlME7u3wKtrRWU/NDj5Bmf9BpVtBfonZthaLzGqRliIfU6Vpgn8WcamMloKbEIiP8TOpyReC1NQGGNhHE1qDCoArQzqsonAYz/3/uu+diHWn5W3FrqS7YWEEEIIsb585clTfPHrRzkyUeGLXz/KV548tdYhXTVJdMWq2PFzD3XKfcdmGZoNKFdctmUZ57tvYWJhlgvPPkscXn+vbteezQDomsNIbZLZrohCM2FHmKP6/edeU3fkg3d3ynMTPewqnOeSV2fCq3fedy+m+OWIWbtBo0+R1EMsFaKBsMdG2aC0g9Fg0CQWZImFDg1JaqMTC8tYqNgiSVOMHbD7RMrm8eVrHf+Tf7zu6xZi3Vm97YXETUop5SulnlZKvbi0DeBvr3VMQgghrt+j94yiFDz85adQqv36ZiOJrlgVdt6ntHMEAC+ICN0WJjjPHckl7uhtcb57kecX88ycq133uYbfdwBoj2BWl+bRfkaaOkQ5h0MH/+k1dXNDvWC3v+YDz9TxA5gvNKgVIlDtB3CTtCAGZSDtz6FrCVWdYKWGQhYxEnZTCrs45pfJVEKQFHBaFj1z4AU5CqFDOctwg5ggDomUw7bplOG55X100yi+7usWYr0xQGrrFQ+xoYXAQ8aY/cAB4MNKqQfXOCYhhBDXabDk84WH24upfuHhWxks+Wsc0dWTJxCxanb+wgeXSor4FDxTPsPm+Uv0dim022JenWb2Qg1jzJu2s5KRD9zdmac70cozvNCkWNvPV7vu4w9zVTKTdermR/rQS3M
2021-12-07 19:08:54 +01:00
"text/plain": [
2021-12-10 00:06:07 +01:00
"<Figure size 1152x288 with 3 Axes>"
2021-12-07 19:08:54 +01:00
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
2021-12-14 16:17:42 +01:00
"N_deltas = 0*required_N_samples + np.arange(1, 50*required_N_samples+1, 1) # set !={0, required_N_sample} for imperfect FT\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 17:39:00 +01:00
"\n",
"# flags\n",
"unfold_phases = True\n",
"N_sub_max = 2 # how many peaks to get\n",
"enable_single_max = True\n",
2021-12-14 16:17:42 +01:00
"diffs_on_single_max = True\n",
2021-12-10 17:39:00 +01:00
"keep_single_max_phase = True\n",
"\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 00:06:07 +01:00
"\n",
"retrieved_phases = np.empty_like(N_deltas)\n",
2021-12-14 16:17:42 +01:00
"submax_amplitudes = np.empty((len(N_deltas), N_sub_max))\n",
"\n",
2021-12-10 17:39:00 +01:00
"fig, axes = plt.subplots(1,2, figsize=(16,4))\n",
2021-12-10 00:06:07 +01:00
"\n",
"for i, N_delta in enumerate(N_deltas):\n",
2021-12-10 17:39:00 +01:00
" # require N_sub => 2\n",
" N_sub_max = max(2, N_sub_max)\n",
" idx_max = np.empty(N_sub_max, dtype=np.int)\n",
2021-12-14 16:17:42 +01:00
" fft_power_max = np.zeros(N_sub_max)\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 17:39:00 +01:00
" time = np.arange(required_N_samples+N_delta) / sample_rate\n",
" \n",
" f_delta = sample_rate/(required_N_samples+N_delta)\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 17:39:00 +01:00
" fft, ft_freqs = ft_spectrum(signal_func(2*np.pi*f*time), sample_rate)\n",
2021-12-07 19:08:54 +01:00
" fft_power = np.abs(fft)**2\n",
2021-12-10 17:39:00 +01:00
" \n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 17:39:00 +01:00
" idx_single_max = None\n",
2021-12-07 19:08:54 +01:00
" for sub in range(len(idx_max)):\n",
" idx = np.argmax(fft_power)\n",
" idx_max[sub] = idx\n",
2021-12-14 16:17:42 +01:00
" fft_power_max[sub] = fft_power[idx]\n",
2021-12-10 17:39:00 +01:00
" fft_power[idx] = 0 # mask current fft power\n",
" \n",
" if f_delta < np.abs(ft_freqs[idx] - f):\n",
" idx_single_max = idx\n",
" continue\n",
"\n",
2021-12-14 16:17:42 +01:00
" submax_amplitudes[i] = fft_power_max\n",
"\n",
2021-12-10 17:39:00 +01:00
" # No use to interpolate when the max-amplitude frequency\n",
" # is within the frequency resolution of f\n",
" if enable_single_max and idx_single_max is not None:\n",
" freqs = ft_freqs[idx_single_max]\n",
" angles = np.angle(fft[idx_single_max])\n",
" \n",
" if keep_single_max_phase:\n",
" retrieved_phases[i] = angles\n",
" \n",
" l = axes[0].plot(freqs, angles, '1', label=r'$\\Delta N = {}$'.format(N_delta))\n",
" \n",
" axes[1].plot(N_delta/required_N_samples, angles, '1', color=l[0].get_color())\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-14 16:17:42 +01:00
" if not diffs_on_single_max:\n",
2021-12-10 17:39:00 +01:00
" continue\n",
" \n",
" freqs = ft_freqs[idx_max]\n",
" angles = np.angle(fft[idx_max])\n",
2021-12-10 00:06:07 +01:00
"\n",
" # fold angles down for higher submax frequencies\n",
2021-12-10 17:39:00 +01:00
" if unfold_phases:\n",
2021-12-10 00:06:07 +01:00
" folds = 0\n",
2021-12-10 17:39:00 +01:00
" for j in range(len(freqs) - 1):\n",
" if freqs[j] < freqs[j+1] and angles[j] < angles[j+1]:\n",
2021-12-10 00:06:07 +01:00
" folds += 1\n",
2021-12-10 17:39:00 +01:00
" angles[j+1] += - 2*np.pi*folds\n",
2021-12-10 00:06:07 +01:00
"\n",
" if False:\n",
2021-12-10 17:39:00 +01:00
" print(freqs[j], freqs[j+1], \"\\t|\", folds, \"\\t|\", angles[j], angles[j+1])\n",
"\n",
"\n",
2021-12-10 00:06:07 +01:00
"\n",
" # plot frequencies and angles\n",
2021-12-10 17:39:00 +01:00
" axes[0].plot(freqs, angles, '--', alpha=0.5, label=r'$\\Delta N = {}$'.format(N_delta))\n",
" sc = axes[0].scatter(freqs, angles, c=np.arange(len(freqs),0, -1), cmap='Spectral')\n",
2021-12-10 00:06:07 +01:00
" \n",
" # find interpolation between peaks to get the original phase\n",
2021-12-10 17:39:00 +01:00
" dphi_df = (angles[0]-angles[1])/(freqs[0]-freqs[1])\n",
" offset = angles[1] - dphi_df * freqs[1]\n",
2021-12-10 00:06:07 +01:00
" \n",
2021-12-10 17:39:00 +01:00
" angle_at_f = dphi_df * f + offset\n",
2021-12-10 00:06:07 +01:00
"\n",
" # modulo phase\n",
2021-12-10 17:39:00 +01:00
" if not unfold_phases:\n",
" angle_at_f = phase_modulo(angle_at_f)\n",
" \n",
" retrieved_phases[i] = angle_at_f\n",
" axes[0].plot(f, angle_at_f, 'g^')\n",
2021-12-10 00:06:07 +01:00
" \n",
" # Try to fix the midpoints of each line\n",
2021-12-10 17:39:00 +01:00
" if False:\n",
" freq_midpoint = (freqs[0]-freqs[1])/2 + freqs[1]\n",
" angle_midpoint = (angles[0]-angles[1])/2 + angles[1]\n",
" interp_angle_midpoint = dphi_df*freq_midpoint + offset\n",
" \n",
" # modulo phase\n",
" if not unfold_phases:\n",
" angle_midpoint = phase_modulo(angle_midpoint)\n",
" interp_angle_midpoint = phase_modulo(interp_angle_midpoint)\n",
" \n",
" l = axes[0].plot(freq_midpoint, angle_midpoint, '+')\n",
" axes[0].plot(freq_midpoint, interp_angle_midpoint, 'x', color=l[0].get_color())\n",
2021-12-07 19:08:54 +01:00
" \n",
2021-12-10 00:06:07 +01:00
" \n",
"# plot retrieved phases\n",
2021-12-10 17:39:00 +01:00
"axes[1].plot(N_deltas/required_N_samples, phase_modulo(retrieved_phases), '.--')\n",
2021-12-10 00:06:07 +01:00
"\n",
"cbar = fig.colorbar(sc, ax=axes[0])\n",
"cbar.ax.set_ylabel(\"Power ordering\")\n",
2021-12-10 17:39:00 +01:00
"cbar.set_ticks([sc.colorbar.vmin, sc.colorbar.vmax])\n",
"\n",
2021-12-10 00:06:07 +01:00
"\n",
"## horizontal lines\n",
"hlines = [\n",
" (0, None),\n",
" (-np.pi/2, r'$\\frac{-\\pi}{2}$'),\n",
"]\n",
2021-12-07 19:08:54 +01:00
"\n",
2021-12-10 00:06:07 +01:00
"ytrans = axes[0].get_yaxis_transform()\n",
"for location, label in hlines:\n",
" axes[0].axhline(location, alpha=0.1, color='k')\n",
" axes[0].text(-0.06, location, label, transform=ytrans, verticalalignment='center')\n",
" \n",
"axes[0].plot(f, phase_to_retrieve, 'r*')\n",
2021-12-07 19:08:54 +01:00
"axes[0].set_xlabel(\"$f$\")\n",
"axes[0].set_ylabel(r\"$\\varphi_f$\")\n",
"axes[0].axvline(f, alpha=0.1)\n",
2021-12-10 00:06:07 +01:00
"\n",
"\n",
"\n",
2021-12-10 17:39:00 +01:00
"axes[1].set_xlabel(r\"$\\Delta N / N_\\mathrm{required}$\")\n",
2021-12-10 00:06:07 +01:00
"axes[1].set_ylabel(r\"$\\varphi_f$\")\n",
2021-12-10 17:39:00 +01:00
"axes[1].axhline(phase_to_retrieve, alpha=0.1, color='r')\n",
2021-12-10 00:06:07 +01:00
"\n",
2021-12-10 17:39:00 +01:00
"# zooming\n",
"if True:\n",
" x_res = 100\n",
" y_min = -2\n",
" y_max = 0.5\n",
" \n",
2021-12-14 16:17:42 +01:00
" if False:\n",
2021-12-10 17:39:00 +01:00
" x_res = 0.3\n",
" y_min = phase_to_retrieve - 0.1\n",
" y_max = phase_to_retrieve + 0.1\n",
2021-12-10 00:06:07 +01:00
"\n",
2021-12-10 17:39:00 +01:00
" axes[0].set_xlim(f-x_res, f+x_res)\n",
" axes[0].set_ylim(y_min, y_max)\n",
" \n",
2021-12-10 00:06:07 +01:00
"\n",
2021-12-10 17:39:00 +01:00
"plt.show()"
2021-12-07 19:08:54 +01:00
]
2021-12-14 16:17:42 +01:00
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7cAAAEbCAYAAAABJXO/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9a7Q1V1km+rx1WWvv/X0JCQmgEFCGIt129xEV1La1B+IVETnaw4P0YTR4UI6nvfZBWlrRJgreULSPirQoRtKCXHREhqBo24CAQYlggEQDIYTcky/ffd/Wqqo5z4+qWWvWba2q+c69Z631zWeMjOy9vzWrZs1VNWs+83ne9yUpJTw8PDw8PDw8PDw8PDw81hmB6w54eHh4eHh4eHh4eHh4eHDhya2Hh4eHh4eHh4eHh4fH2sOTWw8PDw8PDw8PDw8PD4+1hye3Hh4eHh4eHh4eHh4eHmsPT249PDw8PDw8PDw8PDw81h6e3Hp4eHh4eHh4eHh4eHisPTy59fDw8PC4pEFETyeiexjtX0dEP2WzT8cNIpJE9IXFz1avh4iuI6JX2jqeh4eHh4dHFzy59fDw8PBYexDRnUR0QES7RPRAQahOHsF5XkhEH9D/JqX8finlz9o+lyvo18Ml/h4eHh4eHscJT249PDw8PDYFz5ZSngTwFABfCuC/OO6Ph4eHh4eHxzHCk1sPDw8Pj42ClPIBAO9GTnIBAEQ0JaJfJqK7iOjBwnq73daeiF5GRJ8mootEdCsRfUfx938O4HUA/nWhEJ8r/l7abonoH4no27RjRUT0MBF9WfH7VxHR3xDROSK6mYie3nUdXf0o/u2FRPRBIvrV4lh3ENFXF3+/m4geIqIXaJ+/rrjmvyyO9z4i+ryO815HRK8kohMA/gzAY4vr3SWix9ZtxnV1l4i+lIg+UpznLQC2asf/NiL6h6Lff0NE/5v2bz9ORPcWbW8joq/vGh8PDw8PD486PLn18PDw8NgoENE1AJ4J4Hbtz78I4IuQE94vBPA4AD/dcYhPA/haAI8AcC2A/0FEnyul/EcA3w/gRinlSSnlFS1t3wzgedrv3wzgYSnlR4jocQDeCeCVAB4J4McA/BERPWpIP7R//0oAHwNwFYA3AfhDAE8rru/5AH6jZs3+PwH8LICrAfwDgD/oOC8AQEq5h3wc7yuu96SU8r5lbYhoAuAGANcX1/g2AP9O+/cvA/AGAP930e//DuAdxebDkwH8IICnSSkvQz52dy47n4eHh4eHhw5Pbj08PDw8NgU3ENFFAHcDeAjAfwUAIiIA3wfgP0kpz0gpLwL4OQDf3XYQKeXbpJT3SSmFlPItAD4F4Ct69uFNAL6diHaK3/998TcgJ5zvklK+qzj2XwK4CcC3GvbjM1LK35NSZgDeAuDxAH5GSjmTUv4FgDlyoqvwTinlX0spZwB+ErkC/fie19UXXwUgBvBrUspESvl2AB/W/v37APx3KeXfSikzKeXvA5gV7TIAUwBfTESxlPJOKeWnLffPw8PDw2OD4cmth4eHh8em4H8vFL+nA/hnyBVKAHgUgB0Af19YYc8B+PPi7w0Q0X/QbLPnAPxL7VhLIaW8HcA/Anh2QXC/HQty+3kAvksdtzj21wD43LZj9ejHg9rPB8X563/Tldu7tX7uAjgD4LF9rmsAHgvgXiml1P72We3nzwPwktoYPB7AY4ux+1EArwDwEBH9IRHZ7p+Hh4eHxwbDk1sPDw8Pj42ClPJ9AK4D8MvFnx5GTvT+hZTyiuK/RxTJpyoo4lBfj9wee1VhPf4EAFKH79EFZU1+DoBbC9IG5OTyeq0PV0gpT0gpf8GgHyYoVdrCrvxIAEttxmi/3j3kmwUKn6P9fD+AxxVqucITtJ/vBvCq2hjsSCnfDABSyjdJKb8GOQmWyO3kHh4eHh4eveDJrYeHh4fHJuLXAHwjET1FSimQE8VfJaJHAwARPY6Ivrml3QnkpOpU8bnvQa6YKjwI4JoitrQLfwjgmwD8P1iotgDwP5Arut9MRCERbRXJmK4x6IcJvpWIvqbo+88C+Fsp5d0r2jwI4CoieoT2t38ojvVIIvoc5Gqrwo0AUgA/XCTT+k5UrdSvB/D9RPSVlOMEET2LiC4joicT0TOIaArgEPmGRMa7ZA8PDw+PSwme3Hp4eHh4bByklKcAvBHATxV/+nHkCaY+REQXAPxPAE9uaXcrgF9BTtIeBPCvAHxQ+8j/AnALgAeI6OGOc99ftP9q5LGw6u93I1dzfwI5ab0bwEvR8i7u0Q8TvAl5HPIZAF+OPMHUUkgp/wm5En1HYSN+LPJkUTcjT/b0F6he4xzAdwJ4IYCzAJ4L4I+1f78JedztbxT/fnvxWSCPt/0F5Er7AwAejXysPDw8PDw8eoGqYTEeHh4eHh4emwYiug7APVLKl7vui4eHh4eHx1HBK7ceHh4eHh4eHh4eHh4eaw9Pbj08PDw8PDw8PDw8PDzWHt6W7OHh4eHh4eHh4eHh4bH28Mqth4eHh4eHh4eHh4eHx9rDk1sPDw8PDw8PDw8PDw+PtUfkugNcBEEgt7e3XXfDw8PDw8PDw8PDw8PDwzL29/ellLKXKLv25HZ7ext7e3uuu+Hh4eHh4eHh4eHh4eFhGUR00Pez3pbs4eHh4eHh4eHh4eHhsfbw5NbDw8PDw8PDw8PDw8Nj7eHJrYeHh4eHh4eHh4eHh8faw5NbDw8PDw8PDw8PDw8Pj7XHsZFbInoDET1ERJ/o+Hciov+PiG4noo8R0ZcdV988PDw8PDw8PDw8PDw81hvHqdxeB+Bblvz7MwE8qfjvxQB+6xj65OHh4eHh4eHh4eHh4bEBODZyK6X8awBnlnzkOQDeKHN8CMAVRPS5x9M7Dw8PDw8PDw8PDw+PSwtSSgghXXfDGsYUc/s4AHdrv99T/M3Dw2PkeMENL8AP/9kPu+6Gh4eHh4eHh4dHTwgh8fefPYvbT+267oo1jIncUsvfWrcRiOjFRHQTEd2UpukRd8tj7Hjth1+LO8/d6boba4233vJW3HvhXuP2b7z5jfj1v/t14/YfvOuDeMV7X2Hc3sPDw8PDw8PDYxgeujjDuf0E23HouivWMCZyew+Ax2u/XwPgvrYPSil/W0r5VCnlU6MoOpbOeYwT5w/P4wfe9QP4+jd+veuurC1m6QzPfftz8Yw3PsNZH77m974G177vWuP2p/ZO4UV/8iIcJAcWe+Xh4eHh4eHhsbk4vTdDFBKuuXLbdVesYUzk9h0A/kORNfmrAJyXUt7vulNjx3vvfC+e/8fPh5Sb45UfAiEFAODMwbJw7s3G7Wdux2U/fxk+febTRu3VGH723GdtdutY8f673o83/MMbcOupW43apyLFs9/8bPztPX9r3IdMZOVYenh4eHh4eHiMHef3E1y5MwFRm4F2PXGcpYDeDOBGAE8monuI6EVE9P1E9P3FR94F4A4AtwN4PYD/eFx9W2d84/XfiD/4+B8gEYmT87/3zveCriU8uPugUfuzB2fxPX/yPbg4u2jUXrY71y8pXH/z9did7+L6j11v1D6TGYD1HstM5NegrmUozh6cxZ9+8k9x4z03GvfhmX/wTLz0L15q3P6u83fhQ/d8yLi9h4eHh8dm4aG9h1ibpneeuxOn908bt7/7/N34yP0fMW5/7vAc3njzG43bCynwqr9+Fc4enDU+xmtufA3e/9n3G7d/w0ffgFf+9SuN2//57X+OL3ndlyDJzNbpnzn7GdC1xPoe6FrCT/7VTzb+LoTEQZLhxHSzXLDHmS35eVLKz5VSxlLKa6SUvyulfJ2U8nXFv0sp5Q9IKb9ASvmvpJQ3HVff1hlKsXWl3L7mxtcAgPGi/Bc+8Au47h+uw2s//Fqj9jaUshv+6Qb81R1/xT6OKV530+sw+dmJ8bVwd9vUeceg/pv2QZFaRXKPuz2QLyLuPH+ncftf+uAv4blvf65x+9sevg3f9qZvw2F6aHyMjz34MaTC5zHw8HANKSVrTj69f5pFah7efxgfuOsDxu0P00P8yt/8ivF8IqXET/zVTxi7cQDgJe9+CX7xA79
"text/plain": [
"<Figure size 1152x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Relative Amplitudes of the maximum and submaximum amplitudes plot\n",
"fig2, axs2 = plt.subplots(1,1,figsize=(16,4))\n",
"\n",
"relative_amplitudes = np.empty_like(submax_amplitudes)\n",
"for i, amps in enumerate(submax_amplitudes):\n",
" relative_amplitudes[i] = amps[1:] / amps[0]\n",
" \n",
"\n",
"axs2.plot(1+N_deltas/required_N_samples, relative_amplitudes[:,1], 'g')\n",
"\n",
"axs2_submax = axs2.twinx()\n",
"axs2_submax.plot(1+N_deltas/required_N_samples, submax_amplitudes[:,0], alpha=0.3)\n",
"axs2_submax.plot(1+N_deltas/required_N_samples, submax_amplitudes[:,1], alpha=0.3)\n",
"axs2_submax.set_yticks([],[])\n",
"\n",
"axs2.set_title(\"Relative amplitudes\")\n",
"axs2.set_xlabel(r\"$N_\\mathrm{samples} / N_\\mathrm{required}$\")\n",
"axs2.set_ylabel(r\"$A_1/A_2 \\; (\\delta_1)$\")\n",
"\n",
"\n",
"\n",
"plt.show()"
]
2021-12-07 19:08:54 +01:00
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
