m-thesis-introduction/02_fourier_phase.ipynb

{
 "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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"
   ]
  },
  {
   "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",
   "execution_count": 4,
   "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",
    "# set signal_func to exp(i*phi)\n",
    "if False:\n",
    "    signal_func = lambda phase: np.exp(1j*phase)\n",
    "\n",
    "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",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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
      "text/plain": [
       "<Figure size 1152x576 with 4 Axes>"
      ]
     },
     "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",
    "plot_submax = True\n",
    "\n",
    "\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",
    "if plot_submax:\n",
    "    phi_f_sub = np.empty(N_delta_f)\n",
    "    f_submax  = np.empty(N_delta_f)\n",
    "\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",
    "    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",
    "    \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",
    "for hline in [0, -np.pi/2]:\n",
    "    ax1.axhline(hline, color='k', alpha=0.5)\n",
    "\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",
    "# 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",
    "if False:\n",
    "    res = 50\n",
    "    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",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\varphi_f$ vs $\\varphi_t$ and the effect of $f/f_\\mathrm{sample}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample_rate = 1/1e-4 # Hz\n",
    "\n",
    "phase_offsets = np.linspace(-np.pi, np.pi, 500, endpoint=True)# rad\n",
    "\n",
    "frequencies = sample_rate * np.array([0.5, 0.49, 0.45, 0.3, 0.28, 0.25, 0.05, 0.001])\n",
    "\n",
    "signal_func = lambda phase: np.sin(phase)\n",
    "\n",
    "\n",
    "# Precreate figure and axis\n",
    "fig, (ax1) = plt.subplots(1,1, figsize=(16,4))\n",
    "\n",
    "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",
    "\n",
    "    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",
    "    \n",
    "    \n",
    "    ax1.plot(phase_offsets, phi_f, '.--', label=\"$f/f_\\\\mathrm{{sample}} = {}$\".format(f/sample_rate))\n",
    "\n",
    "    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",
    "\n",
    "        \n",
    "ax1.set_title(\"Frequencydomain phase of max-power frequency $\\\\varphi_f$ \\n vs Timedomain phase $\\\\varphi_t$ with varying $f / f_\\\\mathrm{{sample}}$\")\n",
    "ax1.set_xlabel('$\\\\varphi_t$')\n",
    "ax1.set_ylabel('$\\\\varphi_f$')\n",
    "ax1.legend(loc='lower right')\n",
    "\n",
    "# grid lines\n",
    "## vertical lines\n",
    "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",
    "## horizontal lines\n",
    "hlines = [\n",
    "    (1, ''),\n",
    "    (-2, ''),\n",
    "    (-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()"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# reconstruct phase from off-frequency ft\n",
    "\n",
    "Require atleast $f_\\mathrm{sample} \\geq 4 f$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Required signal length is: 0.001s\n",
      "Required number of samples: 10.0\n",
      "Phase to be retrieved: -1.5707963267948966\n"
     ]
    }
   ],
   "source": [
    "sample_rate = 1/1e-4 # Hz\n",
    "f = 100 # Hz\n",
    "required_N_samples = sample_rate/f\n",
    "\n",
    "phase_to_retrieve = phase_modulo(-np.pi/2)\n",
    "signal_func = lambda phase: np.cos(phase + phase_to_retrieve)\n",
    "\n",
    "print(\"Required signal length is: {}s\".format(1/f))\n",
    "print(\"Required number of samples: {}\".format(required_N_samples))\n",
    "print(\"Phase to be retrieved: {}\".format(phase_to_retrieve))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8cAAAEPCAYAAABvK7aqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXybd5Xo/895tFqS9zWO4yRO0jRLk6Z1d7pBaWkptHRoKQOlLKUwDAz39WPuBYbh3hk6My8uszCXZZYODJRtmLIUysBAF7pvNGmbdEnTNmkWx46zeF+1nd8fjyRLtuQ4xJbs5Lzz0svSc57lK1mOdM73+3wfUVWMMcYYY4wxxpiTmVPqBhhjjDHGGGOMMaVmybExxhhjjDHGmJOeJcfGGGOMMcYYY056lhwbY4wxxhhjjDnpWXJsjDHGGGOMMeakZ8mxMcYYY4wxxpiTniXHxhhjFhwRWSIiD4jIdhF5UUQ+Weo2GWOMMWZhE7vOsTHGmIVGRBYBi1T1GREpB7YA16rqSyVumjHGGGMWKOs5NsYYs+CoapeqPpO6PwhsBxaXtlXGGGOMWcgsOTbGGLOgicgyYBPwVGlbYowxxpiFzFvqBswndXV1umzZslI3w5zk4kn3VAevIyVuiTnZbdmy5bCq1s/W/hY3btDx8cEZrXukf/eLwFjWottV9fbJ64lIBPgJ8D9UdWBWGmrmFftsNsYYM5um+35T0uRYRK4H/gJYA5ytqpunWdcDbAb2q+rVqWWPAOWpVRqA36nqtSJyCfBz4PVU7Keq+oWjtWfZsmVs3lywCcYURc9wFICasL/ELTEnOxHZM5v7G4sOcdUbb5vRut+766YxVW2fbh0R8eEmxt9X1Z/OQhPNPGSfzcYYY2bTdN9vSt1z/AJwHfCvM1j3k7jnlFWkF6jqhen7IvIT3IQ47ZF0Em2MMWZ+0FkaESEiAnwT2K6q/zArOzXGGGPMSa2k5xyr6nZV3XG09USkBXgr8I0C8XLgjcDPZreFxhhjZo1A0iMzus3ABcBNwBtF5LnU7aq5fQLGGGOMOZGVuud4pv4R+F9MDKGe7B3A/ZPONztPRLYCncCfquqLc9xGY4wx01AgOUs9x6r6KGAn5htjjDFm1sx5ciwi9wFNeUKfU9Wf51k+efurgYOquiV1LnE+7ya3V/kZYKmqDqV6En4GrCqw/1uBWwFaW1uP1hxjjDG/L5m95NgYY4wxZrbNeXKsqpcd5y4uAN6eSnKDQIWIfE9V3wsgIrXA2bi9x+ljDmTd/5WI/JOI1Knq4Tztux24HaC9vV2Ps63GGGMKUBESPruC4IlMRN4C/D/AA3xDVb84KR4AvgOcCRwB3qWqu4vdTmOMMSafef8tRVU/q6otqroMuBH4bToxTrke+C9VzVzyQ0SaUpO1ICJn4z7PI0VstjHGmDySjszoZhae1FUlvg5cCawF3i0iayet9iGgV1VXAl8G/m9xW2mMMcYUVtLkWETeISIdwHnAL0XkN6nlzSLyqxnu5kbgPyYteyfwQuqc468AN6qq9QobY0wJqUDScWZ0MwvS2cBrqrpLVaPAD4FrJq1zDXBH6v6PgTeli9nGGGNMIVv29PL1B15jy57eOT1OSSfkUtW7gLvyLO8Epsw6qqoPAg9OWnZJnvW+BnxtlpppjDFmVsisXcrJzEuLgX1ZjzuAcwqto6pxEekHaoEppz0ZY4wxAE/uOsz7/v1p4okkfq/D9285lzOXVs/Jsaw8b4wxpjhm91JOZv7J94ubPGprJusgIreKyGYR2Xzo0KFZaZwxxpiF6Rdbu4jGkyQVYvEkT+6au7NlLTk2xhhTFAokvM6MbmZB6gCWZD1uwb2cYt51RMQLVAI9k3ekqreraruqttfX189Rc40xxiwE/tT3Ao+Az+twblvtnB1roVzn2BhjzEJnl3I60T0NrBKR5cB+3DlB/nDSOncDNwNP4M4P8lubE8QYY8x0eoejVId83HJhG+e21c7ZkGqw5NgYY0yRKDZk+kSWOof448BvcC/l9O+q+qKIfAHYrKp3A98Evisir+H2GN9YuhYbY4xZCLZ19NO+rIY/vnTlnB/LkmNjjDHFYT3HJzxV/RXwq0nL/nfW/THcSzAaY4wxRzU0Huf1I8Ncd8biohzPkmNjjDFFY7NVG2OMMWamIgEvz33+cnTq3I1zwpJjY4wxRaHWc2yMMcaYY1QZ8hXtWJYcG2OMKQoVIe6zmaiNMcYYMzNff+A1Ksp83HTu0qIcz76lGGOMKRp1ZEY3Y4wxxpjvPbmHp1+fcsW/OWM9x8YYY4rGhlUbY4wxZiYODo7R1T/GhpbKoh3TkmNjjDFFYeccG2OMMWamtu3rB2DjkqqiHdOSY2OMMUWjdp1jY4wxxszAto4+HIF1zRVFO6Ylx8YYY4pDxHqOjTHGGDMjY/EkG5dUEfIXL2W15NgYY0xxCDj2qWOMMcaYGfizq9agWpzrG6fZ1xRjjDFFISiOp7gfcsYYY4xZuESKO+LMLuVkjDGmOAQcR2d0M8YYY8zJ69cvHODtX3uUrv7Roh7Xeo6NMcYUjSW+xhhjjDmaLXt6ePnAILXhQFGPa8mxMcaYohDBhlUbY4wx5qi2dvSzdlEFfm9xBzqXfFi1iFwvIi+KSFJE2qdZb7eIPC8iz4nI5qzlNSJyr4i8mvpZnVouIvIVEXlNRLaJyBnFeD7GGGPyExSvNzmjmzHGGGNOTomk8sL+fja2VBb92CVPjoEXgOuAh2ew7qWqerqqZifRnwHuV9VVwP2pxwBXAqtSt1uBf569JhtjjDlmds6xMcYYY45i56EhRqIJNrRUFf3YJU+OVXW7qu44jl1cA9yRun8HcG3W8u+o60mgSkQWHcdxjDHGHCfHozO6GWOMMebkpApXrm/ijKXVRT92yZPjY6DAPSKyRURuzVreqKpdAKmfDanli4F9Wet1pJYZY4wpAbGe4xNWoVOcJq1zuog8kTqVapuIvKsUbTXGGDO/rW4q55/feybL68JFP3ZRkmMRuU9EXshzu+YYdnOBqp6BO1z6j0XkoqMdNs+yKd+4RORWEdksIpsPHTp0DM0xxhhzrCw5PmEVOsUp2wjwPlVdB7wF+EcRKf6YOWOMMfNa/2isZMcuSnKsqpep6vo8t58fwz46Uz8PAncBZ6dC3enh0qmfB1PLO4AlWbtoATrz7Pd2VW1X1fb6+vpjf3LGGGNmRGRmibElxwtSoVOcMlT1FVV9NXW/E/fz2j54jTHGZETjSc76q/v4yv2vluT4C2JYtYiERaQ8fR+4HHciL4C7gZtT928Gfp61/H2pWavPBfrTw6+NMcYUnwj4fMkZ3cyCU+gUp7xE5GzAD+wsQtuMMcYsEC8fGCCaSLKiPlKS45f8Osci8g7gq7jV41+KyHOqeoWINAPfUNWrgEbgLhEBt80/UNVfp3bxReBOEfkQsBe4PrX8V8BVwGu4Q7k+UKznZIwxJj/rFV64ROQ+oClP6HPHuJ9FwHeBm1U1byUkNbfIrQCtra3H2FJjjDEL1daOfgA2lOAyTjAPkmNVvQt3mPTk5Z24yS2qugvYWGD7I8Cb8ixX4I9ntbHGGGN+byLYTNQLmKpeVigmIt0iskhVuyad4jR5vQrgl8Cfp64kUehYtwO3A7S3t9ubxhhjThLb9vVRE/bTUl1WkuMviGHVxhhjTgyOM7ObWXAKneKUISJ+3GL4d1T1R0VsmzHGmAVia0cfG1sqSY0YLrqS9xwbY4w5SYhNtnUCy3uKk4i0Ax9V1VuAG4CLgFoReX9qu/er6nMlaK8xxph56CMXraA67CvZ8S05NsYYUxSCnXN8oprmFKfNwC2p+98DvlfkphljjFlA/uDMlpIe35JjY4wxRSECXpuJ2hhjjDF5vNI9iCqsbiovWRvszC5jjDFFY9c5NsYYY0w+X7n/VT747adL2gbrOTbGGFMcds6xMcYYYwrY1tGf9xJOX93TzaFonOubqvnRgV7q/V4+sbRxTtpgPcfGGGOKQnAv5TSTmzHGGGNOHr3DUfb2jLChpWpK7F1NNQjw5s2vIKn
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_deltas = np.arange(1, required_N_samples+1, 1) # set !={0, required_N_sample} for imperfect FT\n",
    "\n",
    "\n",
    "# flags\n",
    "unfold_phases = True\n",
    "N_sub_max = 2 # how many peaks to get\n",
    "enable_single_max = True\n",
    "continue_on_single_max = True\n",
    "keep_single_max_phase = True\n",
    "\n",
    "\n",
    "\n",
    "retrieved_phases = np.empty_like(N_deltas)\n",
    "fig, axes = plt.subplots(1,2, figsize=(16,4))\n",
    "\n",
    "for i, N_delta in enumerate(N_deltas):\n",
    "    # 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",
    "\n",
    "    time = np.arange(required_N_samples+N_delta) / sample_rate\n",
    "    \n",
    "    f_delta = sample_rate/(required_N_samples+N_delta)\n",
    "\n",
    "    fft, ft_freqs = ft_spectrum(signal_func(2*np.pi*f*time), sample_rate)\n",
    "    fft_power = np.abs(fft)**2\n",
    "    \n",
    "\n",
    "    idx_single_max = None\n",
    "    for sub in range(len(idx_max)):\n",
    "        idx = np.argmax(fft_power)\n",
    "        idx_max[sub] = idx\n",
    "        tmp = fft_power[idx]\n",
    "        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",
    "    # 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",
    "\n",
    "        if continue_on_single_max:\n",
    "            continue\n",
    "    \n",
    "    freqs = ft_freqs[idx_max]\n",
    "    angles = np.angle(fft[idx_max])\n",
    "\n",
    "    # fold angles down for higher submax frequencies\n",
    "    if unfold_phases:\n",
    "        folds = 0\n",
    "        for j in range(len(freqs) - 1):\n",
    "            if freqs[j] < freqs[j+1] and angles[j] < angles[j+1]:\n",
    "                folds += 1\n",
    "                angles[j+1] += - 2*np.pi*folds\n",
    "\n",
    "                if False:\n",
    "                    print(freqs[j], freqs[j+1], \"\\t|\", folds, \"\\t|\", angles[j], angles[j+1])\n",
    "\n",
    "\n",
    "\n",
    "    # plot frequencies and angles\n",
    "    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",
    "    \n",
    "    # find interpolation between peaks to get the original phase\n",
    "    dphi_df = (angles[0]-angles[1])/(freqs[0]-freqs[1])\n",
    "    offset = angles[1] - dphi_df * freqs[1]\n",
    "    \n",
    "    angle_at_f = dphi_df * f + offset\n",
    "\n",
    "    # modulo phase\n",
    "    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",
    "    \n",
    "    # Try to fix the midpoints of each line\n",
    "    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",
    "    \n",
    "    \n",
    "# plot retrieved phases\n",
    "axes[1].plot(N_deltas/required_N_samples, phase_modulo(retrieved_phases), '.--')\n",
    "\n",
    "cbar = fig.colorbar(sc, ax=axes[0])\n",
    "cbar.ax.set_ylabel(\"Power ordering\")\n",
    "cbar.set_ticks([sc.colorbar.vmin, sc.colorbar.vmax])\n",
    "\n",
    "\n",
    "## horizontal lines\n",
    "hlines = [\n",
    "    (0, None),\n",
    "    (-np.pi/2, r'$\\frac{-\\pi}{2}$'),\n",
    "]\n",
    "\n",
    "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",
    "axes[0].set_xlabel(\"$f$\")\n",
    "axes[0].set_ylabel(r\"$\\varphi_f$\")\n",
    "axes[0].axvline(f, alpha=0.1)\n",
    "\n",
    "\n",
    "\n",
    "axes[1].set_xlabel(r\"$\\Delta N / N_\\mathrm{required}$\")\n",
    "axes[1].set_ylabel(r\"$\\varphi_f$\")\n",
    "axes[1].axhline(phase_to_retrieve, alpha=0.1, color='r')\n",
    "\n",
    "# zooming\n",
    "if True:\n",
    "    x_res = 100\n",
    "    y_min = -2\n",
    "    y_max = 0.5\n",
    "    \n",
    "    if True:\n",
    "        x_res = 0.3\n",
    "        y_min = phase_to_retrieve - 0.1\n",
    "        y_max = phase_to_retrieve + 0.1\n",
    "\n",
    "    axes[0].set_xlim(f-x_res, f+x_res)\n",
    "    axes[0].set_ylim(y_min, y_max)\n",
    "    \n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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
}