m-thesis-introduction/fourier/01_fourier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": "markdown",
   "metadata": {},
   "source": [
    "The Fourier transform (FT) of a quantity $f(t)$ is defined as\n",
    "\n",
    "$$\n",
    "\\hat{f}(\\omega) = F(f(t)) = \\frac{1}{\\sqrt{2\\pi}} \\int_{-\\infty}^{\\infty} \\mathrm{d}t\\, f(t)\\, \\exp{ \\left(-2 \\pi i f t \\right)}\n",
    ".\n",
    "$$\n",
    "\n",
    "In signal processing, taking the integral from $t = -\\infty$ to $t = \\infty$ is not practical. Together with a finite sample space $t$, the integral can be turned into a (discrete) sum, giving a Discrete Fourier Transform (DFT)\n",
    "\n",
    "$$\n",
    "F(t) = \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N-1} f(t_k)\\, \\exp{ \\left(-i \\omega t_k \\right)}\n",
    "\\label{eq:DFT} \\tag{1}\n",
    ",\n",
    "$$\n",
    "where the sum iterates over each sample $t_k$ up to the end ($k=N-1$) of the signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interdependence $\\Delta t$, $T$, $\\Delta f$, $f_{min}$, $f_{max}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The highest frequency $f_\\mathrm{max}$ (also known as the Nyquist frequency) that can be (reliably) measured within a signal is dependent only on the sampling rate ($1/\\Delta t$) of the signal with\n",
    "$$\n",
    "    f_\\mathrm{max} = \\frac{1}{2} f_\\mathrm{sampling} = \\frac{1}{2} \\frac{1}{\\Delta t}\n",
    "    .\n",
    "$$\n",
    "Higher frequencies will fall in between two adjacent samples. This means that these frequencies are then 'measured\\` with differing phases in such a way that they are 'folded' around the Nyquist frequency and appear as aliased frequencies below it.\n",
    "\n",
    "Likewise, the lowest frequency depends on the length of the signal $T$, such that over the whole signal one oscillation can be found\n",
    "$$\n",
    "    f_\\mathrm{min} \\sim \\frac{1}{T}\n",
    "    .\n",
    "$$\n",
    "Even lower frequencies will no longer be detected as a frequency. Instead, they will introduce a bias for the $0 Hz$ component in the spectrum.\n",
    "\n",
    "\n",
    "$$\n",
    "\\delta f = \\frac{1}{T}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FTs of ideal signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "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\n",
    "    \n",
    "\n",
    "def plot_signal( ax, signal, sample_rate = 1):\n",
    "    ax.set_title(\"Signal\")\n",
    "    ax.set_xlabel(\"t (s)\")\n",
    "    ax.set_ylabel(\"A(t)\")\n",
    "\n",
    "    ax.plot(np.arange(len(signal)) / sample_rate, signal)\n",
    "    \n",
    "    return ax\n",
    "\n",
    "\n",
    "def plot_signal_and_spectrum( signal, sample_rate, title=None,\n",
    "                             ft_kwargs = {}, spectrum_kwargs = {}):\n",
    "    \"\"\"Plot a signal and its spectrum in one go.\"\"\"\n",
    "    fig = plt.figure(figsize=(16, 4))\n",
    "    \n",
    "    if title:\n",
    "        fig.suptitle(title)\n",
    "\n",
    "    # setup plot layout\n",
    "    gs0 = gridspec.GridSpec(1, 2, figure=fig)\n",
    "    gs00 = gs0[0].subgridspec(1, 1)\n",
    "    gs01 = gs0[1].subgridspec(2, 1, height_ratios=[3,1], hspace=0)\n",
    "      \n",
    "    # plot the signal\n",
    "    ax1 = fig.add_subplot(gs00[0, 0])\n",
    "    plot_signal(ax1, signal)\n",
    "    \n",
    "    # plot spectrum\n",
    "    signal_fft, freqs = ft_spectrum(signal, sample_rate, **ft_kwargs)\n",
    "    _, (ax2, ax3) = plot_combined_spectrum(signal_fft, freqs, spectrum_kwargs = spectrum_kwargs, fig=fig, gs=gs01)\n",
    "\n",
    "    # return the axes\n",
    "    axes = np.array([ax1, ax2, ax3])\n",
    "    \n",
    "    return (fig, axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_rate = 1/1e-3 # s\n",
    "n_samples = int(1e3)\n",
    "t = np.arange(n_samples) / sample_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FTs: constant signal\n",
    "\n",
    "A constant signal (also termed offset or bias) results in only a single component in the Fourier transform, at the zero frequency.\n",
    "\n",
    "Since often such a component is included in a signal (for example when a signal is not on average zero), it can help to mask out the zero frequency when plotting the spectrum.\n",
    "Doing this for a constant signal, we get a totally flat spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7wAAAEjCAYAAAAR7A2IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5xddX3v/9c7FwiXBEJmRCBcAkxEFC+YcqkMam0VqErVakFbhaOlnmrtafX0p21/6qGnP3usrdbKUVGo2nrwQtVixaJHUYIVJVG5XxICkXDLhQABDOTy+f2xV3QYZpKZzOzsmT2v5+OxH9lrre9e+72XG9d89ve7vitVhSRJkiRJ3WZapwNIkiRJktQOFrySJEmSpK5kwStJkiRJ6koWvJIkSZKkrmTBK0mSJEnqSha8kiRJkqSuZMErSVIXSPL6JN/cBe/zwiSr2v0+kiSNBwteSVLXSPK6JEuSPJzkniTfSHJSG99vXIu/JGcluXJnXltVn6uql4xXFkmSuoEFrySpKyT5U+DDwP8H7A8cAvxv4PRO5pIkSZ1jwStJmvSS7AOcC7y1qr5cVY9U1aaq+lpV/femze5JPpzk7ubx4SS7N9temGRVknckWd30Dp89YP+nJbkxyYYkdyV5Z5K9gG8ABzY9yg8nOTDJcUl+kOSBZj8fTbLbgH1VkrckWZZkfZLz0vJ04OPAic2+Hhjms56VZEWT5fYkrx+w/soB7V6S5JYkDyb530m+l+TNA9sm+WCT4fYkpw547dlJbmreY0WSPxi//7UkSdp1LHglSd3gRGAW8JXttPkL4ATgOcCzgeOAvxyw/anAPsBBwJuA85LMbbZdAPxBVc0Gngl8p6oeAU4F7q6qvZvH3cAW4E+AnibXi4E/HJTlZcCvNDleC7y0qm4C3gL8oNnXvoM/QFNkfwQ4tcnyq8BPh2jXA1wMvBuYB9zStB3o+GZ9D/AB4IIkabatbjLOAc4GPpTk2MHvI0nSRGfBK0nqBvOAtVW1eTttXg+cW1Wrq2oN8D+A3xuwfVOzfVNVXQo8DDxtwLajk8ypqvVV9ePh3qSqllbVVVW1uaruAD4BvGBQs7+pqgeq6mfA5bSK8JHaCjwzyR5VdU9V3TBEm9OAG5re7s20iuR7B7VZWVWfrKotwGeAA2gNBaeqvl5Vt1XL94BvAv2jyChJ0oRgwStJ6gbrgJ4kM7bT5kBg5YDllc26X+xjUMH8KLB38/zVtIrIlc3Q4BOHe5MkC5P8e5J7kzxE65rinkHNBhafA99nu5pe5d+h1RN8T5KvJzlqiKYHAncOeF0BgyfXunfA9kebp3s3n+HUJFclub8ZWn3aEJ9BkqQJz4JXktQNfgBsBH5rO23uBg4dsHxIs26HqurqqjodeArwVeCL2zYN0fxjwM1AX1XNAf4cyBDthnyrEWS5rKp+g1aP7M3AJ4dodg8wf9tCM1R5/hDtnqS5rvlfgQ8C+zdDqy9l5J9BkqQJw4JXkjTpVdWDwHtoXXf7W0n2TDKz6an8QNPsIuAvk/Q217i+B/iXHe07yW7NPW73qapNwEO0rtMFuA+Y10yatc3sps3DTe/rfx3FR7kPmD9wkqtBWfZP8ormWt7HaA273jJE068DxzTHYgbwVlrXKI/EbsDuwBpgczOZlbc7kiRNSha8kqSuUFV/D/wprYmo1tAa0vs2Wj2yAP8TWAJcC1wH/LhZNxK/B9zRDFF+C/C7zXveTKuQXtHMynwg8E7gdcAGWr2vXxjFx/gOcANwb5K1Q2yfBryDVs/0/bSuDR48IRZVtRZ4Da3JqNYBR9P67I/tKEBVbQDeTqsXe33zWS4ZxWeQJGnCSOuyHkmS1K2STKN1De/rq+ryTueRJGlXsYdXkqQulOSlSfZtrsnddh3xVR2OJUnSLmXBK0lSdzoRuA1YC7wc+K2q+nlnI0mStGs5pFnqEkleD7yxqto6uUySFwL/UlUjmvFVkiRJ6hR7eKVJJslJSf4zyYPNPTK/n+RXqupz7S52JUmaKoY737bx/e5I8uvt2r80Vc3odABJI5dkDvDvtG5z8kVatw/pZwQzr0qSpJGZiOfbJDOqanOn3l+arOzhlSaXhQBVdVFVbamqn1fVN6vq2iRnJblyW8MkL0lyS/PL9P9O8r0kb262nZXkyiQfTLI+ye3NvTa3vfbsJDcl2ZBkRZI/2PUfVZKkjtnR+fb7Sf6xOcfenOTF216YZJ8kFyS5J8ldSf5nkukDtv/+gHPsjUmOTfLPwCHA15I8nOTPkhyWpJK8KcnPgO8keWGSVQODDuwZTvK+JF9K8i/N/q9LsjDJu5OsTnJnEkeDaUqx4JUml1uBLUk+k+TUJHOHapSkB7gYeDcwD7gF+NVBzY5v1vfQulfnBUnSbFsNvAyYA5wNfCjJseP9YSRJmqB2dL49HlhB6xz6XuDLSfZrtn0G2AwcCTwXeAmw7Qfn1wDvA95A6xz7CmBdVf0e8DPg5VW1d1V9YMB7vQB4OvDSEWZ/OfDPwFzgJ8BltP7mPwg4F/jECPcjdQULXmkSqaqHgJOAAj4JrElySZL9BzU9Dbihqr7cDH/6CHDvoDYrq+qTVbWF1sn5AGD/5n2+XlW3Vcv3gG/SGsolSVLXG8H5djXw4araVFVfoPUD8m82208F/ltVPVJVq4EPAWc0r3sz8IGquro5xy6vqpU7iPO+Zl8jnWV9cVVd1pz/vwT0An9TVZuAzwOHJdl3hPuSJj0LXmmSqaqbquqsZpbkZwIHAh8e1OxA4M4Brylg1aA29w7Y/mjzdG+A5tfsq5pJOh6gVUD3jO8nkSRp4trB+faueuKtTlY22w8FZgL3JHmgOYd+AnhK0+5gWrcLG407d9zkCe4b8PznwNrmx+1ty9Cc76WpwIJXmsSq6mbg07ROxAPdA/zitkHNUOUR3UYoye7AvwIfBPavqn2BS4Fs94WSJHWpIc63Bw24DAha19/eTas4fQzoqap9m8ecqnpG0+5O4Ijh3mYE6x8B9ty20Fwb3DuazyJNNRa80iSS5Kgk70gyv1k+GDgTuGpQ068DxyT5rSQzgLcCTx3h2+wG7A6sATY3k1k5wYUkacoYwfn2KcDbk8xsrst9OnBpVd1D6zKgv0syJ8m0JEckeUHzuk8B70zyvLQcmeTQZtt9wOE7iHYrMCvJbyaZCfwlrXO2pGFY8EqTywZaE2X8MMkjtE681wPvGNioqtYCr6E1GdU64GhgCSO4nUJVbQDeTus2DOuB1wGXjN9HkCRpwtvR+faHQB+wFvhr4Leral2z7Q20fjy+kdZ59GJa82RQVV9q2v+f5j2+Cmyb7Or9wF82Q6HfOVSoqnoQ+ENahfNdtHp8B1+yJGmAPPHyA0ndKMk0WifE11fV5Z3OI0nSZJXkLODNVXVSp7NI2jF7eKUuleSlSfZtrsn9c1rX4A4e+ixJkiR1LQteqXudSGsmyLW07sn3W6O4pYEkjVmSczqdoRt4HMfOYzh2HsOx8xiO3c4cQ4c0S5KktkiypKoWdTrHZOdxHDuP4dh5DMfOYzh2O3MM7eGVJEmSJHWlrunh7enpqcMOO6zTMSRJXWLp0qVrq8r7W47B3nvvXUcddVSnY0x6a9asobfXr+JYeAzHzmM4dh7DsVu6dOnDVTV7NK+Z0a4wu9phhx3GkiVLOh1DktQlkqzsdIbJ7qijjvLcLEkaN0luGe1rHNIsSZIkSepKFrySJEmSpK5kwStJkkhyYZLVSa4fZnuSfCTJ8iTXJjl2V2eUJGm0uuYaXkmSNCafBj4KfHaY7acCfc3jeOBjzb+SpA7atGkTq1atYuPGjZ2OMm5mzZrF/PnzmTlz5pj3ZcErSZKoqiuSHLadJqcDn63W7R2uSrJvkgOq6p7hXnDfQxtZunI9zzt07jinlSRts2rVKmbPns1hhx1Gkk7HGbOqYt26daxatYoFCxaMeX8OaZYkSSNxEHDngOVVzbphrd7wGK/75FUsXbm+rcEkaSrbuHEj8+bN64piFyAJ8+bNG7ceawteSZI0EkP9JVVPapSck2RJkiUAj2/eylU
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7wAAAEjCAYAAAAR7A2IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5xdVX338c8XEohKEAwRgQABxSLiBZ0iVFS8FIGi1HuQVqEq+qiltdoWWx+xtD62VotaqYqCVwSRoqJG0UcBxUeUBEXuNgYj4SIhhJsKEvJ7/th79DDMNTOTc+bM5/16zWvO3mvtvX9nz4GV31lrr5WqQpIkSZKkfrNZtwOQJEmSJGk6mPBKkiRJkvqSCa8kSZIkqS+Z8EqSJEmS+pIJryRJkiSpL5nwSpIkSZL6kgmvJElDJDkyyTc2wXUOTLJ6jDqnJ/nT6Y5lyDXPT/LqKTjP4iSVZM5UxLWRMWyZ5OokD+9WDJKk7jHhlSQNK8nLkyxLcleSG5N8LckB03i9MZO/CZ7vqCQXbsyxVXVaVR00VbFsrCSPB54AfKnbscxUVXUPcCrw992ORZK06ZnwSpIeIMnfAO8D/g+wPbAL8F/A4d2MaxZ6LXBaVVW3A+lV4+w9/izwyiRbTnc8kqTeYsIrSbqfJA8FTgDeUFVnV9WvqureqvpyVf1tW2fLJO9LckP7877BZGKwpzbJm5Pc3PYOH91x/kOTXJnkziTXJ3lLkocAXwN2bHuU70qyY5J9k3w/yW3teT6YZIuOc1WS1yX5nyTrkpyUxmOADwP7t+e6bYT3elSSlW0s1yY5smP/hR31DkpyTZLbk/xXkgsGh/wO1k3ynjaGa5Mc0nHs0Umuaq+xMslrJ/DnOAS4YEi830tyYntPVib5o3b/de39fmVH/T9J8qMkd7Tl7+gom5fkM0nWtue6OMn2w9yjHZL8JMlb2u2HJjml/Xtcn+Rfkmzelm3e3odbkqwE/mSkN5bkZR1/67uS3JPk/LZsy/Y8v0jyyyQfTvKgtmzw8/X3SW4CPt7uf02SFUluTXJOkh0Hr1VVq4F1wH4TuPeSpD5gwitJGmp/YB7whVHq/CNN8vBEmiG3+wJv6yh/BPBQYCfgVcBJSbZty04BXltV84G9gW9X1a9okrsbqmqr9ucG4D7gTcB2bVzPBl4/JJbDgD9s43gp8Nyqugp4HfD99lzbDH0DbZL9AeCQNpY/An48TL3tgLOAtwILgGvaup2e0u7fDng3cEqStGU3tzFuDRwNnJjkSUOvM0J8u7XnHXqtn7SxfBY4o33/jwL+DPhgkq3aur8CXgFsQ5N8/q/8/nngV9L8jXZuz/U64DdDYlhMk3B/sKre0+7+JLC+vd4+wEHA4PO+r2nf6z7AAPDikd5fVX1u8G8N7AisBE5vi/8NeDTN5+tRNJ+jt3cc/gjgYcCuwDFJngW8i+bvvwOwqr0vna6i+YxIkmYRE15J0lALgFuqav0odY4ETqiqm6tqDfBPwJ93lN/blt9bVUuBu4A/6CjbK8nWVbWuqi4Z6SJVtbyqLqqq9VX1c+AjwDOGVPvXqrqtqn4BnEeTJI3XBmDvJA+qqhur6oph6hwKXNH2dq+nSZJvGlJnVVV9tKruo0kId6AZCk5VfbWqflaNC4BvAE8bR2yDSfqdQ/ZfW1Ufb6/1OZqE9YSquqeqvgH8liZJpKrOr6rLqmpDVf2EJqEcvH/30vytH1VV97X3+o6O6+wFnA8cX1UnA7Q9wIcAf932/N8MnAgsaY95KfC+qrquqm6lSUJHlWQzmsT9/Kr6SPtFwWuAN1XVrVV1J83Q+iUdh21o47qnqn5D83k8taouaZ/ZfStN7/7ijmPu7LinkqRZwoRXkjTUWmC7jP5s5I40vWiDVrX7fneOIQnzr4HBXscX0SSRq9qhwfuPdJEkj07ylSQ3JbmDJvHZbki1zuSz8zqjanuVX0bTs3ljkq8m2XOYqjsC13UcV8DQybVu6ij/dftyq/Y9HJLkonao7W00733oexjO4DDs+UP2/7Lj9W/aaw7dN3jtpyQ5L8maJLfTvNfBa38aOBc4I82w9HcnmdtxniOB62l6twftCsyluV+3te/nI8DgDMj3u1fc/zMykne27/HYdnsh8GBgecc1vt7uH7Smqu7u2L7f57Gq7qL5HO/UUWc+v7+nkqRZwoRXkjTU94G7gdGWwrmBJvkZtEu7b0xVdXFVHU6TJH0ROHOwaJjqHwKuBvaoqq2BfwAyTL1hLzWOWM6tqj+m6ZG9GvjoMNVuBBYNbrQ9kIuGqfcAaZ5r/m/gPcD27dDqpYzjPbQJ+c9ohvZurM8C5wA7V9VDaZ5rTnv+e6vqn6pqL5oh2ofRDH8e9A7gFuCzg8/o0iSz9wDbVdU27c/WVfXYtvxGmh7nQbuMFlySJcARwIur6t529y00SftjO67x0Hbo86Chf9v7fR7b4eALaBL2QY8BLh0tHklS/zHhlSTdT1XdTvO85ElJ/jTJg5PMbXsq391WOx14W5KF7TOubwc+M9a5k2yRZo3bh7YJzh00z+lC03O5IM2kWYPmt3Xuantf/9cE3sovgUXpmORqSCzbJ3l+mxzdQzPs+r5hqn4VeFx7L+YAb6B5hnQ8tgC2BNYA69NMZjWR5Y6W8sAh3BMxH7i1qu5Osi/w8sGCJM9M8rg2mb2DZohz5/u/F3gJ8BDg00k2q6obaYZkvzfJ1kk2S/LIJIMxngkcm2RR+8z2cSMFlmQf4D+BP22HxQNQVRtovng4Me3auUl2SvLcUd7nZ4Gjkzyx/ZLh/wA/aIfBk2Qnmmd+LxrjfkmS+owJryTpAarqP4C/oZmIag1Nz94baXpkAf4FWEYzedJlwCXtvvH4c+Dn7RDl19FMtERVXU2TSK9sh7LuCLyFJkm7kyYJ+twE3sa3gSuAm5LcMkz5ZsCbaXoHb6VJLIdOiEVV3UKT+L2bZpjsXjTv/Z6xAmifPz2WJhFc176XcybwHk4GjuyYAGuiXg+ckOROmi8lzuwoewTNcOU7aCZ0uoAhX1pU1W+BF9L0xp/aPm/7CppE/kqa93QWTQ85NH+jc2l6Ui8Bzh4ltsOBbYELO2Zq/lpb9vfACuCi9nPyf/n9M+APUFXfAv43TW/6jcAjuf8zvy8HPtk+3ytJmkXi0n6SJI1fm/StBo6sqvM2wfU+C5xZVV8cs7IeoO3xvRR4ejvJliRpFjHhlSRpDO1w2h/QPFv6tzTDmndvZwiWJEk9yiHNkiSNbX+aCaRuAZ5H89ypya4kST3OHl6pTyQ5EnhlVU1kQpyNuc6BwGeqalyz1EqSJEndYg+vNMMkOSDJ/0tye7uu5/eS/GFVnTbdya4kSbPFSO3tNF7v50meM13nl2arOd0OQNL4Jdka+ArN0ixn0syU+jTGMVusJEkan15sb5PMqar13bq+NFPZwyvNLI8GqKrTq+q+qvpNVX2jqn6S5KgkFw5WTHJQkmvab6b/K8kFSV7dlh2V5MIk70myLsm17fqgg8ceneSqJHcmWZnktZv+rUqS1DVjtbffS/KfbRt7dZJnDx6Y5KFJTklyY5Lrk/xLu971YPlrOtrYK5M8KcmngV2AL7dLdP1dksVJKsmrkvwC+HaSA5Os7gy0s2c4yTuSfD7JZ9rzX5bk0UnemuTmJNclcTSYZhUTXmlm+SlwX5JPJjkkybbDVUqyHc3amG8FFgDXAH80pNpT2v3b0awvekrHWp83A4cBWwNHAycmedJUvxlJknrUWO3tU4CVNG3o8cDZSR7Wln0SWA88CtgHOAgY/ML5JcA7aNaz3hp4PrC2qv4c+AXwvKraqqre3XGtZwCPAZ47ztifB3yaZp3rH9Gsjb0ZsBNwAvCRcZ5H6gsmvNIMUlV3AAcABXwUWJPknCTbD6l6KHBFVZ3dDn/6AHDTkDqrquqjVXUfTeO8A7B9e52vVtXPqnEB8A2aoVySNGFJTm17ly6fovN9PcltSb4yZP93k/y4/bkhiWsXa6OMo729GXhfVd1bVZ+j+QL5T9ryQ4C/rqpftWs/nwgsaY97NfDuqrq
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "signal_constant = np.ones(int(n_samples))\n",
    "\n",
    "fig, axes = plot_signal_and_spectrum(signal_constant, sample_rate, \"Constant signal\")\n",
    "axes.flat[1].set_xlim(0,10);\n",
    "\n",
    "fig, axes = plot_signal_and_spectrum(signal_constant, sample_rate, \"Constant signal (masked zero)\", ft_kwargs={'mask_bias':True})\n",
    "axes.flat[1].set_xlim(0,10);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FTs: noise (normal distributed)\n",
    "\n",
    "A FT of normal distributed noise will not have easily identifiable peaks in the spectrum.\n",
    "Instead, the spectrum is on average a bit flat, with only low power in any individual peak."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEjCAYAAAAyt1xdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeZgdZZX/v+f2kiZka8KahHQI+y4kQHQQHVdURkVACIwKggzOjOhPHQcdF8RdEBBFIAEEBcKugixhC0tIOksHsi8knXTS2dPpdJZOd9/l/P6oeuu+VbfWe+tu3efzPP30vXWr3vfU9r7vec95zyFmhiAIgiAIgiAIgiBUKolyCyAIgiAIgiAIgiAIfojiKgiCIAiCIAiCIFQ0orgKgiAIgiAIgiAIFY0oroIgCIIgCIIgCEJFI4qrIAiCIAiCIAiCUNGI4ioIgiAIgiAIgiBUNKK4CoIgCCWDiF4goq8Uodz/IKLb4y63GBARE9ExIfe9kYgeMj+PJaK9RFQTkxx3E9GPzM8fJqL2OMrNU5YHiOjn5ufTiGhWuWQRBEEQKhNRXAVBEITQENE6ItpKRAdq264hotfDHM/Mn2LmB2OWqR7ADwHcHGe5lQYzr2fmIcyc9tuPiK4kopkhyruOmX8Wh2xRlPEgmHkRgF1E9G9xlCcIgiD0D0RxFQRBEKJSC+Cb5RZC43MAVjDzxrgLJqLauMusBOKy2haRhwH8R7mFEARBECoHUVwFQRCEqNwM4LtENMLtRyL6ABHNI6Iu8/8HtN9eJ6JrzM/HENEb5n47iOgxbb8TiOhlItpJRCuJ6Is+8nwKwBvaseNMC+BXiGi9Wfb/ab8PIqLbiWiT+Xc7EQ0yf/swEbUT0f8S0RYAf9a2fY+IthHRZiL6PBF9mohWmTL+QCv/bCKaTUS7zH3/aFqFAyGio8xrsoeIXgZwsMt51ZrfrySiVnPftUR0BRGdCOBuAO833Yp3mfs+QER3EdHzRLQPwL/q7rlaHT8wr9c6IrrC7b5pdc80P79pbl5o1nmpuf0CInrXvA6ziOg07fgziGiBKftjABocl+J1AB9V90UQBEEQRHEVBEEQojIfhmLxXecPRHQQgOcA3AFgJIBbATxHRCNdyvkZgJcANAIYA+APZhkHAngZwCMADgUwGcCfiOhkD3lOBbDSZfu5AI4H8FEAPzaVOgD4PwCTALwPwOkAzobhaqw4HMBBAJoAXKttawAwGsCPAUwF8O8AJgD4oFn+eHPfNID/B0PpfL9Z/396yO7kEQAt5rE/A+C6Hti8RncA+BQzDwXwAQDvMvNyANcBmG26FeuTC5cD+AWAoQDcXIkPN+sdbdY7hYiODxKYmc8zP55u1vkYEZ0J4H4YVtORAO4B8Iw5aVAP4O8A/grjOj8B4CJHmRsBJGHcP0EQBEEQxVUQBEHIix8D+AYRHeLY/hkA7zHzX5k5xczTAKwA4LZeMQlDORzFzD3MrJSpCwCsY+Y/m2UsAPAUgIs9ZBkBYI/L9p8y835mXghgIQwlFQCuAHATM29j5u0AfgrgS9pxGQA/YeZeZt6vyfoLZk4CeBSGgvd7Zt7DzEsBLAVwGgAwcwszN5uyr4OhtH3IQ3YLIhoL4CwAPzLrfhPAsz6HZACcQkQHMPNmUw4//sHMbzNzhpl7PPZRdb8BYwLCz9Ltx9cA3MPMc5g5ba5r7oUxYTAJQB2A25k5ycxPApjnUsYeGPdWEARBEERxFQRBEKLDzEsA/BPADY6fRgFoc2xrg2HFc/I9AARgLhEtJaKvmtubAJxjupjuMt1dr4BhEXSjE4YV0ckW7XM3gCEeMraZ2xTbXRS7Di0oklJmt2q/71flE9FxRPRPItpCRLsB/BKay68PowB0MvM+h2w5mPtcCsO6upmIniOiEwLK3xDwu1vdo7x2DqAJwHcc9/BIs7xRADYyMzvqcjIUwK486xcEQRD6GaK4CoIgCPnyExiWNV0p3QRDadEZCyAncBIzb2HmrzHzKBgupX8yI9NuAPAGM4/Q/oYw89c95FgE4LgIcjtlHGtus0SLUJYbd8GwMh/LzMMA/ACGgh7EZgCNpEVsNmVzhZmnM/PHARxh1jdV/eR1SED9bnWr67IPwGDtN69JBMUGGBZq/R4ONi3wmwGMJiL9mtjOk4hGAaiHuwu4IAiCMAARxVUQBEHIC2ZeDeAxANdrm58HcBwRXU5EtWagnpNgWGdtENElRDTG/NoJQ7FKm/seR0RfIqI68+8sbY2qk+cRwhVXYxqAHxLRIUR0MAy354ciHB/EUAC7Aew1raBeCrcNZm6DsX74p0RUT0Tnwt3FGkR0GBF91lQ0ewHshXHtAMMSPCZsQCgHqu4PwnDZfsLc/i6ALxDRYHNy4WrHcVsBjNe+TwVwHRGdQwYHEtFniGgogNkAUgCuN5+RL8BYZ6zzYQCvMXNvHucgCIIg9ENEcRUEQRAK4SYAlpWOmTtgKDzfAdABwx34Ambe4XLsWQDmENFeAM8A+CYzr2XmPQA+AeAyGBa/LQB+A8ArwuyzAE4wrXRh+DkMBXERgMUAFpjb4uK7MAIh7YGhwD3mv7uNywGcA2AnDIv2Xzz2S8C4xpvMfT+EbACo12Csud1CRG7X3YstMCYQNsFIR3MdM68wf7sNQB8MBfVB83edGwE8aLoFf5GZ58Owxv/RLHM1gCsBgJn7AHzB/N4Jw+X5aUd5V8CIjiwIgiAIAACyLzERBEEQhOqDiK4FcBIzf6vcsgiFQUSnApjCzO8vtyyCIAhC5SCKqyAIgiAIgiAIglDRiKuwIAiCIAiCIAiCUNGI4ioIVQIRXUFEL5Wgng8TUXux6xEEQRAEQRCEsIjiKggVBhGdS0SziKiLiHYS0dtEdBYzP8zMnyi3fIIgCILQX/Dqc4tY3zoi+lixyheE/kxtuQUQBCELEQ2DkQrk6wAeh5HH8IMw0l0IgiAIghATldjnElEtM6fKVb8gVDJicRWEyuI4AGDmacycZub9zPwSMy8ioiuJaKbakYg+QUQrzVniPxHRG0R0jfnblUQ0k4huIaJOIlpLRJ/Sjr2KiJYT0R4iaiWi/yj9qQqCIAhCWQnqc98moj+Y/ewKIvqoOpCIhhPRfUS0mYg2EtHPiahG+/1rWj+7jIjOJKK/AhgL4Fki2ktE3yOicUTERHQ1Ea0H8Jrbkh3dUktENxLRE0T0kFn+YiI6joi+T0TbiGgDEYmHltDvEMVVECqLVQDSRPQgEX2KiBrddiKigwE8CeD7AEYCWAngA47dzjG3HwzgtwDuIyIyf9sGI9fmMABXAbiNiM6M+2QEQRAEoYIJ6nPPAdAKox/9CYCniegg87cHAaQAHAPgDBi5p9Xk8SUwcht/GUY/+1kAHcz8JQDrAfwbMw9h5t9qdX0IwIkAPhlS9n8D8FcAjQDeATAdxrh+NIz82veELEcQqgZRXAWhgmDm3QDOBcAApgLYTkTPENFhjl0/DWApMz9tuhTdAWCLY582Zp7KzGkYHewRAA4z63mOmdewwRsAXoLhHiUIgiAIA4IQfe42ALczc5KZH4MxGfwZ8/dPAfgWM+9j5m0AbgNwmXncNQB+y8zzzH52NTO3BYhzo1nW/pDiv8XM080xwBMADgHwa2ZOAngUwDgiGhGyLEGoCkRxFYQKg5mXM/OVzDwGwCkARgG43bHbKAAbtGMYgDMS8Bbt927z4xAAMGeWm81AFLtgKMIHx3smgiAIglDZBPS5G83+VdFm/t4EoA7AZiLaZfaj9wA41NzvSABrIoqyIXgXG1u1z/sB7DAnqtV3wOzzBaG/IIqrIFQwzLwCwAMwOlOdzQDGqC+mC/AYhICIBgF4CsAtAA5j5hEAngdAvgcKgiAIQj/Gpc8drS2xAYz1qZtgKJm9AA5m5hHm3zBmPtncbwOAo72qCbF9H4DB6ou5dvaQKOciCP0RUVwFoYIgohOI6DtENMb8fiSAyQCaHbs+B+BUIvo8EdUC+C8Ah4esph7AIADbAaTMoE0SxEEQBEEYUITocw8FcD0R1ZnrVk8E8Dwzb4axxOZ3RDSMiBJ
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "signal_normal_noise = rng.normal(size=n_samples)\n",
    "\n",
    "plot_signal_and_spectrum(signal_normal_noise, sample_rate, \"Noise (normal distributed)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FTs: sine wave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(95.0, 105.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA74AAAEjCAYAAAAVGd21AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZxeRZU+/pxshITsCWEJJAESIWYDIoYhYZEBwuKCoBJRlhFRWaIzIuIIyjgwKD9GEb/KCLIpGBAEVJaAgiELBJJAIPuekJCtu7Pvyfue3x/3vu9by6l675vu5O1OzvP5JN11b9Wpc+vUvV2nnlNVxMxQKBQKhUKhUCgUCoVif0WzaiugUCgUCoVCoVAoFArF3oQ6vgqFQqFQKBQKhUKh2K+hjq9CoVAoFAqFQqFQKPZrqOOrUCgUCoVCoVAoFIr9Gur4KhQKhUKhUCgUCoViv4Y6vgqFQqFQKBQKhUKh2K+hjq9CoVAoFAqFQqFQKPZrqOOrUCgUCoVCoVAoFIr9Gur4KhQKhaIsiGgmEZ3ZQLKWENG/7kG5jxHRe0S0iYhGNYQuBxqIqBsR/Z2I1hHRQ9XWR6FQKBSKfQV1fBUKhUIBACCiYUT0JhFtIKK1RDSRiD4BAMz8cWYeW2UVbwYwlpnbMfN9Vdal0SHjhMIPAMxn5k7M/LV61NWZiJ4joi1EtJSIvryn+cvJynD/BiKaQkQ7iOhR594SItpJRF2d69OIiImol5P3X518VxHRhIzNolAoFIpGDHV8FQqFQgEiag/gBQC/AtAZwJEA/gvAjmrq5aAngJnSDSJqsY91aar4VwBPN4CcXwPYCaA7gMsB3E9EH9/D/OVklbu/AsAdAB4O1L0YwMhCgogGADg4wzMqFAqFYj+COr4KhUKhAIC+AMDMo5k5x8zbmPlVZv4AsNmw9PebiOiDlB1+iohaFwQR0UlGSPLT6f07pEqJ6Agi+jMR1RDR4lAIMxG9DuAsAP+PiDYTUd9Uj+8T0QcAthBRi3LyiOhEIno31e0pInqyoFvKAB5n5H3U1DsmO0ObHEVEz6Zl64jo/6XXv0dEf3Z0/BUR3Su0wS1EtDDVfRYRXWzc+wOAowH8LW2fm52yrYhoA4ABaZ7pUjtnARG1BXAJgNuYeTMzTwDwVwBfrTR/OVlZ6mLmZ5n5eQB1AZX/AOAKI30lgN/vwXN/KW3bwr8dRDS2UjkKhUKhqA7U8VUoFAoFAMwDkCOix4jofCLqVCb/FwGMANAbwEAAVwGJgwXgOQCPImGORwO4WBJARM0A/A3A+0gY5rMBfIeIznPzMvOnAIwHcAMzH8LM89JbIwFcCKAjgHxMXqrb80gcoc5ImM9LyjxnJbqG2qQ5EjZ9KYBeafkn0zKPAxhBRB3TvC0AfCnV0cVCAMMBdEDCxj9ORIen7fNVAB8C+HTaPnc77bcTwKkA1qT3BzjP9wIRrQ/8e8HRoy+AnGEDpO0SYnxj+cvJqrQuCZMAtCeiE1JbfAlJu1cEZn4qbbtDABwBYBGS/q1QKBSKJgB1fBUKhUIBZt4IYBgABvAggBoi+isRdQ8UuY+ZVzDzWiQO4eD0+lAALdL7u5j5WQDvBGR8AkA3Zv4JM+9k5kVp3ZdVoPp9zLyMmbdlkDcUQEsA96a6PQNgcsZ6sugaapNTkDhK32PmLcy8PWUuwcwrAYwD8IU07wgAtcw81VWAmZ9O5eeZ+SkA81PZWTEYidPogZkvYuaOgX8XOdkPAbDBubYBQLtAvbH85WRVWlcIBdb3HABzAHwUyPe86fQD+I2bIZ0E+SOS9ea/rVAPhUKhUFQJuiZKoVAoFAAAZp6NEkt5PBJW7F4Y6yMNrDJ+34rEsUP68yNmZuP+skCVPQEckToYBTRHwuxmhSm7nDxJt6UZ68mia6hNjgKwlJl3B2Q/BuBbSBzpr0Bme0FEVwD4DySsMZA4hV2lvAEEHd8KsRlAe+daewCb9iB/OVmV1hXCH5BMMPRGPMz5c8z8j0KCiK4CcI2T504kjrfuLK5QKBRNCMr4KhQKhcIDM89BEq7cv8KiKwEcSURkXDsqkHcZgMUOu9iOmS+oRNUK5Em6HW38vhVAGyN9WAPpugzA0RTegOt5AAOJqD+AiwA84WYgop5IHOMbAHRh5o4AZgAwn4Xdcg4GIeD4EtHLzvpV89/LTvZ5AFoQUR9HtrjxWJn85WRVWpcIZl6KZJOrCwA8W0lZE0R0GZKJoEuZedeeylEoFArFvoc6vgqFQqEAER1PRN8loh5p+igkA/xJFYp6C0AOwA3pZlOfRTgc9x0AG9MNqg4mouZE1J/SI5T2AOXkvQVgN4BRqW6fd3SbBuDLabkRAM5oIF3fQeJ0/5SI2hJRayI6rXCTmbcDeAZJ+Ow7zPyhIKMtEse2BgCI6Gr4kxKrARwT0SPo+DLz+YX1q8K/8528W5A4jz9Jn+c0AJ9FgKmO5S8nK0tdqS1bI2Hgm6ftK00yfA3Ap1KZFYOITkSy6/nnmLlmT2QoFAqFonpQx1ehUCgUQBI6+kkAbxPRFiQO7wwA361ESLqJ0ueROBnrkYTuvgDhWCRmzgH4NJIQ3MUAagH8DsnmTRWjnDxDt6sArEOyyZHJ/n07Lb8eybE5zzeErkbZ45BsQLU8rdvEY0h2XA45j7MA/C8S5311mneik+0uALem61NvMm8Q0WEAOiFZ39oQuA7JkUBrkGzw9C1mLrKwKYP8nxnzR2VluH8rgG0AbkHS37al1yww80JmnrLHT5w43J0ATIiw4QqFQqFopCB7qZNCoVAoFA0LInobwP8x8yPV1sUFET0KYDkze47SPtbjaCRO6WHpRmMKhUKhUCgaEMr4KhQHCIjociJ6dR/UcyYRLd/b9SgaL4joDCI6LA1BvRLJ0T5jqq1XY0W6S/B/AHhSnV6FQqFQKPYO1PFVKPYzENEwInqTiDYQ0VoimkhEn2DmJ5j53Grrpzgg8DEka0k3IAmVvjQ9tkfhgIjaAtiI5JidH1dZHYVCUQFCf2/3Yn1LiOhf95Z8hWJ/hx5npFDsRyCi9kjWU34LwJ8AtAIwHML6SoVib4GZHwDwQLX1yAJmvqrK9W9BciyRQqFoQmiMf2+JqEXk2DSF4oCHMr4Kxf6FvgDAzKOZOcfM25j5VWb+gIiuIqIJhYxEdC4RzU1nqn9DRG8Q0TXpvauIaAIR3UNE64hoMRGdb5S9mohmE9EmIlpERN/Y94+qUCgUCkXVUO7v7UQi+lX6N3YOEZ1dKEhEHYjoISJaSUQfEdEdRNTcuP9142/sLCI6iYj+gOT4tb+lG6vdTES9iIiJ6GtE9CGA16XlRiZTTES3E9HTRPR4Kn86EfUloh8Q0RoiWkZEGh2m2C+hjq9CsX9hHoAcET1GROcTUScpExF1RXJ8yg8AdAEwF8C/ONk+mV7vCuBuAA8RFc8/XYPkvNH2AK4G8AsiOqmhH0ahUCgUikaKcn9vPwlgEZK/oT8G8CwRdU7vPYbkaLXjAJwI4FwAhYnnLwC4HcAVSP7GfgZAHTN/Fcmu8J9Ojxm726jrDAAnADgvo+6fRrKDfCcA7wF4BYlPcCSAnwD4bUY5CkWTgjq+CsV+hHRjnGFIzvt8EEANEf2ViLo7WS8AMJOZn03Dou4DsMrJs5SZH0yPYnkMwOEAuqf1vJgeDcLM/AaAV5GEeCkUCoVCsd8jw9/bNQDuZeZdzPwUkonkC9P75wP4DjNvYeY1AH4B4LK03DUA7mbmyenf2AXMvLSMOrensrZlVH88M7+S/v1/GkA3AD9l5l0AngTQi4g6ZpSlUDQZqOOrUOxnYObZzHwVM/cA0B/AEQDudbIdAWCZUYaRnC1qYpVxf2v66yEAkM5uT0o381iPxJHu2rBPolAoFApF40WZv7cfsX1m6NL0fk8ALQGsTM/cXo+EYT00zXcUgIUVqrKsfBYLq43ftwGoTSe5C2lA9x5Q7IdQx1eh2I/BzHMAPIrkD7KJlQB6FBJpCHMPZAARHQTgzwDuAdCdmTsCeAkARQsqFAqFQrGfQvh7e6SxPAhI1ueuQOKk7gDQlZk
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = 1e2 # Hz\n",
    "if f > sample_rate/2:\n",
    "    print(\"Sampling a frequency above the sample_rate/2 gives problems\")\n",
    "\n",
    "signal_sine = np.sin( 2*np.pi*f*t )\n",
    "\n",
    "if False:\n",
    "    # aliased peak\n",
    "    f_alias = (1.1* sample_rate/2)\n",
    "    signal_sine += 1/2*np.sin( 2*np.pi* f_alias *t)\n",
    "\n",
    "fig, axes = plot_signal_and_spectrum(signal_sine, sample_rate, \"Single frequency at $f = {:g}$MHz\".format(f/10**6))\n",
    "axes.flat[1].axvline(sample_rate/2, color='r', label='Nyquist Frequency');\n",
    "axes.flat[1].legend();\n",
    "axes.flat[1].set_xlim(f-5, f+5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sine wave will give very clear peaks at the frequency of the sine wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA80AAAHgCAYAAACfG480AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdaVyTd77//9eVhCUkCAiCIiIqBBVUUFCxblVRqwhVSRftaF2mc6adOZ52Zs50Tud05j+nZ9pZzuOc2lrn8Zvu1W4X7kup1qVabVXc60JAXBAXXCkJa5LrfyOWal1qFUjAz/NOG7iSfEKr5J33dX2/iqZpCCGEEEIIIYQQ4no6bw8ghBBCCCGEEEL4KgnNQgghhBBCCCHETUhoFkIIIYQQQgghbkJCsxBCCCGEEEIIcRMSmoUQQgghhBBCiJuQ0CyEEEIIIYQQQtyEwdsDNKeIiAgtLi7O22MIIYQQTa6qzglAkP899ateCCHEPW7nzp3nNU1r15iPeU/9Jo2Li6OgoMDbYwghhBBNbufxiwD069zWy5MIIYQQzUdRlOON/ZhyerYQQgghhBBCCHETEpqFEEIIIYQQQoibkNAshBBCCCGEEELchIRmIYQQQgghhBDiJiQ0CyGEEEIIIYQQN3FPrZ4thBBC3CuKK77m8OU96I1DSYlM8fY4QgghRIslTbMQQgjRyqw9tpY/7/lX8o7+kxn5M9hYutHbIwkhhBAtljTNQgghRCtQ765nU+kmVJvKllNbGr7u1JzMWT+HEbEjsFqsDIweiE6Rz8yFEEKI2yWhWQghhGjBTtlPsahoEUuKlnCu+hxRQVFMip/EipKVuNxODHoDo2JH8eWpL/nsxGfEmGOYbJnMg/EPEmGM8Pb4QgghhM+T0CyEEEK0ME63k80nN6PaVL4o+wKAITFDeN7yPIM7DsagM9AzeCSHL+8hp7vnmuY6Vx2fHf8M1aby8q6XmbdnHiM6jcCaaKV/+/7SPgshhBA3oWia5u0Zmk1aWppWUFDg7TGEEEKIO3LGcYbFRYtZXLSYs1VnaWdsx6SESUxOmEwHc4drjt15/CIA/Tq3ve5xSipKyLPlsfzIcipqK4gNjsVqsZITn0NYYFizvBYhhBCiKSiKslPTtLRGfUwJzUIIIYTvcrldbDm1BbVQZVPZJjRNY1D0IKyJVobFDMOgu/FJY7cKzd+qddWy5tga8mx57CrfhZ/Oj1GdR2G1WEmLSkNRlCZ5TUIIIURTaYrQLKdnCyGEED6ovKq8oVU+7ThNeGA4M5NnMjlhMjHBMY3yHAH6ACZ0m8CEbhMovlRMXpGnff7k6Cd0CelCbkIuOfE5hASENMrzCSFES1VfX8/Jkyepqanx9ijiisDAQGJiYvDz82vy55KmWQghhPARbs3Nl6e+RLWpbCzdiEtzMbDDQKwWK/d3uh8//e2/MbidpvlGqp3VrDm2ho9tH7Pv3D78df6MjhuN1WIlNTJV2mchxD3p6NGjBAcHEx4eLn8P+gBN07hw4QKVlZV06dLlmu9J0yyEEEK0Querz7O0eCl5tjzK7GW0DWzLtKRp5CbkEtsmtllnMRqM5MTnkBOfQ+HFQlSbyqqSVawsWUm3kG5YE61kdc2S9lkIcU+pqakhLi5OArOPUBSF8PBwzp071zzPJ02zEEII0fzcmpttp7eh2lQ2nNiAU3PSv31/rBYrI2JH4K/3v6vHv9Om+Uaq6qvIP5aPWqjy9YWvCdAHMCZuDFaLlT7t+sibSCFEq3fo0CF69Ojh7THE99zov4s0zUIIIUQLd6H6AsuOLCPPlkdpZSkhASFM7TGVXEsucSFx3h7vhoL8gpiUMIlJCZM4dOFQQ/u8/MhyEsISsFo87XOwf7C3RxVCiFZr7ty5zJ8/n759+7Jw4cI7fpz169fz61//mrq6Ovr168cbb7yBwWDg0qVLzJw5kyNHjhAYGMibb75JcnIyAGazGbvd3vAYb7/9NgUFBbz66qt3/bpaAgnNQgghRBPTNI0dZ3ag2lQ+O/EZTreTvpF9eTLlSTI7ZxKgD/D2iLetR3gPns94nl+l/YrVR1ejFqr8eduf+d+d/8vYuLFYLVaSI5KlfRZCiEb22muv8cknn1x3De+P4Xa7mT59OuvWrcNisfD888/zzjvvMGvWLP785z+TkpLCkiVLOHz4ME899RTr1q1rxFfQckloFkIIIZrIpZpLLD+ynDxbHse+OUawfzCPJD5CriWXbqHdvD3eXTH5mbBarFgtVg6cP4BqU1l9dDVLipfQvW13rBYr47uOx+Rn8vaoQgjR4v3Lv/wLJSUlZGdnM3PmTJ5++uk7epwLFy4QEBCAxWIBIDMzkxdffJFZs2Zx8OBBfve73wHQvXt3jh07xtmzZ4mKirrlY6akpDT8e2FhIfn5+QwbNuyO5vNVEpqFEEKIRqRpGrvKd/Fx4cesPb6Wenc9Ke1S+O/B/83ozqMJNAR6e8RGlxSRRFJEEr9O+zWrSlah2lT+66v/4n8K/odxXcdhtVjpGd7T22MKIUSjsJ2tpLKmvlEfMzjQD0vUzS9x+cc//kF+fj4bNmwgIiLimu8VFhby8MMP3/B+GzduJDQ0tOF2REQE9fX1FBQUkJaWRl5eHqWlpQD06dOHxYsXM3jwYLZv387x48c5efIkUVFRVFdXXxOOL168SHZ2NgB79uwBYMWKFfz1r39l0KBBd/ZD8GESmoUQQohGUFFbwYojK1BtKiUVJQT7BZNrycVqsZIQluDt8ZqF2d/Mw90f5qHEh9h3fh9qocrKIyvJs+WRFJ6E1WLlgS4PEOQX5O1RhRCi1UhMTGwIrj9EURQ+/PBDnn76aWpraxk9ejQGgycSPvvss8yZM4eUlBR69epFampqw/eMRuM1z/HtNc3fKioq4je/+Q3r169vln2Tm5vPhmZFUcYCLwN64HVN0166yXG5gAqka5omS2MLIYRoNpqmsffcXlSbyqfHPqXWVUvviN78adCfGNtlLEaD0dsjeoWiKPRp14c+7frw7/3/nRVHVpBny+OPX/6RvxX8jayuWVgtVhLbJnp7VCGE+NFu1Qh7w49pmgEyMjLYvHkzAGvWrMFmswHQpk0b3nrrLcDz+61Lly63df20w+HgoYce4p///CfR0dF381J8lk+GZkVR9MA8IBM4CexQFGW5pmkHv3dcMPCvwLbmn1IIIcS96pu6b1h5ZCWqTaX4cjEmPxMPxj8oQfAG2vi3YWqPqUzpPoU95/agFqosKVrCR4Uf0TuiN7mW3Hv6AwYhhLhbP6ZpBigvLycyMpLa2lr+8pe/8NxzzwFw+fJlgoKC8Pf35/XXX2fo0KG0adPmBx9vxowZzJgxgyFDhtzxa/B1Phmagf5AsaZpJQCKonwI5AAHv3fcfwF/BX7dvOMJIYS412iaxv7z+1FtKvlH86lx1ZAUnsQfM/4opxzfBkVRSI1MJTUyld/2/y3LjyxHtak8v/V5/rbjb2R1y7qnTmUXQghv+dvf/sbKlStxu938/Oc/Z8SIEYBnz+Np06ah1+vp2bMnb7zxxg8+1vHjx8nLy8Nms/Hmm28C8Prrr5OW1qjbJHudommat2e4zpVTrsdqmjb7yu2fAAM0TfvFVcekAr/XNG2yoigbgV//0OnZaWlp2tXn3gshhBA/xF5nb1jcqvBSIUaDkfFdx/v84lY7j18EoF/ntl6e5OY0TWPn2Z2oNvWaRdOsidZWu2iaEKJlOnToED169PD2GOJ7bvTfRVGUnZqmNWpq99Wm+UabOzake0VRdMD/Ao//4AMpyhPAEwCxsbGNNJ4QQojW7uptlKqd1fRo24P/HPifso1SI1IUhbT2aaS1T+PZmmcbtud67ovn+Mv2v5DdLRurxUrX0K7eHlUIIcQ9zFdD80mg01W3Y4BTV90OBpKBjYqiALQHliuKkv39tlnTtP8H/D/wNM1NObQQQoiWraq+ilVHV6EWqhy6eAijwcjYuLFYLVaSI5K58jtHNIGwwDCmJ01nWs9p7DizA9Wm8mHhhyw4tIB+Uf2wWqxkds7EX+/v7VGFEELcY3w
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with_spectrum = True\n",
    "phase_offsets = np.linspace(0, 2*np.pi, 8, endpoint=False) # rad\n",
    "f = 99 # Hz\n",
    "\n",
    "if f > sample_rate/2:\n",
    "    print(\"Sampling a frequency above the sample_rate/2 gives problems\")\n",
    "\n",
    "if with_spectrum:\n",
    "    fig, axes = plt.subplots(2,1, sharex=True, figsize=(16,8))\n",
    "else:\n",
    "    fig, ax = plt.subplots(figsize=(16,4))\n",
    "    axes = np.array([None, ax])\n",
    "    \n",
    "if with_spectrum:\n",
    "    axes[0].set_ylabel(\"Amplitude\")\n",
    "    axes[0].axvline(f, alpha=0.3, label='f = {:g}Hz'.format(f))\n",
    "\n",
    "axes[1].axvline(f, alpha=0.3, label='f = {:g}Hz'.format(f))\n",
    "axes[1].set_ylabel(\"Phase\")\n",
    "axes[1].set_xlabel(\"f\")\n",
    "axes[1].set_xlim( f-2, f+2)\n",
    "axes[1].set_ylim( -1 - 0.2, 1+ 0.2)\n",
    "axes[1].grid()\n",
    "axes[1].yaxis.set_major_formatter(tck.FormatStrFormatter('%g $\\pi$'))\n",
    "axes[1].yaxis.set_major_locator(tck.MultipleLocator(base=0.5))\n",
    "\n",
    "for phase in phase_offsets:\n",
    "    signal = np.sin( 2*np.pi*f*t + phase )\n",
    "    \n",
    "    ft_signal, freqs = ft_spectrum(signal, sample_rate)\n",
    "    \n",
    "    if with_spectrum:\n",
    "        axes[0].plot(freqs, np.real(ft_signal),  '.-', label='Re $\\\\varphi = {}\\\\pi$'.format(phase/np.pi))\n",
    "    \n",
    "    axes[1].plot(freqs, np.angle(ft_signal)/np.pi, '.-', label='$\\\\varphi = {}\\\\pi$'.format(phase/np.pi))\n",
    "\n",
    "# trigger legends\n",
    "if with_spectrum:\n",
    "    axes[0].legend(loc='center right');\n",
    "axes[1].legend(loc='center right');\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a sine wave, the phase of the spectrum at it's frequency is $\\theta( \\varphi = 0 ) = - \\dfrac{\\pi}{2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FTs: step function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV9fn/8deVAWETCDMhIWyZMkRU3K2CqwoOHFX8Wv1qq347re23w9r2q63tz9q6V4d7gIoWtbaKkyBLkKkQCCTssEcgybl+f5wTPEZGQnLOfXLyfj4e58E55/7c933dQblz3ddnmLsjIiIiIiIikkxSgg5AREREREREpL4p2RUREREREZGko2RXREREREREko6SXREREREREUk6SnZFREREREQk6SjZFRERERERkaSjZFdERERERESSjpJdERERERERSTpKdkVEJFBm1tfM5prZDjO7OY7nXWhmp8TrfCIiIhJfSnZFRAQzG21mH5nZNjPbbGYfmtkxkW0rzexrMTz9LcA0d2/l7n+O1UmqX4e7D3D3aTE6V4mZHR2LY0eO39TMHjOzoshDgrlmNvYQ7aeZWZmZ7Yy8lsYqtmrnbWdmL5nZrkisl0Vt625mU81si5mtM7N7zSwtHnGJiEjjoGRXRKSRM7PWwGvAX4B2QDbwK2BvnELIAxbG6VwxZ2ZZQEdgcQxPkwasBk4G2gA/B543s+6H2OdGd28ZefWty8nN7DYzu60GTe8D9gGdgMuBB8xsQGTb/cAGoAtwNOFr+XZd4hIREYmmZFdERPoAuPsz7l7p7nvc/V/uPt/MngBygVcjFcFbAMysq5lNMrONZrYiuvtxpIL6EzNbFKna/dXMMg50YjN7GzgVuDdy/D5m5mbWK6rN38zsN1HH/qGZzY9UoZ+LPraZdTOzyZG4Ss3s3sj3X7mO6pVeMzsqUgHdGunifF61azroeaPa9SKchKYApZEY6r1a6e673P02d1/p7iF3fw1YAQw/0mOa2bWRv7NtZva6mXWsS4xm1gIYD/zc3Xe6+wfAFOCbkSb5wPPuXubu64A3gAFR+19vZv80s/vMbJOZrTGzr0dt72Fmr0W2bTOzt+oSr4iIJB8luyIi8hlQaWZ/N7OxZpZZtcHdvwmsAs6NVAR/b2YpwKvAPMJV4NOB75rZmVHHvBw4E+hJOJn+2YFO7O6nAe/zRdXxsxrEezEwhnCyNBiYCGBmqYQr1EVA90hszx7sOqIPaGbpkWv6F+Gq7E3AU2YWXQE94HmrXc8y4IfAi5HztHf3isNdUCRp23qQ12s12L8T4Z/zoSrkd0QSww+t2lhlM/spcD1wHtABKAF+c7jzHkYfoLLa3+k8vkho7wEmmFlzM8sGxhJOeKsMBo4jnCB3BB4Cfhy1/R/A64Srxp2A2+oYr4iIJBkluyIijZy7bwdGAw48Amw0symRBOpAjgE6uPvt7r7P3Qsj+02IanOvu692983Ab4FL6zHkP7v7msixXyXcBRZgJNAV+FGk8lkWqSbWxCigJXBn5JreJpw4R8d9sPNWNwT4JPoLM/utmb1vZi+aWfPqO7j7Oe7e9iCvcw4VeCRRfwr4u7svOUizHwM9CD8AeJhwhbtnZP+OhB9GXOruy9x9H/AY4b/numgJbKv23TagVeT9u4QT3+1AMTALeDmq7WDCfx9vunsIWFTtWD2BVCA18nf9YR3jFRGRJKNkV0REcPfF7j7R3XOAgYSTxj8dpHke0DW6+gj8lHB1rcrqqPdFkePVl3VR73cTTqoAugFFNamkHkBXYHUkqapSRDg5PNx5qzuacAUTADMbCPR09xOBfwP/dQTxHVCkyv4E4XGxNx6snbvPcPcd7r7X3f8OfAicFdl8OtAE+Djq7/MNqiWq0dVn4Fbg1sNUn3cCrat91xrYEYn7TWAy0ALIAjKB30W1HUT4oUKVgXw54b0c+AawxsKTdbU72PWLiEjjpGRXRES+JFId/Bvh5ALCFd9oq4EV1aqPrdz9rKg23aLe5wJrahHCbiC6+tm5hvutBnIPMUa2+nVEWwN0iyRhVXIJd+etscj+A/lyZfdEwt1tifw5+gD7vW5fzJRc/fV69faRfYxwBbYTMN7dy2sRqgMWed8OeKna32cbdz/lSztEVZ+BOwlXXQ9Vff4MSDOz3lHfDSHc1bod4f9G7o0k4KXAX4kk4GaWT3gSruhZo4cS9XN197fd/XSgf+S4E2tx/SIi0ggo2RURaeTMrJ+Z/cDMciKfuxHuvlsQabKecBfYKh8D283sx2bWzMxSzWygRZYqiviOmeVEqm0/BZ6rRUifAJdFjjuG8Cy9NfExsBa408xamFmGmZ0Qtb36dUSbAewCbjGz9MiY1nOJjPmthWaRV/T9NZMvqqTbCCd6X+LuY6NmSq7+OtiSQg8ARxEeh7znYAGZWVszOzPy80gzs8uBkwhXVgHmAKea2bBI+9Zm9o1IMn3E3H0X4crt7ZG/jxMIV2KfcPdNhCfUuiESU1vgKr6oiA8GPq1WaR9atd3MxplZ70iMrQj/jL/UdVxERETJroiI7ACOBWaY2S7CSe4C4AeR7XcAP4t0V/2hu1cSTgSPJpywbAIeJbwETpWnCU/2VBh51Wayo/+JHH8r4a6qLx+6eVhUXL0IT0ZVDFwS1eRL11Ft332EJ2caG7me+4ErDzEG9mAx7AIeBBaZWXHk6y188bNpA2yuzTEPxMzygP8m/HewLqoKfHlk++uRSacA0gn//DcSvrabgPPdfWkk5unA7cAkM9tJuKvwGHc/VCW8pr5NOPnfADwD3ODuVZNojSM84ddGYBlQAXwvsm0wUcmrmbUnXOFfEPlqNOExvzuAqYSrzG/XQ7wiIpJErH7uZSIiImFmthL4lrv/O+hYEoGZDQJ+4u6Xmdl1QFN3/0vQcYmIiCQ7VXZFGhAzu9zM/hWH85wSVZUSkTpw90+BIjN7n/ByTI8HHJKIiEijoGRXJAGZ2Wgz+8jMtpnZ5si6mMe4+1PufkbQ8YlI7bj7T9z9RHcfH+nqLCIJ4mD33Bieb6WZfS1WxxeRLxxsxkoRCYiZtSa8vucNwPOElwQ5EdgbZFwiNeXu3YOOQUSkJhLxnmtmaUe4hJqIVKPKrkji6QPg7s+4e6W773H3f7n7fDObaGYfVDU0szPMbGnkafT9ZvaumX0rsm2imX1gZn8wsy1mtsLMxkbte7WZLTazHWZWaGb/Hf9LFRERCdTh7rkfmtlfIvfZJWZ2etWOZtYmssbzWjMrMbPfmFlq1PZro+6zi8xsmJk9QXhZs1cjk8rdYmbdzczN7BozWwW8faDhRNEVYTO7zcxeMLMnI8f/1Mz6mNlPzGyDma02M/UEk0ZPya5I4vkMqDSzv5vZWDPLPFAjM8sCXgR+ArQnvB7l8dWaHRv5Pgv4PfBY1HIiG4BzgNbA1cDdVUuPiIiINBKHu+ceS3hG+Szgl8BkCy+pBvB3wrOI9yK8NNYZQNUD54uA24ArCd9nzwNK3f2bhGeLPzeytNjvo851MuHlxM6sYeznAk8QXnprLuHlxFKAbMIzrD9Uw+OIJC0luyIJxt23E15Ww4FHgI1mNsXMOlVrehaw0N0nR7o7/RlYV61Nkbs/ElmS5e9AF6BT5Dz/dPflHvYu4WViTozdlYmIiCSWGtxzNwB/cvdyd3+O8APksyPbxwLfdfdd7r4BuBuYENnvW8Dv3X1m5D67zN2LDhPObZFjHXTd7Gred/c3I78DvAB0ILwMVznhNcK7R9awFmm0lOyKJCB3X+zuE909BxgIdAX+VK1ZV2B11D5OeF3RaOuitu+OvG0JEHmCXRCZjGMr4eQ5q36vREREJLEd5p5bUm3N6aLI9jzCa1ivjazdvZVwJbVjpF03YHktQ1l9+CZfsj7q/R5gU+ThdtVniNzzRRorJbsiCc7dlwB/I3wDjrYWyKn6EOmenEMNmFlTYBLwB6CTu7cFpgJ2yB1FRESS2AHuudlRw38gPN52DeHEdC+Q5e5tI6/W7j4g0m410PNgp6nB97uA5lUfImOBO9TmWkREya5IwjG
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXjU5fX38fdJAoSwJywKIcMiiCiKEhEVtVprwbWKIiitWpf6a9W22lrtYn3sorVaq9VacWm1IIhrUXEtYkUJEJRNEFmzgAgk7BCyneePmeAYErJOZjL5vK5rLmbmu51vUCZnzn3fx9wdERERERERkXiSEO0ARERERERERBqbkl0RERERERGJO0p2RUREREREJO4o2RUREREREZG4o2RXRERERERE4o6SXREREREREYk7SnZFREREREQk7ijZFRERERERkbijZFdERGKSmR1uZp+Y2U4zu6kJr/upmX2jqa4nIiIikaFkV0REqmVmI83sIzPbbmaFZvahmR0f2rbOzM6M4OVvBWa5ewd3fyhSF6l8H+5+pLvPitC11pvZ0EicO+was8ysyMx2hR4rInm90DVTzexlM9ttZjlmdlml7X3MbIaZbTWzjWb2sJklRTouERFp2ZTsiohIlcysI/Aa8DcgFegF/D9gXxOFEAA+baJrRZyZdQW6A8ub4HI3uHv70OPwhpzIzO40sztr2O0RoBjoAVwOPGpmR4Zt/zuwCTgUGAqcBvywIXGJiIjURMmuiIhUZyCAu09x9zJ33+vub7v7YjP7N5ABvBqqHt4KYGY9zexFM9tsZmvDhx+HKqi3m9myUIXvn2aWXNWFzWwmcDrwcOj8A83MzeywsH3+ZWa/Dzv3z8xscagK/Vz4uc2st5m9FIqrwMweDr1/wH1UrvSa2RGhaum20BDn8yvdU7XXDdvvMCCP4OduQSiGqFQ2zeza0N/BdjN7w8y6N/B87YAxwG/cfZe7zwamA98N260vMM3di9x9I/AmcGTYOa43s9fN7BEz22JmG8zsW6Ft/czstdD7283snYbEKyIiLYeSXRERqc7nQJmZPW1mo82sS8UGd/8ukAucF6oe3mtmCcCrwCKCVeBvAj8xs2+HnfNy4NtAf4LJ9K+rurC7nwF8wFcVys9rEe9YYBTBxOpo4EoAM0skWKHOAfqEYpta3X2En9DMWoXu6W2CVdkbgclmFl4trfK6le5nFfAz4IXQddLcvbSmGwoleduqebx2kEPvDiWHH4bPPzazXwLXA+cD3YD1wO9riqMGA4GySn9HiwhLZoEHgXFmlmJmvYDRBBPeCkcDJxJMkrsDjwG/CG17BniDYNW4B3BnA+MVEZEWQsmuiIhUyd13ACMBBx4HNpvZdDPrUc0hxwPd3P0udy929zWh48aF7fOwu+e5eyHwB2B8I4b8kLtvCJ37VYLDZQGGAz2Bn7v77lB1cXYtzzkCaA/cE7qnmQQT5/C4q7tuZccAC8PfMLM/mNkHZvaCmaVUPsDdz3X3ztU8zq3mOr8A+hFM6icSrFr3D1Vwfw2Md/dV7l4MPEnw760h2gPbK723HegQ9vp9gsnvDiAfyAZeCdt+NMGf8VvuXg4sC9vWH0gEEkN/dx82MF4REWkhlOyKiEi13H25u1/p7unAUQSTxr9Ws3sA6BlefQR+SbAaVyEv7HlO6HyNZWPY8z0EkzCA3kBObSqpVegJ5IUSsAo5BBPJmq5b2VCCFU8AzOwooL+7nwK8C3y/HvEdwN3nuvtOd9/n7k8DHwJnE6y0twbmhf39vMmBierXKsrAbcBtB6ko7wI6VnqvI7AzdK4E4C3gJaAd0BXoAvwpbP8hBL8oqHAUXyW8lwMXABvM7EkzS63Lz0NERFouJbsiIlIr7v4Z8C+CiQgEK77h8oC1laqPHdz97LB9eoc9zwA21CGEPUB49fOQWh6XB2QcZI5s5fsItwHoHUrYKmQQHP5ba6Hjj+Lrld1TCA7PJfTnyCqOe8O+WlW58uONyvtXwwEjuMjYy5X+fjq5+zcOOCCsogzcQ7DqWl1F+XMgycwGhL13DF8tLpZK8O/94VACXgD8k2ACjpn1BZKA8FWjjyX0s3L3me7+TWBw6LxX1vK+RUSkhVOyKyIiVTKzQWZ2i5mlh173Jjh8Nyu0y5cEh8tWmAfsMLNfmFlbM0s0s6Ms1Koo5Edmlh6qzv0SeK4OIS0ELguddxTBFX1rYx7wBXCPmbUzs2QzOzlse+X7CDcX2A3camatQvNfzyM057cO2oYe4Z+7XfiqqrqdYFL4Ne4+OmxV5cqP0ZX3N7POZvbt0D0mmdnlwKkEK6sfA6eb2XGhfTua2QVmZnW8l8ox7iZYtb0r9PM9mWAl9t+h7VuAtcD/hWLqDFzBV1Xuo4EllarnxwKLzOwiMxsQirFD6Gf2taHgIiIi1VGyKyIi1dkJnADMNbPdBJPcpcAtoe13A78ODW39mbuXEUwEhxJMbrYATwCdws75LMHFntaEHnVZHOnHofNvIzi09ZWD7x4UFtdhBBejygcuDdvla/dR6dhigos5jQ7dz9+B74Wq3LUWSgj/ASwzs/zQ21v56mfTCSisyzmr0Yrgz3RzKN4bge+4+wp3nwPcBbxoZrsIDhMe5e4Hq2zX1g8JJvObgCnA/7l7eNuoiwgu4rUZWAWUAj8NbTuasATWzNIIVu2XEqx2v0/wv8UZBCvMMxshXhERaQGscT7jREREDs7M1gHXuPu70Y4lFpjZEOB2d7/MzK4D2rj736Idl4iISLxQZVekGTGzy83s7Sa4zjfCqk8iEgHuvgTIMbMPCLZjeirKIYmIiMQVJbsiMcjMRprZR2a23cwKQ70yj3f3ye5+VrTjE5HG4e63u/sp7j4mNNRZRJpYdZ+5EbzeOjM7M1LnF5GvVLcypYhEiZl1JNjH8/+AaQRbhZwC7ItmXCIN5e59oh2DiEi4WPzMNbOkerZKE5FKVNkViT0DAdx9iruXufted3/b3Reb2ZVmNrtiRzM7y8xWhL6N/ruZvW9m14S2XWlms83sPjPbamZrzWx02LFXmdlyM9tpZmvM7AdNf6siIiJRVdNn7odm9rfQ5+xnZvbNigPNrFOo9/MXZrbezH5vZolh268N+5xdZmbHmdm/CbYvezXUQuxWM+tjZm5mV5tZLjCzqulE4RVhM7vTzJ43s0mh8y8xs4FmdruZbTKzPDPTSDBp8ZTsisSez4EyM3vazEabWZeqdjKzrsALwO1AGsEelSdV2u2E0PtdgXuBJ8PajGwCzgU6AlcBD1S0JBEREWkhavrMPYHgyvFdgd8CL1mwdRrA0wRXFj+MYLuss4CKL5wvAe4Evkfwc/Z8oMDdv0twVfjzQi3E7g271mnAEQTn8NfGeQRbfHUBPiHYYiwB6EVw5fXHankekbilZFckxrj7DoLtNhx4HNhsZtPNrEelXc8GPnX3l0LDnR4CNlbaJ8fdHw+1XnkaOBToEbrO6+6+2oPeJ9gO5pTI3ZmIiEhsqcVn7ibgr+5e4u7PEfwC+ZzQ9tHAT9x9t7tvAh4AxoWOuwa4193nhz5nV7l7Tg3h3Bk6195ahv+Bu78V+h3geaAbwfZcJQR7gfcJ9bUWabGU7IrEIHdf7u5Xuns6cBTQE/hrpd16AnlhxzjB/qHhNoZt3xN62h4g9A12Vmgxjm0Ek+eujXsnIiIisa2Gz9z1lXpR54S2Bwj2tf4i1KN7G8FKavfQfr2B1XUMJa/mXb7my7Dne4EtoS+3K15D6DNfpKVSsisS49z9M+BfBD+Aw30BpFe8CA1PTqcWzKwN8CJwH9DD3TsDMwA76IEiIiJxrIrP3F5h038gON92A8HEdB/Q1d07hx4d3f3I0H55QP/qLlOL93cDKRUvQnOBu9XlXkREya5IzDGzQWZ2i5mlh173BsYDWZV2fR0YYmbfMbMk4EfAIbW8TGugDbAZKA0tXKWFLEREpEWpxWdud+AmM2sVmod7BDDD3b8gOP3nfjPraGYJZtbfzE4LHfcE8DMzG2ZBh5lZILTtS6BfDaF9DiSb2Tlm1gr4NcHPbRG
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_0 = t[n_samples//4]\n",
    "signal_step = np.heaviside(t - t_0, 1)\n",
    "fig, axes = plot_signal_and_spectrum(signal_step, sample_rate, \"Step function $t_0 = {:g}ns$\".format(t_0 * 10**9), ft_kwargs={'mask_bias':True});\n",
    "axes.flat[1].set_xlim(0, 20);\n",
    "\n",
    "t_0 = t[n_samples//2]\n",
    "signal_step = np.heaviside(t - t_0, 1)\n",
    "fig, axes = plot_signal_and_spectrum(signal_step, sample_rate, \"Step function $t_0 = {:g}ns$\".format(t_0 * 10**9), ft_kwargs={'mask_bias':True});\n",
    "axes.flat[1].set_xlim(0, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FTs: delta peak\n",
    "\n",
    "All frequencies are accounted for with the same amplitude, but the phase will change depending on the position of the delta peak within the signal.\n",
    "\n",
    "Still, the individual phases are linearly related to their corresponding frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZhcZZn///enOztJOishC1nAILLJEgIoS9gDw6IiGHREHBD1q4KjMwozjjCM/ERcUIRRUVBUZJERBQUBgahRtoAgu4QkkM6+dWcP6e7798c53amuVFdX0lWpqu7P67r66jrnPOc+96mqpPqu5zzPUURgZmZmZmZm1p3UlDsBMzMzMzMzs2JzsWtmZmZmZmbdjotdMzMzMzMz63Zc7JqZmZmZmVm342LXzMzMzMzMuh0Xu2ZmZmZmZtbtuNg1MzMzMzOzbsfFrpmZmZmZmXU7LnbNzHoYSfMlnVCkWG+X9DdJayVdXIyY23HsFyVN25nHrATlfv2Kefxy607nYmZm23Kxa2ZWZdI/0DemBUqDpL9K+oSk7f4/vQh/7H8BmBkRgyLiui7EyStXnhGxb0TMLNUxtyeXrO2fljRb0mZJP+mgzRhJ9Z213wnFWN7Xz8WgmZlVMxe7ZmbV6fSIGARMAK4GvgjcVIY8JgAvluG4lWwR8BXg5jxtTgV+vx3tS8Wvn5mZdVsuds3MqlhENEbEPcAHgI9I2g/aeg7/T9JySfNyXaIq6WfAeOBeSeskfSFdf6mk19Oe45ckvTfXsSU9AhwLXJ/uv1e6PiS9LaPdTyR9JX08X9K/Sfq7pEZJd0jql9F2d0m/SvNeKen6PHm29TpKeoekmWlP94uSzsjKNe9xs9p2eP4d5ZL1mvwqIn4NrMwVP3UqcF++9p0c68ACz6XD56Wj168rxy/kfZcR/4uSFqbP86uSjs/Ylu81mC/p39Pjr5d0k6RRku5P2/9B0tCs9pelcVZL+nGe56vg/M3MrPK52DUz6wYi4kmgHjhKyeXM9wLPAWOB44HPSjo5a58PA2+S9BIPjIhr0k2vA0cBdcB/Az+XNDrHMY8D/gx8Ot3/HwWmew4wHZgEHACcDyCpFvgt8AYwMc399jx5ku7XOz3fB4Fdgc8At0p6eyHHzaHD8+8sl0Kk+R4NPJSvXSfH6vRcOnteOnv9tvf4hb7v0tzeDnwaODS9QuFkYH5Gk87eg2cBJwJ7AacD9wP/AYwg+dsmu0j9UHqMPdN9vpQjp4LzNzOz6uBi18ys+1gEDAMOBUZGxJUR8VZEzAV+CMwoJEhE/DIiFkVES0TcAbwGTC1intel8VeRFBcHpuunAmOAf4+I9RGxKSJmFRDvcGAgcHV6vo+QFM3nFnjcdnbC+R8NPBcRa7sQo5BzKfR5Kdbxt+d91wz0BfaR1Dsi5kfE660bC3gNvhsRSyNiIUnB/kRE/C0iNgN3AwdlHe/6iFiQ5ntVB89Bl/7dmJlZ5elV7gTMzKxoxgKrSMZhjpHUkLGtlqQo6JSk84DPkfSuQlIwjShemizJeLyBpMAF2B14IyKatjPeGGBBRLRkrHuD5Pko5Ljt7ITzb7uEuQsKOZdCn5diHb/g911EzJH0WeAKYF9JDwCfi4hFUNBrsDTj8cYcywOzDrkg4/Eb5H6+uvTvxszMKo+LXTOzbkDSoSRFzCySP/TnRcTkAnaNrDgTSHqzjgcei4hmSc8C2o50NgADMpZ3I7nEujMLgPGSeuUoeCPXDqlFwO6SajIKu/FAoZdVtynw/PPlUohTgZzjoHPoyrGK8bxsz/EXUPj7joj4BfALSYOBHwBfAz5cpPdgtt0zHo8neW66lL+ZmVU+X8ZsZlbFJA2WdBpwO/DziHgeeBJYk04A1F9SraT90oI421Jgj4zlXUgKnOVp/I8C+21nWs8CH0yPOx04psD9ngQWA1dL2kVSP0nv7iDPTE8A64EvSOqt5N67p5M8J9urkPPPlwuSeqUTINUCtel59Eq3TQL6RsQrhbTv7FidKMbzsj3HL/h9p+T+vsdJ6gtsIumNbU43F+M9mO1TksZJGkYytveOruRvZmbVwcWumVl1ulfSWpLeqP8EvgV8FCAimkmKmgOBecAK4Eckk/1k+yrwJSWz9f5bRLwEfBN4jKTQ2R/4y3bmdkl6/AaSiYF+XchOGXm/jWRipHqSWaa3yTNrv7eAM4BTSM71f4HzMgvKQhV4/h3mkvoSSfF2KfDP6ePWCZH+iW0vYc7XvrNj5TuXYjwvBR9/O993fUlumbWC5JLoXUmK0EJfg+31C5KJuuamP1/pYv5mZlYFFNHVq7HMzMysEJLuI5ksqatjdq1AkuYDF0bEH8qdi5mZ7Vzu2TWrIpI+JOnBnXCcaZIKGWNpZttnJvBouZMwMzPrCVzsmlUgSUdK+qukRkmrJP1F0qERcWtEnFTu/Mxsx0TENRGxsdx5mNlWHX3mlvB48yWdUKr4ZraVZ2M2qzDpzKS/BT4J3An0AY4CNpczLzOzahQRE8udg1WuSvzM7WBGejPbAe7ZNas8ewFExG0R0RwRGyPiwYj4u6TzJc1qbSjpJEmvpt9G/6+kP0q6MN12vqRZkr4habWkeZJOydj3o5JelrRW0lxJH9/5p2pmZlZWnX3m/kXSd9PP2VckHd+6o6Q6STdJWixpoaSvSKrN2P6xjM/ZlyQdLOlnJLe/ulfSOklfkDRRUki6QNKbwCO5hhNl9ghLukLSLyX9PI3/vKS9JF0maZmkBZJ8JZj1eC52zSrPP4BmSbdIOkXS0FyNJI0A7gIuA4YDrwLvymp2WLp+BHANcJOk1ntVLgNOAwaTzOJ7raSDi30yZmZmFayzz9zDSGbwHgFcDvwqvYUVwC1AE8kM8gcBJwGtXzifDVwBnEfyOXsGsDIiPkwy2/zpETEwIq7JONYxwDuAkwvM/XTgZ8BQ4G/AAyR/248FriS5f7VZj+Zi16zCRMQa4EiS+0z+EFgu6R5Jo7Kangq8GBG/Si93uo7kFh6Z3oiIH6a31LgFGA2MSo/zu4h4PRJ/JLktx1GlOzMzM7PKUsBn7jLg2xGxJSLuIPkC+Z/S7acAn42I9RGxDLgWmJHudyFwTUQ8lX7OzomINzpJ54o0VqHj+v8cEQ+kfwP8EhgJXB0RW0jupz1R0pACY5l1Sy52zSpQRLwcEedHxDhgP2AM8O2sZmNI7rHauk+Q3Jc005KM7RvShwMB0m+wH08n42ggKZ5HFPdMzMzMKlsnn7kLo/19Ot9It08AegOL0/tQN5D0pO6attsdeH07U1nQeZN2lmY83gisSL/cbl2G9DPfrKdysWtW4SLiFeAnJB/AmRYD41oX0suTx1EASX2B/wO+AYyKiCHAfYDy7mhmZtaN5fjMHZsx/AeS8baLSArTzcCIiBiS/gyOiH3TdguAPTs6TAHr1wMDWhfSscAjt+dczMzFrlnFkbS3pM9LGpcu7w6cCzye1fR3wP6S3iOpF/ApYLcCD9MH6AssB5rSias8kYWZmfUoBXzm7gpcLKl3Og73HcB9EbGYZPjPNyUNllQjaU9Jx6T7/Qj4N0mHKPE2SRPSbUuBPTpJ7R9AP0n/JKk38CWSz20z2w4uds0qz1qSCTGekLSe5AP3BeDzmY0iYgVwNsnEUyuBfYDZFHC7hIhYC1xMcpuF1cAHgXuKdwpmZmZVobPP3CeAycAK4Crg/RGxMt12HsmXxy+RfJbeRTI3BhHxy7T9L9Jj/Bpondjqq8CX0suf/y1XUhHRCPw/kqJ5IUlPb/ZQJTPrhNoPQzCzaiWphuSD8EMR8Wi58zEzM6tmks4HLoyII8udi5ntGPfsmlUxSSdLGpKOwf0PkjG32Zc7m5mZmZn1OC52zarbESSzPa4gud/ee7bjlgVmZkUh6aJy51AsPpfK013OA3wulai7nAf4XHLG8WXMZmZm1hWSZkfElHLnUQw+l8rTXc4DfC6VqLucB/hccnHPrpmZmZmZmXU7VdezO2LEiJg4cWK50zAzs27i6aefXhERvn9lFwwcODD23nvvcqdRFMuXL2fkyO7xdugu59J
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde7xVZb3v8c+XhYIGIrdQuQgq5kYtrSXqydQ0FW+Q18BOpWFUW9Nuu3DXVvLkSe1uYmVpedwqoKUbi9S8UFkqgmKKihKgLDBFboqCuOB3/hhj0XQy11oT1hxr3r7v14sXY47xjN/4PZOlc/3mM55nKCIwMzMzMzMzqyVdyp2AmZmZmZmZWam52DUzMzMzM7Oa42LXzMzMzMzMao6LXTMzMzMzM6s5LnbNzMzMzMys5rjYNTMzMzMzs5rjYtfMzMzMzMxqjotdMzMzMzMzqzkuds3MapikxZI+UqJY75H0uKTXJV1Qiphbce15ko7szGtWk3L/O5fy+uVWS30xM6t3LnbNzCpY+ov3urTwWC3pb5I+J2mr//9dgl/ivwbMjIieEXFVB+K0qVCeEbFvRMzM6ppbk8s2xjlf0mxJb0n6dSttdpPU1F77TijG2vx3djFoZmbVwsWumVnlOzkiegK7A5cDXweuK0MeuwPzynDdWrAM+DZwfRttTgDu2or2WfG/s5mZ1QQXu2ZmVSIi1kTEdOBjwKck7QebRwR/I2m5pEWFbj2VdCMwBLhT0lpJX0v3T5T0j3Tk+GlJpxS6tqT7gQ8DV6fn753uD0l75bT7taRvp9uLJX1V0t8lrZE0VVL3nLaDJf02zXuFpKvbyHPzaKKkf5M0Mx3pnidpdF6ubV43r22r/W8tl20REb+NiDuAFW00OwGY0Vb7dnI6oJg+p3EKvoet/Tt35PrF/HzmxP+6pKXpv8d8SUfnHGvr32qxpP9Ir/+GpOskDZD0h7T9vZJ657W/KI2zStKv2vgZKTp/MzOrLC52zcyqTETMApqADym5nflO4AlgIHA08EVJx+Wd8wngRZJR4h4RcWV66B/Ah4BewLeA/5a0a4FrHgX8BTg/Pf+5ItM9ExgFDAPeC5wNIKkB+B3wAjA0zX1KG3mSnrdd2t97gHcDXwBukvSeYq5bQKv9by+XUkr7dTjwx7batZNTUX1u6z1s7995a69f7M9nmtd7gPOBg9I7GY4DFuc0ae9n9TTgGGBv4GTgD8B/Av1Ift/JL1I/nl5jz/ScbxbIqej8zcys8rjYNTOrTsuAPsBBQP+IuDQiNkTEQuAXwNhigkTErRGxLCI2RcRU4HlgZAnzvCqNv5KkaDgg3T8S2A34j4h4IyLWR8SDRcQ7BOgBXJ72936Sonlckdd9h07of7EOB56IiNc7EKOoPlP8e1iK62/Nz+dGoBswQtJ2EbE4Iv7RcrCIf6ufRMTLEbGUpGB/JCIej4i3gNuBA/Oud3VELEnzvayV/nfovy8zMyuvruVOwMzMtslAYCXJ/MrdJK3OOdZA8st+uyR9EvgyyegqJEVQv9KlyT9ztt8kKXABBgMvRETzVsbbDVgSEZty9r1A8n4Uc9136Gj/Jc0Ejmjl8F8j4rAiQ22+hbkDiuozxb+Hpbh+0T+fEbFA0heBScC+ku4GvhwRy6Cof6uXc7bXFXjdI++SS3K2X6Dw+9Wh/77MzKy8XOyamVUZSQeRFCYPkvwCvygihhdxauTF2Z1klOpo4KGI2ChpLqCtSOdNYMec17uQ3GLdniXAEEldCxS8UeiE1DJgsKQuOcXaEKDY26o3K7L/beVCRBy5tddtxQlAwfnShS7bwWt19D3cmusvofifTyLiZuBmSTsBPweuAD5Rop/VfINztoeQvC8dyt/MzCqLb2M2M6sSknaSdBIwBfjviHgSmAW8li7ss4OkBkn7pQVxvpeBPXJev4ukcFmexj8H2G8r05oLnJVedxStj3LmmwW8BFwu6V2Sukv6YCt55noEeAP4mqTtlDx792SS92RrFdP/tnIpmqSu6QJIDUBD2t+u6bFhQLeIeLaY9iXIqaPv4dZcv+ifTyXP9z1KUjdgPclo7Mb0cCl+VvOdJ2mQpD4kc3undiR/MzOrPC52zcwq352SXicZZfoG8APgHICI2EhSqBwALAJeBX5JsohPvu8A31SyAu9XI+Jp4PvAQyQFzP7AX7cytwvT668mWfDnjmJOysl7L5IFj5pIVpneIs+88zYAo4HjSfp6DfDJ3EKxWEX2v9VcttI3SYq3icD/TrdbFkQ6kS1vYW6rfYdyKsF7WPT1t/LnsxvJo7VeJbkl+t0kRWix/1Zb62aSRboWpn++3cH8zcyswiiio3dDmZmZ2baSNINksaSOztm1IklaDJwbEfeWOxczM8uOR3bNqoikj0u6pxOuc6SkYuZdmlnHzQQeKHcSZmZmtcbFrlkFknSYpL9JWiNppaS/SjooIm6KiGPLnZ+ZlU5EXBkR68qdh1m9au0zN8PrLZb0kazim9m/eDVmswqTrkL6O+DzwDRge+BDwFvlzMvMrFZExNBy52CVoRI/c1tZpd7MtoFHds0qz94AEXFLRGyMiHURcU9E/F3S2ZIebGko6VhJ89Nvo6+R9CdJ56bHzpb0oKTvSVolaZGk43POPUfSM5Jel7RQ0mc7v6tmZmZl1d5n7l8l/ST9nH1W0tEtJ0rqJek6SS9JWirp25Iaco5/Judz9mlJ75d0I8mjru6UtFbS1yQNlRSSxkt6Ebi/0HSi3BFhSZMk3Srpv9P4T0raW9JFkl6RtESS7wSzuudi16zyPAdslHSDpOMl9S7USFI/4DbgIqAvMB/4X3nNDk739wOuBK6T1PJcyleAk4CdSFb2/aGk95e6M2ZmZhWsvc/cg0lW6+4HXAL8Nn1cFcANQDPJqvIHAscCLV84nwFMAj5J8jk7GlgREZ8gWYH+5IjoERFX5lzrCODfgOOKzP1k4EagN/A4cDfJ7/YDgUtJnlVtVtdc7JpVmIh4DTiM5JmSvwCWS5ouaUBe0xOAeRHx2/R2p6tIHteR64WI+EX6+IwbgF2BAel1fh8R/4jEn0gewfGh7HpmZmZWWYr4zH0F+FFEvB0RU0m+QD4xPX488MWIeCMiXgF+CIxNzzsXuDIiHk0/ZxdExAvtpDMpjVXsHP6/RMTd6e8AtwL9gcsj4m2S52YPlbRzkbHMapKLXbMKFBHPRMTZETEI2A/YDfhRXrPdSJ672nJOkDyrNNc/c46/mW72AEi/wX44XYxjNUnx3K+0PTEzM6ts7XzmLo13PqfzhfT47sB2wEvpM6dXk4ykvjttNxj4x1amsqT9Ju/wcs72OuDV9MvtlteQfuab1SsXu2YVLiKeBX5N8gGc6yVgUMuL9PbkQRRBUjfgN8D3gAERsTMwA1CbJ5qZmdWwAp+5A3Om/0Ay33YZSWH6FtAvInZO/+wUEfum7ZYAe7Z2mSL2vwHs2PIinQvcf2v6YmYuds0qjqR9JH1F0qD09WBgHPBwXtPfA/tL+qikrsB5wC5FXmZ7oBuwHGhOF67yQhZmZlZXivjMfTdwgaTt0nm4/wbMiIiXSKb/fF/STpK6SNpT0hHpeb8EvirpA0rsJWn39NjLwB7tpPYc0F3SiZK2A75J8rltZlvBxa5Z5XmdZEGMRyS9QfKB+xTwldxGEfEqcAbJwlMrgBHAbIp4XEJEvA5cQPKYhVXAWcD00nXBzMysKrT3mfsIMBx4FbgMOD0iVqTHPkny5fHTJJ+lt5GsjUFE3Jq2vzm9xh1Ay8JW3wG+md7+/NVCSUXEGuDfSYrmpSQjvflTlcysHXrnNAQzq1aSupB8EH48Ih4odz5mZmbVTNLZwLkRcVi5czGzbeORXbMqJuk4STunc3D/k2TObf7tzmZmZmZmdcfFrll1O5RktcdXSZ6399GteGSBmVlJSJpQ7hxKxX2pPLXSD3BfKlGt9APcl4JxfBuzmZmZdYSk2RHRWO48SsF9qTy10g9wXypRrfQD3JdCPLJrZmZmZmZmNafqRnb79esXQ4cOLXcaZmZWI+bMmfNqRPj5lR3
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9e5ykW1nf+33qfulbVffM9G26p2Hf2OARyQajCYpgIqCISTTZoDH4wXAkmsREYyDxJMQjR/RoSIySHBGFiApoPAkm4CUiUQS37C0o7g0b9p6a6st0z/RUVd/qfln5Y73VU9NT3V2XVd3VPc/386lPd72X9a71vG9VrWet33oeMcagKIqiKIqiKIqiKOcJ32lXQFEURVEURVEURVFco86uoiiKoiiKoiiKcu5QZ1dRFEVRFEVRFEU5d6izqyiKoiiKoiiKopw71NlVFEVRFEVRFEVRzh3q7CqKoiiKoiiKoijnDnV2FUVRFEVRFEVRlHOHOruKoiiKoiiKoijKuUOdXUVRlHsMEbkmIl/vqKwHReQzIrIrIv/IRZldXPtJEXnZSV5zGDjt++fy+qfNeWqLoiiKcjfq7CqKopwxvA560XNQtkTkkyLyPSLS9Xe6g87+DwEfN8aMGmN+uo9yjqRdPY0xzzfGfHxQ1+ymLgf2f5+IPC4iZRF57yHHzIrI6oFt94tISUTe3+m1HHDk/VNnUFEURTnLqLOrKIpyNnmNMWYUWATeAfxz4D2nUI9F4MlTuO4wcx34UeAXjjjm1cBvHdj2s8CnB1WpQ9D7pyiKopxb1NlVFEU5wxhjto0xHwb+DvD3ROQFsD9z+F9EZFNEUu0kqiLyS8AC8JsisiciP+Rtf4uIPOvNHD8lIn+j3bVF5GPA1wE/453/gLfdiMh9Lce9V0R+1Pv/moj8oIj8uYhsi8gHRSTScuxlEfkNr94ZEfmZI+q5P+soIs8TkY97M91Pisg3H6jrkdc9cOyh7T+sLgfuyW8YY/4rkGlXvsergY+0lPsosAX8XofXemGHbTnULofdv36u38lz11L+PxeRNc/OT4vIK1r2HXUPronIP/OunxeR94jIJRH5qHf8/xSRxIHj3+qVkxORXzzCXh3XX1EURRl+1NlVFEU5Bxhj/gRYBV4qVs78m8CfAXPAK4DvF5FvOHDO3wWWsbPEI8aYn/B2PQu8FBgH/g3wfhGZaXPNlwN/CHyfd/4XO6zu3wZeCSwB/wfwBgAR8QP/HUgDV7y6f+CIeuKdF/Ta+zvAReAfAr8sIg92ct02HNr+4+rSCV59vwb4Xe/9GPAjwA+0HnfMtY5ty3F2Oe7+dXv9Tp87r24PAt8HvNhTKHwDcK3lkOOewb8F/DXgAeA1wEeBfwFMYfs2B53Ub/eu8VzvnB9uU6eO668oiqKcDdTZVRRFOT9cB5LAi4ELxpgfMcZUjDFXgXcDj3ZSiDHm14wx140xDWPMB4EvAS9xWM+f9srPYp2LF3rbXwLMAv/MGJM3xpSMMZ/ooLy/DIwA7/Da+zGs0/y6Dq97ByfQ/q8B/swYs+u9/7+B9xhjVrooo5O2dGqXXmh3/W6euzoQBh4WkaAx5pox5tnmzg7uwX8wxtwwxqxhHfbHjDGfMcaUgf8f+IoD1/sZY8yKV9+3H2KDvj43iqIoyvAROO0KKIqiKM6YA7LYdZizIrLVss+PdQqORUS+E/in2NlVsA7TlLtqstHyfwHr4AJcBtLGmFqX5c0CK8aYRsu2NNYenVz3Dk6g/fsSZhF5IfD13O2cHUcnbenULr3Q7vodP3fGmGdE5PuBtwHPF5HfBv6pMeY6dHQPbrT8X2zzfuTAJVsHEtK0t1dfnxtFURRl+FBnV1EU5RwgIi/GOjGfwHb0U8aY+zs41RwoZxE7m/UK4FPGmLqIfBaQLqpTAGIt76exEuvjWAEWRCTQxuE17U7wuA5cFhFfi2O3AHQqq96nw/YfVZdOeDXQXIP6MqxDtywiYO+dX0QeNsa8qM9rubBLN9dfofPnDmPMrwC/4sm4/z/gx4G/6+gZPMjllv8XsLbpq/6KoijK8KMyZkVRlDOMiIyJyDcBHwDeb4z5HPAnwI4XACgqIn4ReYHnEB/kBvCclvdxrIOz6ZX/XcALuqzWZ4HXe9d9JfC1HZ73J8A68A4RiYtIRET+yiH1bOUxIA/8kIgExebefQ3WJt3SSfuPqgsiEvACIPmxjmtERALeviUgbIz5gnf4z2HXkb7Qe/0n4H9g15cee61jcGGXbq7f8XMnNr/vy0UkDJSws7F1b7eLZ/Ag3ysi8yKSxK7t/WA/9VcURVHOBursKoqinE1+U0R2sbNR/xL4t8B3ARhj6lin5oVACrgF/Dw22M9Bfgz4YbHRen/QGPMU8FPAp7COzpcBf9Rl3f6xd/0tbGCg/9rJSS31vg8bGGkVG2X6rnoeOK8CfDPwKmxb3wV8Z4tD2TEdtv/Qunj8MNZ5ewvwHd7/zYBI30hLFGZjTMEYs9F8AXtAyRiz2eG1jmqLC7t0fP0un7swNmXWLawk+iLWCe30HnTLr2ADdV31Xj/aZ/0VRVGUM4AY068aS1EURVGUThCRj2CDJX3k2IMVJ4jINeC7jTH/87TroiiKopwsOrOrKGcIEfl2EfmdE7jOy0SkkzWWiqJ0x8eB3z/tSiiKoijKvYA6u4oyhIjIXxWRT4rItohkReSPROTFxphfNsb89dOun6IovWGM+QljTPG066Eoym0O+80d4PWuicjXD6p8RVFuo9GYFWXI8CKT/nfgzcCHgBDwUqB8mvVSFEU5ixhjrpx2HZThZRh/cw+JSK8oSg/ozK6iDB8PABhjftUYUzfGFI0xv2OM+XMReYOIfKJ5oIj8dRF52huNfpeI/C8R+W5v3xtE5BMi8pMikhORlIi8quXc7xKRz4vIrohcFZH/8+SbqiiKoiinynG/uX8kIv/B+539goi8onmiiIyLyHtEZF1E1kTkR0XE37L/77f8zj4lIi8SkV/Cpr/6TRHZE5EfEpErImJE5I0isgx8rN1yotYZYRF5m4j8moi83yv/cyLygIi8VURuisiKiKgSTLnnUWdXUYaPLwJ1EXmfiLxKRBLtDhKRKeDXgbcCk8DTwFcfOOwrve1TwE8A7xGRZq7Km8A3AWPYKL7vFJEXuW6MoiiKogwxx/3mfiU2gvcU8K+B3/BSWAG8D6hhI8h/BfDXgeaA87cBbwO+E/s7+81Axhjzd7HR5l9jjBkxxvxEy7W+Fnget1OPHcdrgF8CEsBngN/G9u3ngB/B5q9WlHsadXYVZcgwxuwAfxWbZ/LdwKaIfFhELh049NXAk8aY3/DkTj+NTeHRStoY824vpcb7gBngkned/2GMedZY/hc2LcdLB9cyRVEURRkuOvjNvQn8O2NM1RjzQewA8jd6+18FfL8xJm+MuQm8E3jUO++7gZ8wxnza+519xhiTPqY6b/PK6nRd/x8aY37b6wP8GnABeIcxporNp31FRCY6LEtRziXq7CrKEGKM+bwx5g3GmHngBcAs8O8OHDaLzbHaPMdg85K2stGyv+D9OwLgjWD/sReMYwvrPE+5bYmiKIqiDDfH/OaumTvzdKa9/YtAEFj38lBvYWdSL3rHXQae7bIqK8cfcgc3Wv4vAre8we3me/B+8xXlXkWdXUUZcowxXwDei/0BbmUdmG++8eTJ83SAiISB/wL8JHDJGDMBfASQI09UFEVRlHNMm9/cuZblP2DX217HOqZlYMoYM+G9xowxz/eOWwGee9hlOtieB2LNN95a4AvdtEVRFHV2FWXoEJGHROQHRGTee38ZeB3wxwcO/R/Al4nIt4hIAPheYLrDy4SAMLAJ1LzAVRrIQlEURbmn6OA39yLwj0Qk6K3DfR7wEWPMOnb5z0+JyJiI+ETkuSLytd55Pw/8oIj8JbHcJyKL3r4bwHOOqdoXgYiIfKOIBIEfxv5uK4rSBersKsrwsYsNiPGYiOSxP7h/AfxA60HGmFvAt2EDT2WAh4HH6SBdgjFmF/hH2DQLOeD1wIfdNUFRFEVRzgTH/eY+Btw
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9e5xkZ1nv+32qalV3VV+qMplM7sNALoQASiCgqGxubiDIxRvsgEfED8L2CLrdXhD2wW2OHzgbOSrKFlQQhINgQLZbggRBBVQUkCDXhARDkkkmyVwzVX2p26qq5/yx1uqu6emuWqtq3abn+X4+85nuqlVvvetdtbre531+7/MTVcUwDMMwDMMwDMMwdhOFrDtgGIZhGIZhGIZhGHFjwa5hGIZhGIZhGIax67Bg1zAMwzAMwzAMw9h1WLBrGIZhGIZhGIZh7Dos2DUMwzAMwzAMwzB2HRbsGoZhGIZhGIZhGLsOC3YNwzAMwzAMwzCMXYcFu4ZhGIZhGIZhGMauw4JdwzCMswwRuUdEfjCmth4pIl8RkVUR+YU42ozw3reKyNPSfM88kPX1i/P9s2Y3nYthGIZxOhbsGoZhnGH4E/S2H6A0RORfRORnRSTy3/QYJvuvBT6rqkuq+rYZ2hnLdv1U1Uer6meTes8ofdny/GtE5BYR6YrIe3c45iIROSQicyLybhE56F/Pr4jIdWHfKwbGXj8LBg3DMIwzGQt2DcMwzkyer6pLwMOANwO/Brw7g348DLg1g/fNMw8AbwTeM+aY5wJ/A5SA+4CnAjXg14EPi8iBZLu4gV0/wzAMY9diwa5hGMYZjKo2VfUm4D8BPyUij4GNzOH/EpFjInL3dhJVEXk/sB/4mIisichr/cdfJyLf8TONt4nIj2z33iLyaeDpwB/4r7/Sf1xF5PKR494rIm/0f75HRH5FRL4uIk0R+ZCIzI8ce6mI/KXf7xMi8gdj+rmRdRSRR4nIZ/1M960i8oItfR37vluO3fH8d+rLlmvyl6r6V8CJ7dr3eS5ws6quq+oNqnqPqg5V9a+Bu4EnTHivx4U8lx3HZafrF/Jct33/MJ+7kfZ/TUTu98f5DhF55shz467BPSLyq/77r/uZ8fNF5BP+8X8nIudsOf71fjsnReRPx4xX6P4bhmEY+ceCXcMwjF2Aqv4rcAh4inhy5o8BXwMuBp4J/KKIPHvLa34SuBcvS7yoqm/xn/oO8BS8TOP/DfyZiFy4zXs+A/gn4DX+678dsrsvBp4DPBz4LuDlACJSBP4aOAgc8Pt+45h+4r/O8c/3U8A+4OeBD4jII8O87zbseP6T+hIGv7//AfjbbZ47H7gSuHXCe008l0njMun6RX3/sJ87v2+PBF4DPNFXKDwbuGfkkEmfwR8D/qM/Vs8HPgH8N2Av3txma5D6E/57XOa/5g3b9Cl0/w3DMIwzAwt2DcMwdg8PAHuAJwLnqepvqmpPVe8C3gVcH6YRVf0LVX3AzzR+CPh34Ekx9vNtfvsP4QUXj/MffxJwEfCrfsazo6qfC9He9wKLwJv98/00XtD8kpDvewopnP9/AL6mqqujD/rB6QeA96nq7RPaCHMuYcdlGrZ7/yifuwEwB1wtIo6f2f5O8GSIa/A/VfWIqt6PF7B/UVW/oqpd4H8D12x5vz9Q1fv8/r5phzGY6b4xDMMw8kcp6w4YhmEYsXEx8BDePsyLRKQx8lwRLyiYiIi8DPglvOwqeAHT3vi6yeGRn1t4AS7ApcBBVe1HbO8i4D5VHY48dhBvPMK87ymkcP7PBW7e8p4F4P1ADy/jOYkw5xJ2XKZhu/cP/blT1TtF5BeBG4BHi8gngV9S1Qcg1DU4MvJze5vfF7e85X0jPx9k+/Ga6b4xDMMw8ocFu4ZhGLsAEXkiXhDzObyJ/t2qekWIl+qWdh6Gl816JvB5VR2IyFcBidCdFlAd+f0CPIn1JO4D9otIaZuAV7d7gc8DwKUiUhgJ7PYDYWXVG4Q8/3F9CcNzgdE9qIJXXOx84Lmq6sb0XnGMS5T3v4/wnztU9YPAB0VkGfhj4LeAn4zpM7iVS0d+3o83NjP13zAMw8g/JmM2DMM4gxGRZRF5HnAj8Geq+g3gX4EVvwBQRUSKIvIYPyDeyhHgESO/L+AFOMf89n8aeEzEbn0VeKn/vs/BqzQchn8FHgTeLCILIjIvIt+/Qz9H+SKwDrxWRBzxvHefjzcmUQlz/uP6goiU/AJIRaDon0fJf+7hwNwWmfIfAo/C2xvbjvJeE4hjXKK8f+jPnXj+vs8QkTmgg5eNHfhPx/EZ3MqrReQSEdmDt7f3Q7P03zAMwzgzsGDXMAzjzORjIrKKl436v4DfBX4aQFUHeEHN4/Aq+x4H/gSv2M9W/gfwBvGq9f6Kqt4G/A7webxA57HAP0fs23/x37+BVxjor8K8aKTfl+MVRjqEV2X6tH5ueV0PeAFwHd65vgN4WYh9r9v1Icz579gXnzfgBW+vA/4P/+egINIPMSJh9rOY/xnvWh0Wr+rxmoj8RMj3GncucYxL6PeP+Lmbw7PMOo4nid6HF4SGvQZR+SBeoa67/H9vnLH/hmEYxhmAqM6qxjIMwzAMIwwicjNesaSbJx5sxIKI3AP8jKr+XdZ9MQzDMNLFMruGcQYhIj8hIp9K4X2eJiJh9lgahhGNzwKfyboThmEYhnE2YMGuYeQQEfkBEfkXEWmKyEMi8s8i8kRV/YCqPivr/hmGMR2q+pZt9uUahpEhO33nJvh+94jIDybVvmEYm1g1ZsPIGX5l0r8G/k/gw0AZeArQzbJfhmEYZyKqeiDrPhj5JY/fuTtUpDcMYwoss2sY+eNKAFX9c1UdqGpbVT+lql8XkZeLyOeCA0XkWSJyh78a/Q4R+QcR+Rn/uZeLyOdE5LdF5KSI3C0i14289qdF5Fsisioid4nIf07/VA3DMAwjUyZ95/6ziPxP/3v2dhF5ZvBCEamJyLtF5EERuV9E3igixZHnXznyPXubiDxeRN6PZ3/1Mb8Y3WtF5ICIqIi8QkTuBT693Xai0YywiNwgIn8hIn/mt/8NEblSRF4vIkdF5D4RMSWYcdZjwa5h5I9vAwMReZ+IXCci52x3kIjsBT4CvB44F7gD+L4th32P//he4C3Au0Uk8Ko8CjwPWMar4vtWEXl83CdjGIZhGDlm0nfu9+BV8N4L/Abwl76FFcD7gD5eBflrgGcBwYLzi4AbgJfhfc++ADihqj+JV23++aq6qKpvGXmvp+LZkD07ZN+fD7wfOAf4CvBJvLn9xcBv4vlXG8ZZjQW7hpEzVHUF+AE8n8l3AcdE5CYROX/Loc8FblXVv/TlTm/Ds/AY5aCqvsu31HgfcCFwvv8+H1fV76jHP+DZcjwluTMzDMMwjHwR4jv3KPB7quqq6ofwFpB/yH/+OuAXVXVdVY8CbwWu91/3M8BbVPVL/vfsnap6cEJ3bvDbCruv/59U9ZP+HOAvgPOAN6uqi+enfUBE6iHbMoxdiQW7hpFDVPVbqvpyVb0EeAxwEfB7Ww67CM9jNXiN4vmSjnJ45PmW/+MigL+C/QW/GEcDL3jeG++ZGIZhGEa+mfCde7+e6tN50H/+YYADPOj7UDfwMqn7/OMuBb4TsSv3TT7kFI6M/NwGjvuL28Hv4H/nG8bZigW7hpFzVPV24L14X8CjPAhcEvziy5MvIQQiMgf8L+C3gfNVtQ7cDMjYFxqGYRjGLmab79yLR7b/gLff9gG8wLQL7FXVuv9vWVUf7R93H3DZTm8T4vF1oBr84u8FPi/KuRiGYcGuYeQOEblKRH5ZRC7xf78UeAnwhS2Hfhx4rIj8sIiUgFcDF4R8mzIwBxwD+n7hKitkYRiGYZxVhPjO3Qf8gog4/j7cRwE3q+qDeNt/fkdElkWkICKXichT/df9CfArIvIE8bhcRB7mP3cEeMSErn0bmBeRHxIRB3gD3ve2YRgRsGDXMPLHKl5BjC+KyDreF+43gV8ePUhVjwMvwis8dQK4GriFEHYJqroK/AKezcJJ4KXATfGdgmEYhmGcEUz6zv0icAVwHHgT8OO
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAEjCAYAAADzMBbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eXxkSXXn+z25KBetmapFJamUpd5pYAztZrUZYxpMg83iBU8Dzzb+YPPswdt4wTCPZzM2jIFnDx7GZjxgbHgsBsz4eQA3Bg+LbbbG3W6g6W666W5VZkklVakyU1vuS7w/4qZKpZJKmcpIKTPrfD+f/Ei6eW/cc0MhZZw4vzhHjDEoiqIoiqIoiqIoSj/hO2wDFEVRFEVRFEVRFMU16uwqiqIoiqIoiqIofYc6u4qiKIqiKIqiKErfoc6uoiiKoiiKoiiK0neos6soiqIoiqIoiqL0HersKoqiKIqiKIqiKH2HOruKoiiKoiiKoihK36HOrqIoiqIoiqIoitJ3qLOrKIrSx4jIaRF5rqO2bhSRe0VkXUR+xUWbLdz7fhF59kHes5c47N+zy/sfNv30LIqiKFc76uwqiqJ0Md7Eu+A5Hisi8hUR+QURafn/t4NJ/OuALxpjho0x72yjnSuyk53GmMcbY77YqXu2Yss+2/klEblbREoi8r5dzpkUkfltx64XkaKIfNC1TVfgir9ndQYVRVGUXkGdXUVRlO7nRcaYYSABvBX4beC9h2BHArj/EO7bD5wF3gz8xRXOeSHw99uO/SnwL50yahf096woiqL0BersKoqi9AjGmFVjzCeAfwf8jIg8ATYjgv9TRJZFZG4n6amIfACYAT4pIhsi8jrv+OtF5FEvcvyAiPzoTvcWkc8DPwj8iXf9Dd5xIyLXbTnvfSLyZu/70yLymyLyLRFZFZGPikh4y7knReRvPLvTIvInV7BzM5ooIo8TkS96ke77ReTF22y94n23nbvr8+9my34wxvyNMeZvgfQVTnshcOeW+98BrACfa9KmJzXzzF47O/bhbr/ndu7fzPjc0v5vi8iC9/t4SERu2/LelX5Xp0Xkt7z750TkvSJyXEQ+7Z3/v0Uktu38N3jtZEXkL68wRpq2X1EUReku1NlVFEXpMYwxXwfmgWeJlTN/EvgmMAXcBvyaiDx/2zU/BaSwUeIhY8zbvbceBZ4FjAL/CfigiJzY4Z7PAf4Z+CXv+oebNPcngduBWeDfAK8CEBE/8CkgCZzybP/IFezEuy7oPe9ngWPALwMfEpEbm7nvDuz6/HvZ4hLvuf4t8A/ezyPA7wG/sfW8PWxq6pmv1Id7/Z5bvX+z49Oz60bgl4CneEqG5wOnt5yy11j9ceB5wA3Ai4BPA/8ROIKd72x3Ul/p3eNa75o37mBT0/YriqIo3Yc6u4qiKL3JWSAOPAU4aoz5PWNM2RjzGPAe4I5mGjHG/LUx5qwxpm6M+SjwXeCpDu18p9d+Bus0PMk7/lRgEvgtY0zOGFM0xnypifaeDgwBb/We9/NYp/nlTd73Eg7g+Zvl3wLfNMasez//PvBeY8yZFtpo6plpvg9bZaf7tzI+a0AIuFlEgsaY08aYRxtvNvG7+m/GmHPGmAWsw36XMeZeY0wJ+P+AJ2+7358YY8549r5ll+dv6+9LURRFOVwCh22AoiiKsi+mgAx2f+WkiKxsec+PnezviYj8NPDr2OgqWCfoiDszWdryfR7r4AKcBJLGmGqL7U0CZ4wx9S3Hktj+aOa+l9Du84vIF4Ef2OXtLxtjvr/JpjYlzCLyJOC5XO6c7UVTz0zzfdgqO92/6fFpjHlERH4NeBPweBH5DPDrxpiz0NTv6tyW7ws7/Dy07ZZbFxKS7Nxfbf19KYqiKIeLOruKoig9hog8BeuYfAk7gZ8zxlzfxKVmWzsJbJTqNuCrxpiaiHwDkBbMyQPRLT9PYCXWe3EGmBGRwA4Or9npAo+zwEkR8W1x1maAZmXVmzT5/FeyBWPMs1u97y68EGjsQX021qFLiQjY37FfRG42xtyyl01N0G4ftnL/MzQ/PjHGfBj4sCfj/h/A24CfcjRWt3Nyy/cz2H5py35FURSlu1AZs6IoSo8gIiMi8iPAR4APGmPuA74OrHmJfSIi4heRJ3gO8XbOAdds+XkQ67gse+3/LPCEFs36BvAK7763s3uUcztfBxaBt4rIoIiEReT7drFzK3cBOeB1IhIUW3v3Rdg+aZVmnv9KtjSNiAS8BEh+rOMaFpGA994sEDLGfMc7/d3YfaRP8l5/Bvwddn+pC5va7cNW7t/0+BRb3/c5IhICithobM1728VY3c5rRWRaROLYvb0fbcd+RVEUpftQZ1dRFKX7+aSIrGOjTP8X8F+AnwUwxtSwjsqTgDngAvDn2CQ+2/kD4I1iM/D+pjHmAeCPgK9iHZgnAl9u0bZf9e6/gk3487fNXLTF7uuwCY/msVmmL7Nz23Vl4MXAC7DP+i7gp7c4ik3T5PPvakuLvBHrvL0e+D+87xsJkX6YLVmYjTF5Y8xS4wVsAEVjzLILmxz0YdP3b3F8hrCltS5gJdHHsE5os7+rVvkwNknXY97rzW3aryiKonQZYky7aihFURRFUfaLiNyJTZZ0554nK04QkdPAzxlj/vdh26IoiqJ0Do3sKkoPISKvFJHPHsB9ni0izey7VBSlfb4IfOGwjVAURVGUfkOdXUXpQkTk+0XkKyKyKiIZEfmyiDzFGPMhY8wPHbZ9iqK4wxjzdmNM4bDtUJSrld0+czt4v9Mi8txOta8oykU0G7OidBleFtJPAb8IfAwYAJ4FlA7TLkVRlH7BGHPqsG1QuoNu/MzdJUu9oij7QCO7itJ93ABgjPkrY0zNGFMwxnzWGPMtEXmViHypcaKI/JCIPOStRr9LRP5RRH7Oe+9VIvIlEflDEcmKyJyIvGDLtT8rIg+KyLqIPCYi/+fBP6qiKIqiHCp7feZ+WUT+m/c5+x0Rua1xoYiMish7RWRRRBZE5M0i4t/y/s9v+Zx9QERuEZEPYEtdfVJENkTkdSJySkSMiLxaRFLA53faTrQ1IiwibxKRvxaRD3rt3yciN4jIG0TkvIicERFVgilXPersKkr38TBQE5H3i8gLRCS200kicgT4OPAGYBx4CHjmttOe5h0/ArwdeK+INOpSngd+BBjBZvZ9h4jc4vphFEVRFKWL2esz92nYbN1HgN8F/sYrVwXwfqCKzSr/ZOCHgMaC88uANwE/jf2cfTGQNsb8FDYD/YuMMUPGmLdvudcPAI/jYpmxvXgR8AEgBtwLfAY7t58Cfg9bq1pRrmrU2VWULsMYswZ8P7am5HuAZRH5hIgc33bqC4H7jTF/48md3okt17GVpDHmPV75jPcDJ4Dj3n3+zhjzqLH8I7YEx7M692SKoiiK0l008Zl7HvhjY0zFGPNR7ALyD3vvvwD4NWNMzhhzHngHcId33c8BbzfG/Iv3OfuIMSa5hzlv8tpqdg//PxtjPuPNAf4aOAq81RhTwdbNPiUiY022pSh9iTq7itKFGGMeNMa8yhgzDTwBmAT+eNtpk9i6q41rDLZW6VaWtryf974dAvBWsL/mJeNYwTrPR9w+iaIoiqJ0N3t85i6YS+t0Jr33E0AQWPRqTq9gI6nHvPNOAo+2aMqZvU+5hHNbvi8AF7zF7cbP4H3mK8rVijq7itLlGGO+A7wP+wG8lUVguvGDJ0+epglEJAT8T+APgePGmDHgTkCueKGiKIqi9DE7fOZObdn+A3a/7VmsY1oCjhhjxrzXiDHm8d55Z4Brd7tNE8dzQLTxg7cX+Ggrz6Ioijq7itJ1iMhNIvIbIjLt/XwSeDnwtW2n/h3wRBF5qYgEgNcCE03eZgAIActA1UtcpYksFEVRlKuKJj5zjwG/IiJBbx/u44A7jTGL2O0/fyQiIyLiE5FrReQHvOv+HPhNEflesVwnIgnvvXPANXuY9jAQFpEfFpEg8Ebs57aiKC2gzq6idB/r2IQYd4lIDvuB+23gN7aeZIy5ALwMm3gqDdwM3E0T5RKMMevAr2DLLGSBVwCfcPc
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "for n in [1, -1, 4, 2, -4]:\n",
    "    index_t0 = min(n_samples//n, n_samples-1)\n",
    "\n",
    "    t_0 = t[index_t0]\n",
    "    signal_step = np.zeros(n_samples)\n",
    "    signal_step[index_t0] = 1\n",
    "\n",
    "    fig, axes = plot_signal_and_spectrum(signal_step, sample_rate, ft_kwargs={'mask_bias':True}, spectrum_kwargs={'plot_complex':True, 'plot_amplitude':True});\n",
    "    fig.suptitle(\"Delta function at ${}/{}$th of the sample\".format(np.sign(n), np.abs(n)))\n",
    "    axes.flat[1].set_xlim(0, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FTs: gaussian with sine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAEjCAYAAAAc8n3PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxb53ng+9+DHeAOLuIuiZRkybJl2ZYX2U7s2E5iN16zOk3apLeZtGl729w2t9POzbRpPjM3mcw0TdrJpE2b3mTqJq7j2I5T27ETJ3biVYutxdZikRTFTVxBcANAbO/94wAUSJHiBhAk9Xw/H35IHByc8xxJ1MGD532fV4wxKKWUUkoppZRS64Ut3wEopZRSSimllFLZpImuUkoppZRSSql1RRNdpZRSSimllFLriia6SimllFJKKaXWFU10lVJKKaWUUkqtK5roKqWUUkoppZRaVzTRVUoppZZJRJ4WkU/kO45MIrJJRIyIOPIdi1JKKbXSNNFVSim15ojIAyLymohMiEh/6uffExHJRzzGmDuNMd/Nx7mVUkopdT5NdJVSSq0pIvInwNeB/w5UAxuA3wVuBFx5DE0ppZRSq4QmukoppdYMESkBvgj8njHmEWPMmLG8YYz5mDFmMrXf+0TkDREZFZFOEflCxjFuEZGuGcdtF5HbUz9fKyIHUq/tE5GvprZ7RORBERkSkaCI7BeRDannnheRT6V+bhaRn6f2GxSRfxWR0hnn+pyIHBGRERH5NxHxzHG9nxSRl0Tk71L7nhCR22aLO/X4CyLy4AWO1SYiYyJyWkQ+lvHc/yEix0VkWESeEZGNC/5LUUoppVYhTXSVUkqtJXsBN/CjefabAH4TKAXeB3xGRO5b4Dm+DnzdGFMMNAMPp7Z/AigBGoByrCpyeJbXC/AloBbYkdr/CzP2+TBwB7AZ2AV88gLxXAe0ARXAXwKPioh/gddiBSRSAPwtcKcxpgi4ATiUeu4+4D8B7wcqgV8B31/M8ZVSSqnVRhNdpZRSa0kFMGiMiac3iMjLqQprWETeCWCMed4Yc9QYkzTGHMFK3G5e4DliwBYRqTDGjBtjXs3YXg5sMcYkjDEHjTGjM19sjGkxxvzUGDNpjBkAvjrLuf/WGNNjjAkAPwZ2XyCefuBrxpiYMebfgJNYyftiJYHLRMRrjDlrjHkrtf13gC8ZY46n/lz/X2C3VnWVUkqtZZroKqWUWkuGgIrMTsLGmBuMMaWp52wAInKdiPxCRAZEZASr+lqxwHP8NrANOJEannxXavu/AM8AD4lIj4h8RUScM18sIlUi8pCIdIvIKPDgLOfuzfg5BBReIJ5uY4zJeHwGq1q8YMaYCeAjWH8OZ0XkSRHZnnp6I/D11IcFQSCAVZWuW8w5lFJKqdVEE12llFJrySvAJHDvPPt9D3gCaDDGlAB/j5W8gTWs2ZfeUUTsWEN2ATDGnDLGfBSoAv4b8IiIFKQqqn9ljLkUa+jvXVjDo2f6EmCAXanhzx/POPdS1M3oJt0I9Mx2LVjNuWZljHnGGPNuoAY4Afxj6qlO4HeMMaUZX15jzMvLiFkppZTKK010lVJKrRnGmCDwV8D/EpEPikihiNhEZDdQkLFrERAwxkRE5Frg1zOeexvwpBpWOYHPY837BUBEPi4ilcaYJBBMbU6IyLtE5PJUYjyKNZQ5MUuYRcA4EBSROuD/XuZlVwF/KCJOEfkQ1rzfp1LPHQIeSD23B/jgbAcQkQ0ick9qru5kKr507H8P/LmI7EztW5I6j1JKKbVmaaKrlFJqTTHGfAX4Y+BPseav9gH/APxHIF2F/D3giyIyBvwF5xpKYYwZST3/T0A3VlU0swvzHcBbIjKO1ZjqAWNMBKta+ghWknsceAFrWPJMfwVcBYwATwKPLvOSXwO2AoPAfwU+aIwZSj33n7EaZg2nzvu9OY5hA/4EqxIcwJoz/HsAxpjHsCrXD6WGWr8J3LnMmJVSSqm8kunTfpRSSim1WojIJ4FPGWNuyncsSiml1FqiFV2llFJKKaWUUuuKJrpKKaWUUkoppdYVTXSVWiNE5GMi8uwKnOcWEemaf0+lVK4ZY76jw5aVUkqpxdNEV6lVRkRuEpGXRWRERAIi8pKIXGOM+VdjzHvyHZ9SSim1Xsx1z83h+dpF5PZcHV8pdY4j3wEopc4RkWLg34HPYHWJdQHvwFoORCmllFJZshrvuSLiMMbE83V+pdYTregqtbpsAzDGfN8YkzDGhI0xzxpjjojIJ0XkxfSOIvIeETmZ+hT6f4nICyLyqdRznxSRF0Xkf4jIsIicFpE7M177WyJyXETGRKRNRH5n5S9VKaWUyqv57rkvicjfpe6zJ0TktvQLU+tNf1tEzopIt4j8l9Qa2+nn/0PGffaYiFwlIv8CNAI/FpFxEflTEdkkIkZEfltEOoCfzzaFKLMSLCJfEJEfiMiDqeMfFZFtIvLnItIvIp0ioiPA1EVPE12lVpe3gYSIfFdE7hSRstl2EpEKrPU8/xwoB04CN8zY7brU9grgK8C3RURSz/UDdwHFwG8BfyMiV2X7YpRSSqlVbL577nVAG9Z99C+BR0XEn3ruu0Ac2AJcCbwHSH/Y/CHgC8BvYt1n7wGGjDG/AXQAdxtjClNrgqfdDOwA3rvA2O8G/gUoA94AnsF6X18HfBFrbXGlLmqa6Cq1ihhjRoGbAAP8IzAgIk+IyIYZu/4a8JYx5tHUEKe/BXpn7HPGGPOPxpgE1g25BtiQOs+TxphWY3kBeBZruJZSSil1UVjAPbcf+JoxJmaM+TesD4/fl3r+TuCzxpgJY0w/8DfAA6nXfQr4ijFmf+o+22KMOTNPOF9IHSu8wPB/ZYx5JvUe4AdAJfBlY0wMeAjYJCKlCzyWUuuSJrpKrTLGmOPGmE8aY+qBy4Ba4GszdqsFOjNeY4CZnZJ7M54PpX4sBEh9cv1qqvFGECtxrsjulSillFKr2zz33O7U/TXtTOr5jYATOCsiwdR99B+AqtR+DUDrIkPpnH+Xafoyfg4Dg6kPttOPIXXPV+pipYmuUquYMeYE8B2sm2+ms0B9+kFqSHI9CyAibuCHwP8ANhhjSoGnALngC5VSSql1bJZ7bl3GlB+w5tf2YCWlk0CFMaY09VVsjNmZ2q8TaJ7rNAvYPgH40g9Sc38rF3MtSilNdJVaVURku4j8iYjUpx43AB8FXp2x65PA5SJyn4g4gN8Hqhd4GhfgBgaAeKpJlTatUEopdVFZwD23CvhDEXGm5t3uAJ4yxpzFmvLz1yJSLCI2EWkWkZtTr/sn4HMicrVYtojIxtRzfUDTPKG9DXhE5H0i4gQ+j3XfVkotgia6Sq0uY1jNL14TkQmsm+2bwJ9k7mSMGQQ+hNVkagi4FDjAApZEMMaMAX+ItZTCMPDrwBPZuwSllFJqTZjvnvsasBUYBP4r8EFjzFDqud/E+uD4GNa99BGsXhgYY36Q2v97qXM8DqSbWH0J+HxqyPPnZgvKGDMC/B5WwtyNVeGdOT1JKTUPmT71QCm1FomIDesm+DFjzC/yHY9SSim1lonIJ4FPGWNuyncsSqml0YquUmuUiLxXREpTc27/E9Yc25lDnJVSSimllLroaKKr1Nq1F6ur4yDWenr3LWJZAqWUygoR+XS+Y8iW9XIt6+U6QK9ltVov17JergP0WmY9jg5dVkoppdRSicgBY8yefMeRDevlWtbLdYBey2q1Xq5lvVwH6LXMRiu6SimllFJKKaXWlTVV0a2oqDCbNm3KdxhKKaXWiYMHDw4aY3R9ymUoLCw027dvz3cYWTEwMEBl5dr/57BergP0Wlar9XIt6+U6YH1dy8GDB8eNMUXLPY4jG8EsR2oR7ANAtzHmrgvtu2nTJg4cOLAygSmllFr3RORMvmNY67Zv3673ZqWUUlkjIiezcZzVMHT5j4Dj+Q5CKaWUUkoppdT6kNdEV0TqgfdhLYitlFJKKaWUUkotW74rul8D/hRIzrWDiHxaRA6IyIGBgYGVi0wppZRSSiml1JqUt0RXRO4C+o0xBy+0nzHmW8aYPcaYPetlgrV
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = 1e1 # Hz\n",
    "if f > sample_rate/2:\n",
    "    print(\"Sampling a frequency above the sample_rate/2 gives problems\")\n",
    "\n",
    "def gaussian(x, mu=0, sigma=1 ):\n",
    "    return 1/(sigma*np.sqrt(2*np.pi)) * np.exp(- (x-mu)**2 / sigma**2 )\n",
    "\n",
    "def dgaussian(x, mu=0, sigma=1, dmu=1, dsigma=0):\n",
    "    return gaussian(x, mu, sigma)   + gaussian(x, mu+dmu, sigma+dsigma)\n",
    "\n",
    "if not True:\n",
    "    signal = gaussian(t, t[2*len(t)//5], 100/len(t)) - gaussian(t, t[3*len(t)//5], 100/len(t))\n",
    "else:\n",
    "    signal = gaussian(t, t[2*len(t)//5], 100/len(t)) * np.sin(2*np.pi*f*t)\n",
    "\n",
    "fig, axes = plot_signal_and_spectrum(signal, sample_rate, \"Gaussian pulse\", spectrum_kwargs={'plot_complex':True, 'plot_amplitude':True})\n",
    "axes.flat[1].axvline(sample_rate/2, color='r', label='Nyquist Frequency')\n",
    "axes.flat[2].axvline(sample_rate/2, color='r')\n",
    "axes.flat[1].legend();\n",
    "\n",
    "if not True:\n",
    "    axes.flat[1].set_xlim(0, 20);\n",
    "else:\n",
    "    axes.flat[1].set_xlim(max(0, f - max(f/5, 10)), f + max(f/5, 10));\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FTs: delta peak: group velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above plots, the phases follow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Peculiarities $F^{-1}( F(A(t)) ) \\sim A(t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA74AAAEWCAYAAABBrkurAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxcdb3/8ddnJpkmk3RJmlCWUlqQHUqBsF0qFBEoXED2RVCLYpWr1OvVq6JeQeVelUUERaFgf0WtZV8FEVQKgiCUrZS2LLK1FGya7tln5vP745xJp22STptkJid5Px+PNGf5njmfyUznzOd8N3N3RERERERERAaqWLEDEBEREREREelLSnxFRERERERkQFPiKyIiIiIiIgOaEl8REREREREZ0JT4ioiIiIiIyICmxFdEREREREQGNCW+IiIiItLvmdl5ZvZIAc4zycyW9PV5wnPNNLPLC3GuYjCzOWZ2YbHjEAElviKRYGYTzezvZrbazFaY2VNmdlAfnu8dM/t4Xz2+iIhIZ7q73rn7LHc/ttgxFoqZTTGzJ4sdR18xs8vMrN3M1uX8fCPcN8fMWjbad+RG625mjTnrHy32c5L+raTYAYhI98xsGPAH4CLgdiABfBRoLWJMJe6eKtb5RURk4OmP17vBxswMMHfPFOiUt7n7+V3s+7K737zRtsrsgpk5sJ+7v9ln0cmAohpfkf5vNwB3n+3uaXdvdvdH3H1eeDf4KTP7eXh3fJGZHZ090MyGm9mvzewDM3vfzC43s3jO/s+b2UIzW2tmC8zsADP7LTAGeCB799XMxoZ3Vj9nZu8Bf+2sKVhuTXF4J/cOM/td+PivmNluZnaJmS0zs8VmNmju3IuIyGZ1eb2DTWtAzexYM3stvP790swezzarzZY1s6vMbKWZvW1mx+cce0HO9e8tM/tCvkGa2bXhNWyNmT2fW9MYXvtuN7PfhI/9qpnV5ezf38xeCPfdBpRtwXnfMbOvm9m88DnfZmZl4b6FZnZiTtkSM1tuZgeE64eGNemrzOxlM5uUU3aOmf2vmT0FNAE7h3+/t8I43zaz83LKfzY830oz+5OZ7ZSz75jwu8hqM/sFYPk+P5G+psRXpP97HUib2S1mdryZVW20/xDgLaAGuBS428yqw323ACngI8D+wLFA9kvBmcBlwKeBYcDJQIO7fwp4DzjJ3Svd/Yqccx0J7Akcl2fsJwG/BaqAF4E/EXzu7AD8ALgxz8cREZGBb3PXuw5mVgPcCVwCjAReA/5to2KHhNtrgCuAX4c1mgDLgBMJrn8XANdkk8Q8PAdMAKqB3wN3ZBPQ0MnArcAI4H7gF2HMCeBegutiNXAHcHqe58w6C5gMjAPGA1PC7bOBc3PKHQcsd/cXzGwH4EHg8vC8XwfuMrPanPKfAqYCQ4F64DrgeHcfSvB3fSl8DqcA3wZOA2qBv4Xnzr4mdwHfJfib/xM4fAufn0ifUeIr0s+5+xpgIuDATUC9md1vZqPCIsuAn7l7u7vfRnCR//dw//HAf7p7o7svA64BzgmPuxC4wt2f88Cb7v7uZsK5LHys5jzD/5u7/ylsFn0HwUXyx+7eTvClYKyZjcjzsUREZADL43qX6wTgVXe/O7zGXAd8uFGZd939JndPE9wI3g4YFZ7rQXf/Z3j9exx4hKBZdT5x/s7dG9w95e5XA0OA3XOKPOnuD4Xn/S2wX7j9UKCU9dfsOwmS6C1xnbsvdfcVwAMECTgECfjJZpYM1z8ZbgM4H3gojCnj7o8Ccwn+hlkz3f3V8G+ZAjLAPmZW7u4fuPurYbkvAD9y94Vh2f8DJoS1vicAC9z9zvA6/zM2fU02dlZYC5392T73ueZsf2HL/kwim1LiKxIB4QVmiruPBvYBtie4oAC87+6eU/zdcP9OBBfYD7IXDoIa1m3CcjsS3I3dEou3sPy/cpabCe4+p3PWIae/joiIDG6bud7l2p6ca1J4Hdx4JOYPc/Y3hYuVAGGN8jMWDKC1iiBpq8knRjP7WtjUd3V47PCNjs1N9pqAMjMrCWPu7Jq9JTZ+7EqAsJ/rQuCkMPk9mfWJ707AmbkJJsENhu1yHiv3b9kInA18keA7xINmtkfOY12b8zgrCJoz70Dnr8nmvjfc7u4jcn6W5uyblrM939p4kS4p8RWJGHdfBMwk+EIAsENO0y0I+ucuJbjYtAI1OReOYe6+d1huMbBLV6fJY3sjkL2zjAV9h2s3OUJERGQrdHK9y/UBMDq7El4HR3dSbhNmNoSgSe5VwCh3HwE8RB79UcP+vN8kaHJcFR67Op9jw5g7u2b3lmxz508Q1LxmB31aDPx2owSzwt1/nHPsBtf9sLXWMQTJ8SKCGvjsY31ho8cqd/e/h89vx+xjhM9zR0T6CSW+Iv2cme0R3l0eHa7vSHBheyYssg0wzcxKw367exI0afqAoOnW1WY2zMxiZraLmR0ZHncz8HUzO9ACH8kZoOJfwM6bCe11grvY/25mpQR9eob01vMWEZHBJY/rXa4HgX3N7JSwNvVLwLZ5nipBcL2qB1IWDHqV72CLQwmaAtcDJWb2PYJ+wvl4Ojx2Wjj41GnAwXkem49bCZ7HRayv7QX4HUFN8HFmFjezMgsGqOz0RoGZjTKzk82sguAG+jog21rrBuASM9s7LDs8/O4BwWuyt5mdFr4m08j/NRHpc0p8Rfq/tQQDdPzDzBoJvgDMB74W7v8HsCuwHPhf4Ax3bwj3fZrgAr8AWEkwEMh2AO5+R1j+9+E57iUY9ALgR8B3w6ZMX+8sKHdfDfwHQQL9PkEN8MbNzERERPK1uetdB3dfDpxJMGhVA7AXQb/VzU595O5rCZKy2wmujZ8kGIQqH38C/khw8/ddoIU8uwG5exvBoFBTwvOeDdyd53nzefwPCJLrfwNuy9m+mKAW+NsECfti4L/pOg+IEfzNlxI0ZT6S4HqPu98D/AS41czWELw+x4f7sq/Jjwlek12Bp3rr+Yn0lG3YzUBEosTMpgAXuvvEYsciIiJSLGYWI7j5ep67P1bseESk/1GNr4iIiIhETth0d0TYZ/fbBP1sO2sWLSKixFdEREREIukwgtkJlhPMG3/KFky3JyKDjJo6i4iIiIiIyICmGl8REREREREZ0EqKHUAh1dTU+NixY4sdhoiIDBDPP//8cnfX/NU9oGuziIj0pq6uzYMq8R07dixz584tdhgiIjJAmNm7xY4h6nRtFhGR3tTVtVlNnUVERERERGRAU+IrIiIiIiIiA5oSXxERERERERnQBlUfXxER6R3t7e0sWbKElpaWYodSEGVlZYwePZrS0tJihzIoDLb3V3+m976IDBRKfEVEZIstWbKEoUOHMnbsWMys2OH0KXenoaGBJUuWMG7cuGKHMygMpvdXf6b3vogMJEVt6mxmM8xsmZnN72L/f5vZS+HPfDNLm1l1uO8dM3sl3KfhIEVECqilpYWRI0cOiqTEzBg5cqRqHwtoML2/+jO990VkICl2H9+ZwOSudrr7le4+wd0nAJcAj7v7ipwiR4X76/o4zg08/PoLfPvRX/Pw6y8U8rQiIv3KYEpKBtNz7S/0N+8f9DqIyEBR1KbO7v6EmY3Ns/i5wOy+iyY/D7/+At/421dwa+fB95LAVUze7YBihyUiIiJSMEtXNbPowzV8bI9RxQ5FRCQvxa7xzYuZJQlqhu/K2ezAI2b2vJlN7ebYqWY218zm1tfX9ziWJ959EbcW3NpxUjzx7os9fkwREdly8XicCRMmsM8++3DSSSexatWqrX6ssWPHsnz58l6MTqKusrKy2CF0a9KkScydu2lPr0mTJrH77rszYcIEJkyYwJ133gms//+S/bnxxhs7lisrKzuO+fSnP53X+U+5/ik+O1M9zUQkOqIyuNVJwFMbNXM+3N2Xmtk2wKNmtsjdn9j4QHefDkwHqKur854GcsRO+/PA4lKgCSPOETvt39OHFBGRrVBeXs5LL70EwGc+8xmuv/56vvOd7xQ5KpH8pNNp4vF4nzz2rFmzqKvbsBdY7v+XrC984QtAkCxfddVVmxzTnWVrWwF
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_0 = t[n_samples//4]\n",
    "signal_step = np.heaviside(t - t_0, 1) + 1\n",
    "\n",
    "signal_fft = ft.rfft(signal_step) / sample_rate\n",
    "signal_ifft = ft.irfft(signal_fft, len(signal_fft))\n",
    "ft_freqs = ft.rfftfreq(n_samples, 1/sample_rate)\n",
    "\n",
    "\n",
    "alpha = 0.8\n",
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,4))\n",
    "ax1.set_title(\"Spectrum\")\n",
    "ax1.set_xlabel(\"f (Hz)\")\n",
    "ax1.set_ylabel(\"Amplitude\")\n",
    "ax1.set_xlim(-1,20)\n",
    "ax1.plot(ft_freqs, np.real(signal_fft), '.-', label='Real', alpha=alpha)\n",
    "ax1.plot(ft_freqs, np.imag(signal_fft), '.--', label='Imag', alpha=alpha)\n",
    "ax1.plot(ft_freqs, np.abs(signal_fft),  '.-', label='Abs', alpha=alpha)\n",
    "ax1.legend()\n",
    "\n",
    "# Signal and Inversed FFT\n",
    "ax2.set_title(\"Signal and Inversed FFT\")\n",
    "ax2.set_xlabel(\"t (s)\")\n",
    "ax2.set_ylabel(\"A(t)\")\n",
    "ax2.plot(np.arange(len(signal_ifft)) / sample_rate * 2, np.abs(signal_ifft)*sample_rate/2, label=\"Inversed FFT\")\n",
    "ax2.plot(np.arange(len(signal_step)) / sample_rate, signal_step, label=\"Original\")\n",
    "ax2.legend()\n",
    "\n",
    "if False:\n",
    "    ax2.set_xlim(t_0*(1-0.2), t_0*(1+0.2))\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DFT\n",
    "\n",
    "A naive implementation of a DFT is simple to setup by implementing Eq. (1)\n",
    "\n",
    "$$\n",
    "\\hat{f}(\\omega) = F(t) = \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N-1} f(t_k)\\, \\exp{ \\left(-i \\omega t_k \\right)}\n",
    ".\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def naive_DFT( signal ):\n",
    "    \"\"\"\n",
    "    Calculate DFT in a naive way, returns a complex DFT.\n",
    "    \"\"\"\n",
    "\n",
    "    N = np.size(signal)\n",
    "\n",
    "    x_k = 2*np.pi/N*np.arange(N)\n",
    "    ft = np.zeros(N, dtype=np.complex128)\n",
    "    for j in range(N):\n",
    "        ft[j] = np.dot(signal, np.exp(- 1j * j * x_k ))\n",
    "\n",
    "    return ft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DFT: Cooley and Tukey (or Divide and Conquer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A smarter DFT (Cooley and Tukey) can be implemented by splitting the sum in Eq. (1) into two parts: $k$ even and $k$ odd.\n",
    "\n",
    "With a uniform sample space $t$, the label $t_k$ can be rewritten to $t_k = k \\Delta t = $.\n",
    "\n",
    "This gives\n",
    "$$\n",
    "F(t) = \n",
    "    \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k})\\, \\exp{ \\left(-i \\omega \\, t_{2k} \\right)}\n",
    "    +\n",
    "    \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k+1})\\, \\exp{ \\left(-i \\omega \\, t_{2k+1} \\right)}\n",
    "    \\\\\n",
    "    =\n",
    "    \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k})\\, \\exp{ \\left(-i \\omega \\, 2k \\Delta t \\right)}\n",
    "    +\n",
    "    \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k+1})\\, \\exp{ \\left(-i \\omega \\, 2k \\Delta t \\right)}\\exp{ \\left(-i \\omega \\, \\Delta t \\right)}\n",
    "    \\\\\n",
    "    =\n",
    "    \\frac{1}{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k})\\, \\exp{ \\left(-i \\omega \\, 2k \\Delta t \\right)}\n",
    "    +\n",
    "    \\frac{\\exp{ \\left(-i \\omega \\, \\Delta t \\right)} }{\\sqrt{2\\pi}} \\sum_{k=0}^{N/2-1} f(t_{2k+1})\\, \\exp{ \\left(-i \\omega \\, 2k \\Delta t \\right)}\n",
    "    \\\\\n",
    "    =\n",
    "    F(t_{even}) + F(t_{odd}) \\exp{ \\left(-i \\omega \\, \\Delta t \\right)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cooley_tukey_DFT( signal ):\n",
    "    \"\"\"\n",
    "    Calculate a DFT according to the Cooley and Tukey algorithm (radix-2)\n",
    "    also known as Divide and Conquer.\n",
    "    \n",
    "    It divides the terms into 2 streams with odd and even index.\n",
    "    \"\"\"\n",
    "    radix = 2\n",
    "\n",
    "    N = np.size(signal)\n",
    "    \n",
    "    if N <= radix:\n",
    "        return naive_DFT(signal)\n",
    "    \n",
    "    if N % radix > 0:\n",
    "        raise ValueError(\"Provide a signal with  N%{radix} == 0 values.\".format(radix=radix))\n",
    "    \n",
    "    # divide and conquer\n",
    "    DFT_even = cooley_tukey_DFT(signal[::radix])\n",
    "    DFT_odd  = cooley_tukey_DFT(signal[1::radix])\n",
    "    \n",
    "    phase_terms = np.exp( -2j * np.pi * np.arange(0, N) / N )\n",
    "    \n",
    "    return np.tile(DFT_even, 2) + np.tile(DFT_odd,2) * phase_terms\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA74AAAEjCAYAAAAVGd21AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOx9ebxfVXXvd90ARhEkTCLzIIMICoqIT1qsdcB56gBaFavlaUux1dZin1VL9TnUPq0tRZmcBQQRARmkxaCgAcIcQgIZSHIJkJEkDCHJPev9cc75/fZee619zrk3yZ3Wl0+495y999prrb3PuXud7x6ImeFwOBwOh8PhcDgcDsdExcBoK+BwOBwOh8PhcDgcDseWhAe+DofD4XA4HA6Hw+GY0PDA1+FwOBwOh8PhcDgcExoe+DocDofD4XA4HA6HY0LDA1+Hw+FwOBwOh8PhcExoeODrcDgcDofD4XA4HI4JDQ98HQ6Hw+FwOBwOh8MxoeGBr8PhcDgcDofD4XA4JjQ88HU4HA5HI4joPiJ6zWaS9RARvW4Y5Q4lojuJaB0Rnb45dJlsIKLdiOh6IlpNROePtj4Oh8PhcGwteODrcDgcDgAAER1PRL8lojVEtIqIbiaiVwAAM7+YmaePsoqfAjCdmXdg5m+Osi5jDi0/KHwawIPMPI2ZPzyCunYmop8R0ZNEtIiI3jvc/E2yWqSfRkQziegZIvquSHuIiDYQ0a7i/l1ExES0v8j7OpHvFCK6qaVbHA6HwzGG4YGvw+FwOEBEOwK4CsB/ANgZwF4A/hnAM6Opl8B+AO7TEohom62sy3jF6wBcshnknAVgA4DnA3gfgLOJ6MXDzN8kqyl9KYAvALjAqHshgJPrCyI6EsCzW9jocDgcjgkED3wdDofDAQCHAAAzX8jMQ8z8NDP/kpnvAWI2rPr974jonoodvpiIptaCiOhlwZTkS6r0L2iVEtGeRPRTIlpORAutKcxEdAOAPwDwn0T0BBEdUunxD0R0D4AniWibJnlEdDQR3VHpdjERXVTrVjGALwzyfjfUOye7hU/2IaLLqrIrieg/q/t/T0Q/FTr+BxF9Q/HBGUQ0v9J9NhG9K0j7AYB9AVxZ+edToux2RLQGwJFVnns1P7cBEW0P4D0A/omZn2DmmwBcAeD9XfM3yWpTFzNfxsyXA1hpqPwDAB8Irj8I4PvDsPtPK9/W/54hould5TgcDodjdOCBr8PhcDgA4AEAQ0T0PSJ6ExFNa8j/JwBOBHAAgJcAOAUoAywAPwPwXZTM8YUA3qUJIKIBAFcCuBslw/yHAP6GiN4o8zLzawH8BsBpzPxcZn6gSjoZwFsA7ASgyMmrdLscZSC0M0rm8z0NdnbR1fLJFJRs+iIA+1flL6rK/BDAiUS0U5V3GwB/WukoMR/A7wF4Hko2/odE9ILKP+8HsBjA2yr/fFX4bwOAVwFYVqUfKey7iogeN/5dJfQ4BMBQ0Aao/GIxvrn8TbK61qVhBoAdiehFVVv8KUq/dwIzX1z57rkA9gSwAGX/djgcDsc4gAe+DofD4QAzrwVwPAAGcC6A5UR0BRE93yjyTWZeysyrUAaER1X3jwOwTZW+kZkvA3CrIeMVAHZj5jOZeQMzL6jqPqmD6t9k5iXM/HQLeccB2BbANyrdLgVwW8t62uhq+eRYlIHS3zPzk8y8vmIuwcyPAPg1gD+u8p4IYAUz3y4VYOZLKvkFM18M4MFKdlschTJoTMDMb2XmnYx/bxXZnwtgjbi3BsAORr25/E2yutZloWZ9Xw9gDoCHjXyXh0E/gP+SGaqPID9Gud782x31cDgcDscowddEORwOhwMAwMz3o89SHoaSFfsGgvWRAR4Nfn8KZWCH6ufDzMxB+hKjyv0A7FkFGDWmoGR22yKU3SRP021Ry3ra6Gr5ZB8Ai5h5kyH7ewA+hjKQ/jPobC+I6AMAPoGSNQbKoHBXLa8BM/DtiCcA7Cju7Qhg3TDyN8nqWpeFH6D8wHAA8tOc38nM/11fENEpAD4i8nwRZeDtO4s7HA7HOIIzvg6Hw+FIwMxzUE5XPqJj0UcA7EVEFNzbx8i7BMBCwS7uwMxv7qJqB3mabvsGvz8F4DnB9R6bSdclAPYlewOuywG8hIiOAPBWAD+SGYhoP5SB8WkAdmHmnQDMAhDawrKcwEthBL5EdI1Yvxr+u0ZkfwDANkR0sJCtbjzWkL9JVte6VDDzIpSbXL0ZwGVdyoYgopNQfgj6I2beOFw5DofD4dj68MDX4XA4HCCiw4jok0S0d3W9D8oB/oyOon4HYAjAadVmU++APR33VgBrqw2qnk1EU4joCKqOUBoGmuT9DsAmAKdXur1b6HYXgPdW5U4EcMJm0vVWlEH3l4loeyKaSkSvrhOZeT2AS1FOn72VmRcrMrZHGdguBwAi+hDSjxKPATgwo4cZ+DLzm+r1q8q/N4m8T6IMHs+s7Hk1gHfAYKpz+ZtktamrasupKBn4KZV/tY8MHwbw2kpmZxDR0Sh3PX8nMy8fjgyHw+FwjB488HU4HA4HUE4dfSWAW4joSZQB7ywAn+wipNpE6d0og4zHUU7dvQrKsUjMPATgbSin4C4EsALAeSg3b+qMJnmBbqcAWI1yk6OQ/ft4Vf5xlMfmXL45dA3KvhDlBlSDVd0hvodyx2UreJwN4N9QBu+PVXlvFtm+BOAz1frUvwsTiGgPANNQrm/dHPhLlEcCLUO5wdPHmLnHwlYM8j+2zJ+V1SL9MwCeBnAGyv72dHUvAjPPZ+aZw7a4DLinAbgpw4Y7HA6HY4yC4qVODofD4XBsXhDRLQC+xczfGW1dJIjouwAGmTkJlLayHvuiDEr3qDYaczgcDofDsRnhjK/DMUlARO8jol9uhXpeQ0SDW7oex9gFEZ1ARHtUU1A/iPJon2tHW6+ximqX4E8AuMiDXofD4XA4tgw88HU4JhiI6Hgi+i0RrSGiVUR0MxG9gpl/xMxvGG39HJMCh6JcS7oG5VTpP6qO7XEIENH2ANaiPGbnc6OsjsPh6ADr7+0WrO8hInrdlpLvcEx0+HFGDscEAhHtiHI95ccA/ATAdgB+D8r6SodjS4GZzwFwzmjr0QbMfMoo1/8kymOJHA7HOMJY/HtLRNtkjk1zOCY9nPF1OCYWDgEAZr6QmYeY+Wlm/iUz30NEpxDRTXVGInoDEc2tvlT/FxHdSEQfqdJOIaKbiOhrRLSaiBYS0ZuCsh8iovuJaB0RLSCi/731TXU4HA6HY9TQ9Pf2ZiL6j+pv7Bwi+sO6IBE9j4jOJ6JHiOhhIvoCEU0J0v8i+Bs7m4heRkQ/QHn82pXVxmqfIqL9iYiJ6MNEtBjADdpyo5ApJqLPE9ElRPTDSv69RHQIEX2aiJYR0RIi8tlhjgkJD3wdjomFBwAMEdH3iOhNRDRNy0REu6I8PuXTAHYBMBfA/xLZXlnd3xXAVwGcT9Q7/3QZyvNGdwTwIQBfJ6KXbW5jHA6Hw+EYo2j6e/tKAAtQ/g39HIDLiGjnKu17KI9WeyGAowG8AUD94fmPAXwewAdQ/o19O4CVzPx+lLvCv606ZuyrQV0nAHgRgDe21P1tKHeQnwbgTgDXoYwJ9gJwJoBvt5TjcIwreODrcEwgVBvjHI/yvM9zASwnoiuI6Pki65sB3MfMl1XTor4J4FGRZxEzn1sdxfI9AC8A8Pyqnl9UR4MwM98I4Jcop3g5HA6HwzHh0eLv7TIA32Dmjcx8McoPyW+p0t8E4G+Y+UlmXgbg6wBOqsp9BMBXmfm26m/sPGZe1KDO5ytZT7dU/zfMfF319/8SALsB+DIzbwRwEYD9iWinlrIcjnEDD3wdjgkGZr6fmU9h5r0BHAFgTwDfENn2BLAkKMMozxYN8WiQ/lT163MBoPq6PaPazONxlIH0rpvXEofD4XA4xi4a/t4+zPGZoYuq9P0AbAvgkerM7cdRMqy7V/n2ATC/oypLmrNEeCz4/WkAK6qP3PU14HsPOCYgPPB1OCYwmHkOgO+i/IMc4hEAe9cX1RTmvdECRPQsAD8F8DUAz2fmnQBcDYCyBR0Oh8P
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dft_sample_rate = 1/1e-3 # s\n",
    "dft_n_samples = 2**10\n",
    "dft_t = np.arange(dft_n_samples) / dft_sample_rate\n",
    "\n",
    "f = 1e2 # Hz\n",
    "if f > dft_sample_rate/2:\n",
    "    print(\"Sampling a frequency above the sample_rate/2 gives problems\")\n",
    "\n",
    "signal_sine = np.sin( 2*np.pi*f*dft_t )\n",
    "\n",
    "fig, axes = plot_signal_and_spectrum(signal_sine, dft_sample_rate, \"Single frequency at $f = {:g}$MHz\".format(f/10**6), ft_kwargs={'fft': cooley_tukey_DFT})\n",
    "axes.flat[1].axvline(dft_sample_rate/2, color='r', label='Nyquist Frequency');"
   ]
  }
 ],
 "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
}