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": 218,
   "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": 3,
   "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": 4,
   "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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAEuCAYAAAC3Tv7YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xc9ZXw/8+5M6NRL65ykwvYxgaDbXpsugmGEGoSICRkEwjJbnZDwpNN3yckG3bDPslvSSDZhMACKQRCCSU4oRnbFGPj3uQqy5blKqu3qd/fH/eOPJJGMyNppNHI581LL8/MbefeuSPm6HyLGGNQSimllFJKKaWGOivdASillFJKKaWUUsnQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppZRSSimlVEbQBFYppYYREXlCRH6chuNuFZFLB2C/lSKyKNX7TRURmSki60WkSUS+msY4BuT6Z5LoeyXV16M/n6uhco8opdRw4U53AEop1V8iUgmMBUJRL88wxhxMT0QnH2PM6emOIU2+CSwzxsxLZxAn8fWPKfp6OL8f7jLGvJmmcIbEPaKUUsOFVmCVUsPFx40x+VE/nZJXEdE/2KmBMBnYmu4g1JDW4z2iv5eUUqr3NIFVSg1bTpPCb4nIJqBFRNwiMl5EnheRYyKyt2uTPhGZJyLrnOZ+z4jI05GmgyJiROTUqHU7NSuMt28nlm+IyCYRaXD2nR21fJKIvOBse1xEHnZe/1cReb5LjA+JyIOx4gWi9zlLRJaJSL3TpPK6GNfnX52YWkTkMREZKyJ/c/b3poiUOOt+W0T2OK9vE5EbY+xrUdTjHs81xnbfcfZZJyKPx1h3bpzr1mNczntf7SzbISJXJHqfYsTX4zUUkaXAZcDDItIsIjN6OL9+X2MROUVEakVkftQ51IjTTDbG9U/qmM76ie7rXu2vy/nHO6fexpnMvRK97iIR+T1QBrzivEffTPKc432uenP/dLtHpJe/l7rGIlG/kxKdS6JYJc5nVfr4O6mH6/AZJ447nffuNRHJ72l9pZSKyxijP/qjP/qT0T9AJbCoh9c3AJOAHOw/2q0F/i+QBUwDKoCrnPWzgH3A1wEP8AkgAPzYWW6AU6P2/0TUskT7rgRWA+OBEUA58GVnmQvYCPw3kIf9ZXmhs2wc0AIUO8/dwFHg7HjxOs93A9911rscaAJmdrk+H2A3v57g7HcdMA/wAkuBHzjrftKJ3QJucWIaF+s9iHeuPbxHW5z3aATwXuSaJrOvnuICZgJVwHhnvSnAKYnepy6xJXMNl2E3T413b6bqGn/ROf9c4DXgp3Guf1LHTHRf92V/Xc6/x3PqQ5zJ3CuL4j1O5pyJ/7lK+v7p6R6hF7+X4sWS6FySiZUePl/08XdSD+d/LVANzAHqnH93AHen8/8b+qM/+pO5P1qBVUoNFy+KXSWrF5EXo17/hTGmyhjTBpwLjDbG/MgY4zfGVAC/BW511r0A+0vig8aYgDHmOeDDJI+faN+RWA4aY2qBV4C5zuvnYX+B/FdjTIsxpt0Y8y6AMeYQsAI7EQBYDNQYY9YmiPcCIB/4iRPPUuCvwG1d4n7IGHPEGFMNvAOsMsasN8b4gL9gJxMYY551Yg8bY54Bdjlx96Snc43lYec9qgXujxFjj/uKE1cIOwmaLSIeY0ylMWYPyb1PEclew0RSco2NMb91XluFnUR8r7/HTPU5dJXEfdPb/Sa6V1Ih3ueqN/dPPMn+Xhro30mRWLp+vvr6OymWH2L3A7aAFmPMZuB97AQep+r+rogsSfK8lFInOe17oZQaLm4wsQdpqYp6PBkYLyL1Ua+5sL84g/2FrdoYY6KW70vy+In2DXA46nGrczywv8jtM8YEe9j3k8A/Yn/5/Azw+yTiHQ9UGWPCXZZN6LLvI1GP22I8zwcQkTuAe7ErmTivj+ohXuj5XGOJfo/2xVi3x331FJcxZreIfA24DzhdRF5z1kvmfYpI9homkspr/FvgZezqla+/x+yFPu0viXPq7X4T3SupEO9z1Zv7J55kfy8N9O8kiP356uvvpE5EZAwwHzsxvh678gwwGrvaDvYfIR50knOllEpIK7BKqeEu+otfFbDXGFMc9VNgjLnGWX4ImCAiErVNWdTjVuzmmxGlvdh3PFVAmfQ8oMuLwJkicgZ2c7w/JhHvQWCSiFhdllUnEU8nIjIZ+4vqPwMjjTHF2E05Je6GyZsU9bgMO/Z+x2WMecoYsxD7i7wBHqB371PKrmF/z8VZJx94EHgMuE9ERqTo8PHu6z4boPumL/eKifFavHOO97nqz+e8p5ji7TPR76R455KO30ldjQKC2E3vzwI2ikgecDGwREQuwm6i/3URuTaJuJRSShNYpdRJZTXQ6AygkiMiLhE5Q0TOdZavxP6y9VVnYJWb6NzccQPwaWe7xcAlvdh3orgOAT8RkTwRyRaRBZGFxph24DngKWC1MWZ/EvGuwu6n9k0R8Yg92M/HgaeTu1Sd5GF/4T4GICKfB87ow3568hURmegkZN8FnulvXGLPvXm5iHiBduxqXojevU+pvIZ9PpcoPwfWGmPuAl4Ffp2iY8e7r/tjIO6bvtwrR7D7f0aLd87xPlf9+Zz3JN4+E/1Oincu6fid1FUl4MOuvp6F3Qf3f4E/GGP2GGPeATYDFxlj/ppEXEoppQmsUurkYYwJYScgc4G9QA3wKFDkLPcDNwH/gD3YyC3AC1G7uMfZvh64HbsKkdS+k4zrVGA/cMA5drQnsQc/+X3Udj3G6yy7DrjaieVXwB3GmO2J4okR3zbgZ9hfpo84cbzX2/3E8RTwOvaX2wrsAWj6G5cX+An2uR8GxgDf7c37lMpr2M9zQUSux+5r+GXnpXuB+SJyewoO3+N93R8DdN/05V75T+D7Tv/4bzivxfssx/tc9flz3pN4+0zid1KP55KO30kx9tMK3AU8BFyB/X7tBb4GEGnd0KWZvlJKxSWdu1UopZSKJiJPAAeMMd9PcxxlwHag1BjTmM5YUklEKrFHaI3Vf1mpDnqv2DLxd5KIjAP2AHnR/XlFZDr2QFF3D2iwSqlhRSuwSik1xDlVinuBp4dT8qqUykx9+J10GrDLdK+anIU9XY9SSiVNRyFWSqkhzBnw5Aj2yKOL0xyOUuok18ffSTOBnV1f1JGHlVJ9oU2IlVJKKaWUUkplBG1CrJRSSimllFIqI2gCq5RSSimllFIqI2gCq5RSSimllFIqI2gCq5RSSimllFIqI2gCq5RSSimllFIqI2gCq5RSJyER2Soilw7CcZ4QkR8P9HFiHHfAzk9EKkVk0UDse6hIdP2G6jUY6PtaREaLyBsiUicij/V2uVJKqf7TeWCVUmqYEZHmqKe5gA8IOc+/ZIz5ozHm9MGPbPAM9/MbaF2vn4hUAncZY95MT0TJGYT3/TvALmPMlQAiMh5YbYyZGGu5Ukqp1NMKrFJKDTPGmPzID7Af+HjUa39Md3xK9YWIDIU/ui8Cno16fg3w9zjLlVJKpZgmsEopNYSIyLdF5Lkur/1cRH7hPP6WiFSLSJOI7BCRK/p4nI4moM7jfxWRTSLSIiKPichYEfmbc5w3RaQkatvxIvK8iBwTkb0i8tWoZfNEZJ2z3TNAdpfjzhKRZSJS7zT3vK5LTL2J49sissdZtk1Ebox1flHPv+Hsu0FEnhGRTrF1Wfc7zj7rROTxGOvOjbWveDE5y2O+f/GuaQ8xfsZZ/04nxtdEJD/BNp8XkVeinu8WkT9HPa8Skbldr5+I/B4oA14RkWYR+Wa8axDjuInu6UTv47dEZBPQ4twfz3fZ10Mi8mCMuOO+5yIyX0TWO8d91lkes7m7iGSJSAMwx7kOm51F1wBL4ixXSimVasYY/dEf/dEf/RkiP8BkoBUodJ67gEPABcBMoAoY7yybApySYH+VwKJ4rzuPPwDGAhOAo8A6YB7gBZYCP3DWtYC1wP8FsoBpQAVwlfN8H/B1wAN8AggAP3a
      "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 maximum amplitude 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Required signal length is: 0.0004s\n",
      "Required number of samples: 4.0\n",
      "Phase to be retrieved: -1.5707963267948966\n"
     ]
    }
   ],
   "source": [
    "sample_rate = 1/1e-4 # Hz\n",
    "f = 2500 # 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": 282,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3000.0 -2.199114857512855\n",
      "3000.0 -0.9424777960769379\n",
      "3000.0 -3.141592653589793\n",
      "2500.0 -1.5707963267948961\n",
      "2500.0 -1.5707963267948966\n",
      "2500.0 -3.141592653589793\n",
      "2142.857142857143 -0.8975979010256547\n",
      "2142.857142857143 0.8975979010256534\n",
      "2142.857142857143 0.0\n",
      "1875.0 0.3755285390725924\n",
      "1875.0 -2.628699773508142\n",
      "1875.0 -2.25317123443555\n",
      "-1.5707963267948966\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAEMCAYAAAD9KlrlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXzddZ3v8dfnLNmbNHv3fV/olkIpCC2LsokoICCg4gg4I4KO3lHHmdGZzp3LDI7jhkARio5YVBAFRMSFpaUt0H2hLZSu6Zak6ZKlyUnO+dw/TlpL6ZK2SX4nyfvp4/foOb/zzfm9Q2LO+ZzvZu6OiIiIiIiISKoLBR1AREREREREpDVUwIqIiIiIiEinoAJWREREREREOgUVsCIiIiIiItIpqIAVERERERGRTkEFrIiIiIiIiHQKKmBFRKTdmFl/M3vJzNaa2Rozu+cYbczMvm9mG8xspZlNDiKriIiIpL5I0AFERKRLawa+7O5LzawHsMTM/ujubx3R5nJgeMtxDvBAy78iIiIi76EeWBERaTfuvtPdl7bcrgHWAn2PavYR4KeetAjoaWa9OziqiIiIdAKdsge2qKjIBw0aFHQM6SKaEw5AJGQBJxE5PUuWLKly9+K2er6+pWd5Y2NNq9ru2b95DdBwxKnZ7j77WG3NbBAwCXj96EsC2464X95ybmfrEksq0GuziIi0peO9v+mUBeygQYNYvHhx0DGki6iuiwFQkJ0WcBKR02NmW9ry+Roba7hyxr+1qu1Pf/vJBncvO1k7M8sBngK+6O4Hjn74GF/irQogKUOvzSIi0paO9/6mUxawIiLSftyMRBuOSDCzKMni9XF3//UxmpQD/Y+43w/Y0WYBREREpMvQHFgREXkvg3g01KrjpE9lZsAjwFp3/85xmj0DfLJlNeJpwH531/BhEREReR/1wIqIyHs4tGUP7HnArcAqM1vecu4fgQEA7v4g8DxwBbABqAdua6uLi4iISNeiAlZERN7L2q6Adff5HHuO65FtHPh8m1xQREREujQVsCIichTDtSq3iIiIpCAVsCIi8h5ukAirgBUREZHUowJWRETepy1XIRYREZGub8mWvSzauIdpQwqZMjC/3a7T7QrYRHOcigVriDfGKD1vHJGsjKAjiYikFDcjHgkHHUNEREQ6iSVb9nLzjxcRa06QFgnx+GentVsR260K2Mo31/Gnq75BvCEGBh5PcP6j/8Dg6y8MOpqISOrQEGIRERFpJXfn2RU7iDUnSDg0NSdYtHGPCtgz1XywkRc/+A/E9te95/y8T99L4eTh5A7tE1AyEZHU0sbb6IiIiEgXVdPQxDeeXs2zK3YQjYSIxxNEIyGmDSlst2t2mwK2/HeL8IQfvh9JS6c51kiiOc47j73AlFmfCTCdiEhq0SrEIiIiciLLtu7l7ieWsWNfA1/+4AjOGVLIG5uqNQe2rcT21eLxBABpWdn0PWsye7dtYd/2rTRWHQg4nYhICmnDfWBFRESk65n96rv81wvrKc3N4Jd3TmPKwAIApg4qaPdrd5sCtvdFk/BEsoCN1ddTV11FwaAhEIH+V54TcDoRkdThZjRHQ0HHEBERkRQVT8CHxvbiPz42nrzMaIdeu9sUsD2G9GHU313N+oeeo7mugYq31xFPNJM1ooR+V6iAFRE5koYQi4iIyJFeWlcBBjNHlnDnBUMwA7OOf7/QbQpYgKn3fY4+l0zh7R//jub6RoZ+4mIG33gRJCC2s5Zon+xAfggiIqnENYRYREREWjQ2x/mvF9bzyPxNTB9ayIwRxYQCfJ/QrQpYM6PfZWfT77Kz33O+ccsBGt7eS+JgE+lDe6qIFZFuTwWsiIiIbKys5Qtzl7FmxwE+PX0QX7t8VOC1UrcqYI8nbUAPEnVNNG46AA7pw1TEikj35WbaB1ZERKSb21xVx1U/mE96JMSPP1nGJWNKg44EqIAFkj2zGaMLIGQ0bm4pYoeriBWR7ktzYEVERLond8fMGFiYxednDuPayf3olZcRdKzDtMxkCzMjY2Q+aQN60LSrDo/Fg44kIhIIN2iOhFp1SOdjZo+aWYWZrT7O4zPMbL+ZLW85/qWjM4qISDCWbd3L5d+bx8bKWsyMz88cllLFK6RAD6yZ9Qd+CvQCEsBsd/9eQFnIGJGPD8ollB7B3Q+fFxHpTlxDiLuyx4AfknztPZ557n5Vx8QREZGgJRLOQ69u5L9fTO7tWtvYHHSk4wq8gAWagS+7+1Iz6wEsMbM/uvtbQYQxM6yleG18Zx8eT5AxqkBFrIh0H2ZaxKkLc/dXzWxQ0DlERCQ1VBxo4O9/uYL5G6q4cnzvQPZ2PRWBF7DuvhPY2XK7xszWAn2BQArY9whBbEstOGSMVhErIt2ICtju7lwzWwHsAL7i7muO1cjM7gDuABgwYEAHxhMRkbbyyPxNLN5Szb0fG88NU/unfM0TeAF7pJZPhCcBrx/jsQ59kTQz0of2BIzGTfvBnYzRhZje1IlIV2cQCnnQKSQ4S4GB7l5rZlcAvwGGH6uhu88GZgOUlZXpl0ZEpJNobI6ze38jAwqz+NKlI/j41P4MLc4JOlarpMwKHGaWAzwFfNHdDxz9uLvPdvcydy8rLi7uqExkDOtJ+pA8YjvqaFhb3SHXFREJkpkTiSZadUjX4+4H3L225fbzQNTMigKOJSIibWRjZS0f+9ECPvno68SaE2REw52meIUU6YE1syjJ4vVxd/910HmOljG0JxYyQlkp8Z9LRKTdqQe2+zKzXsBud3czO5vkh917Ao4lIiJnyN15ckk533xmDWmREPddN4G0TrijQOAVmSUHWT8CrHX37wSd53jSB+cdvt28r5FwbpqGE4tIl2QGobAK2K7KzOYCM4AiMysHvglEAdz9QeA64G/NrBk4CNzoh5blFxGRTqk+1szXnlrFMyt2MG1IAd+9YVLKbY/TWoEXsMB5wK3AKjNb3nLuH1uGLaWcRH0T9Yt3ESnJInNckYpYEemS1APbdbn7TSd5/Ickt9kREZEuIi0coqKmga98cAR/O2MY4U5cwwRewLr7fKDT/BcMZUVJH55Pw9t7wavIHK8iVkS6FsNVwIqIiHRyiYTz2ILNXD2xD0U56Tz+2WmdunA9JPACtjNKH5gLBg3r98LKSjLHF2Phzv/LICICaBViERGRTu7IvV1j8QSfu3BolyheIYVWIe5s0gfkkjmqgKbKgzTtqgs6jpyB3bW7uXru1eyq3RV0FJGUYIZWIRYREemkXlpfweXfm3d4b9c7LxgSdKQ2pQL2DKT170H21F5E+2QHHUXOwH8v/Davb1/ErFdmBR1FJGWEQt6qQ0RERFLHk0vKuW3OmxT3SOe5L5zPjWcPILlmbtehAvYMRXqmY2bE65o4uLoKj6tHojPZWbOTn6+aS8KdOcvnqBdWhOQ+sCpgRUREOo9Di8VfNKqEv5sxlN98/jyGlfQIOFX7UAHbRhI1MZp21VG/vBJvVhHbWcx6dRYJkj+vuMfVCyvSIhz2Vh0iIiISnEN7u976yBs0xRMUZKfxD5eNIiMaDjpau1EB20aivbLJGFtEfG+DithOYmfNTuYsn0NTvAmAWDymXlgRWvaBVQ+siIhISqtpaOKLv1jOV361gqZ4gtqG5qAjdQgVsG0orXc2meOKiO9roH55hYrYFDfr1Vkk/L0/I/XCiiSpgBUREUldy7ft48rvz+e5lTv58qUj+Pnt08jPTgs6VofQNjptLNorGwxi5bVBR5GTWFi+kFg8RoiMw+di8RgLyhcEmEokeMlViFWcioiIpKJEwvnqkyuJJ5xf3DGNskEFQUfqUCpg20G0NJtISRZmluyFdbCoOrtTzbI7lwFQXRcDoKCbfGolclKm3lUREZFUU1HTQI/0KJlpYR64ZTKF2enkZUWDjtXhVFW1EzPD3alfUUndst14UzzoSCIirWJoCLGIiEgqeWl9BZd/dx7/8fxaAIYU53TL4hVUwLYrMyNtQA8SNU3ULa1QESsinYNBKOytOkRERKT9NDbHmfXcW4f3dv3kuQODjhQ4FbDtLFqcReaEYhK
      "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",
    "fig, axes = plt.subplots(1,2, figsize=(16,4))\n",
    "\n",
    "\n",
    "N_sub = 2 # how many peaks to get\n",
    "retrieved_phases = np.empty_like(N_deltas)\n",
    "\n",
    "for i, N_delta in enumerate(N_deltas):\n",
    "    time = np.arange(required_N_samples+N_delta) / 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",
    "\n",
    "    idx_max = np.empty(N_sub, dtype=np.int)\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\n",
    "\n",
    "    # select only the top two frequencies\n",
    "    idx_top = idx_max[:2]\n",
    "\n",
    "    plot_freqs = freqs[idx_top]\n",
    "    angles = np.angle(fft[idx_top])\n",
    "\n",
    "    # fold angles down for higher submax frequencies\n",
    "    if True:\n",
    "        folds = 0\n",
    "        for j in range(len(plot_freqs)-1):\n",
    "            if plot_freqs[j] < plot_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(plot_freqs[j], \"\\t\", plot_freqs[j+1], \"\\t|\", folds, \"\\t|\", angles[j], angles[j+1])\n",
    "\n",
    "    # plot frequencies and angles\n",
    "    axes[0].plot(plot_freqs, angles, '--', alpha=0.5, label=r'$\\Delta N = {}$'.format(N_delta))\n",
    "    sc = axes[0].scatter(plot_freqs, angles, c=np.arange(len(plot_freqs),0, -1), cmap='Spectral')\n",
    "    \n",
    "    # find interpolation between peaks to get the original phase\n",
    "    dphi_df = (angles[0]-angles[1])/(plot_freqs[0]-plot_freqs[1])\n",
    "    offset = angles[0] + dphi_df * plot_freqs[0]\n",
    "    \n",
    "    phase_at_f = dphi_df * f + offset\n",
    "\n",
    "    # modulo phase\n",
    "    phase_at_f = phase_modulo(phase_at_f)\n",
    "    \n",
    "    retrieved_phases[i] = phase_at_f\n",
    "   \n",
    "    axes[0].plot(f, phase_at_f, 'g^')\n",
    "    \n",
    "    # Try to fix the midpoints of each line\n",
    "    freq_midpoint = (plot_freqs[0]-plot_freqs[1])/2 + plot_freqs[1]\n",
    "    angle_midpoint = phase_modulo((angles[0]-angles[1])/2 + angles[1])\n",
    "    axes[0].plot(freq_midpoint, angle_midpoint, '+')\n",
    "    \n",
    "    print(freq_midpoint, angle_midpoint)\n",
    "    print(freq_midpoint, phase_modulo(dphi_df*freq_midpoint + offset))\n",
    "    print(freq_midpoint, phase_modulo(angle_midpoint - dphi_df*freq_midpoint + offset))\n",
    "    \n",
    "# plot retrieved phases\n",
    "axes[1].plot(N_deltas, retrieved_phases, '.--')\n",
    "\n",
    "cbar = fig.colorbar(sc, ax=axes[0])\n",
    "cbar.ax.set_ylabel(\"Power ordering\")\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",
    "print(phase_to_retrieve)\n",
    "axes[0].set_xlabel(\"$f$\")\n",
    "axes[0].set_ylabel(r\"$\\varphi_f$\")\n",
    "axes[0].axvline(f, alpha=0.1)\n",
    "#fig.legend(loc='center left')\n",
    "\n",
    "\n",
    "\n",
    "axes[1].set_xlabel(r\"$\\Delta N$\")\n",
    "axes[1].set_ylabel(r\"$\\varphi_f$\")\n",
    "axes[1].axhline(phase_to_retrieve, alpha=0.1)\n",
    "\n",
    "\n",
    "if False:\n",
    "    res=200\n",
    "    axes[0].set_xlim(f-res, f+res)\n",
    "    axes[0].set_ylim(-2,+0.5)\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "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
}