PCP_09_dft.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<a href=\"https://www.audiolabs-erlangen.de/fau/professor/mueller\"><img src=\"data_layout/PCP_Teaser.png\" width=100% style=\"float: right;\" alt=\"PCP Teaser\"></a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unit 9: Discrete Fourier Transform (DFT)\n",
    "\n",
    "<ul> \n",
    "    <li><a href='#learn'>Overview and Learning Objectives</a></li>    \n",
    "    <li><a href='#inner'>Inner Product</a></li>\n",
    "    <li><a href='#dft'>Definition of DFT</a></li>\n",
    "    <li><a href='#phase'>DFT Phase</a></li>\n",
    "    <li><a href='#dftmatrix'>DFT Matrix</a></li>\n",
    "    <li><a href='#fft'>Fast Fourier Transform (FFT)</a></li>\n",
    "    <li><a href='#exercise_freq_index'>Exercise 1: Interpretation of Frequency Indices</a></li>\n",
    "    <li><a href='#exercise_missing_time'>Exercise 2: Missing Time Localization</a></li>    \n",
    "    <li><a href='#exercise_chirp'>Exercise 3: Chirp Signal</a></li>   \n",
    "    <li><a href='#exercise_inverse'>Exercise 4: Inverse DFT</a></li>    \n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='learn'></a> \n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<h2>Overview and Learning Objectives</h2>\n",
    "\n",
    "    \n",
    "The <strong>Fourier transform</strong> is one of the most important tools for a wide range of engineering and computer science applications. The general idea of <strong>Fourier analysis</strong> is to decompose a given signal into a weighted superposition of sinusoidal functions. Since these functions possess an explicit physical meaning regarding their frequencies, the decomposition is typically more accessible for subsequent processing steps than the original signal. Assuming that you are familiar with the Fourier transform and its applications in signal processing, we review in this unit the discrete variant of the Fourier transform known as <strong>Discrete Fourier Transform</strong> (DFT). We define the inner product that allows for comparing two vectors (e.g., discrete-time signals of finite length). The DFT can be thought of as comparing a given signal of finite length with a specific set of exponential signals (a complex variant of sinusoidal signals), each comparison yielding a complex-valued Fourier coefficient. Then, using suitable visualizations, we show how you can interpret the amplitudes and phases of these coefficients. Recall that one can express the DFT as a complex-valued square matrix. We show how separately plotting the real and imaginary parts leads to beautiful and insightful images. Applying a DFT boils down to computing a matrix&ndash;vector product, which we implement via the standard NumPy function <code>np.dot</code>. Since the number of operations for computing a DFT via a simple matrix&ndash;vector product is quadratic in the input length, the runtime of this approach becomes problematic with increasing length. This issue is exactly where the fast Fourier transform (FFT) comes into the game. We present this famous divide-and-conquer algorithm and provide a Python implementation. Furthermore, we compare the runtime behavior between the FFT implementation and the naive DFT implementation. We will further deepen your understanding of the Fourier transform by considering further examples and visualization in the exercises. In <a href='#exercise_freq_index'>Exercise 1</a>, you will learn how to interpret and plot frequency indices in a physically meaningful way. In <a href='#exercise_missing_time'>Exercise 2</a>, we discuss the issue of loosing time information when applying the Fourier transform, which is the main motivation for the <a href='https://www.audiolabs-erlangen.de/resources/MIR/FMP/C2/C2_STFT-Basic.html'>short-time Fourier transform</a>. In <a href='#exercise_chirp'>Exercise 3</a>, you will apply the DFT to a <strong>chirp signal</strong>, which yields another illustrative example of the DFT's properties. Finally, in <a href='#exercise_inverse'>Exercise 4</a>, we will invite you to explore the relationship between the DFT and its inverse. Again, an overarching goal of this unit is to apply and deepen your Python programming skills within the context of a central topic for signal processing. \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='inner'></a> \n",
    "## Inner Product\n",
    "\n",
    "In this notebook, we consider [discrete-time (DT) signals](PCP_08_signal.html) of finite length $N\\in\\mathbb{N}$, which we represent as vector \n",
    "\n",
    "$$\n",
    "x=(x(0),x(1),...,x(N-1))^\\top\\in\\mathbb{R}^N\n",
    "$$ \n",
    "\n",
    "with samples $x(n)\\in\\mathbb{R}^N$ for $n\\in[0:N-1]$. Note that $\\top$ indicates the transpose of a vector, thus converting a row vector into a column vector. Furthermore, note that we start indexing with the index $0$ (thus adapting our mathematical notation to Python conventions). A general concept for comparing two vectors (or signals) is the **inner product**. Given two vectors $x, y \\in \\mathbb{R}^N$, the inner product between $x$ and $y$ is defined as follows:\n",
    "\n",
    "$$ \n",
    "\\langle x | y \\rangle := \\sum_{n=0}^{N-1} x(n) y(n).\n",
    "$$\n",
    "\n",
    "The absolute value of the inner product may be interpreted as a measure of similarity between $x$ and $y$. If $x$ and $y$  are similar (i.e., if they point to more or less the same direction), the inner product $|\\langle x | y \\rangle|$ is large. If $x$ and $y$ are dissimilar (i.e., if $x$ and $y$ are more or less orthogonal to each other), the inner product $|\\langle x | y \\rangle|$ is close to zero.\n",
    "\n",
    "One can extend this concept to **complex-valued** vectors $x,y\\in\\mathrm{C}^N$, where the inner product is defined as \n",
    "\n",
    "$$\n",
    "\\langle x | y \\rangle := \\sum_{n=0}^{N-1} x(n) \\overline{y(n)}.\n",
    "$$\n",
    "\n",
    "In the case of real-valued signals, the complex conjugate does not play any role and the definition of the complex-valued inner product reduces to the real-valued one. In the following code cell, we give some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<strong>Note:</strong>\n",
    "One can use the NumPy function <code>np.vdot</code> to compute the inner product. However, opposed to the mathematical convention that conjugates the second argument, this function applies complex conjugation on the first argument. Therefore, for computing $\\langle x | y \\rangle$ as defined above, one has to call <code>np.vdot(y, x)</code>.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following, we generate and visualize three signals $x_1$, $x_2$, $x_3$. Then, we compute and discuss different inner products using the signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import libpcp.signal\n",
    "%matplotlib inline\n",
    "\n",
    "Fs = 64\n",
    "dur = 1\n",
    "x1, t = libpcp.signal.generate_example_signal(Fs=Fs, dur=dur)\n",
    "x2, t = libpcp.signal.generate_sinusoid(dur=dur, Fs=Fs, amp=1, freq=2, phase=0.3)\n",
    "x3, t = libpcp.signal.generate_sinusoid(dur=dur, Fs=Fs, amp=1, freq=6, phase=0.1)\n",
    "\n",
    "def plot_inner_product(ax, t, x, y, color_x='k', color_y='r', label_x='x', label_y='y'):\n",
    "    \"\"\"Plot inner product\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        ax: Axis handle\n",
    "        t: Time axis\n",
    "        x: Signal x\n",
    "        y: Signal y\n",
    "        color_x: Color of signal x (Default value = 'k')\n",
    "        color_y: Color of signal y (Default value = 'r')\n",
    "        label_x: Label of signal x (Default value = 'x')\n",
    "        label_y: Label of signal y (Default value = 'y')\n",
    "    \"\"\"\n",
    "    ax.plot(t, x, color=color_x, linewidth=1.0, linestyle='-', label=label_x)\n",
    "    ax.plot(t, y, color=color_y, linewidth=1.0, linestyle='-', label=label_y)\n",
    "    ax.set_xlim([0, t[-1]])\n",
    "    ax.set_ylim([-1.5, 1.5])\n",
    "    ax.set_xlabel('Time (seconds)')\n",
    "    ax.set_ylabel('Amplitude')\n",
    "    sim = np.vdot(y, x)\n",
    "    ax.set_title(r'$\\langle$ %s $|$ %s $\\rangle = %.1f$' % (label_x, label_y, sim))\n",
    "    ax.legend(loc='upper right')    \n",
    "\n",
    "plt.figure(figsize=(8, 5))\n",
    "ax = plt.subplot(2, 2, 1)\n",
    "plot_inner_product(ax, t, x1, x1, color_x='k', color_y='k', label_x='$x_1$', label_y='$x_1$')\n",
    "ax = plt.subplot(2, 2, 2)\n",
    "plot_inner_product(ax, t, x1, x2, color_x='k', color_y='r', label_x='$x_1$', label_y='$x_2$')\n",
    "ax = plt.subplot(2, 2, 3)\n",
    "plot_inner_product(ax, t, x1, x3, color_x='k', color_y='b', label_x='$x_1$', label_y='$x_3$')\n",
    "ax = plt.subplot(2, 2, 4)\n",
    "plot_inner_product(ax, t, x2, x3, color_x='r', color_y='b', label_x='$x_2$', label_y='$x_3$')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example, one can make the following observations:\n",
    "\n",
    "* The signal $x_1$ is similar to itself, leading to a large value of $\\langle x_1 | x_1 \\rangle=40.0$.\n",
    "* The overall course of the signal $x_1$ strongly correlates with the sinusoid $x_2$, which is reflected by a relatively large value of $\\langle x_1 | x_2 \\rangle=29.9$.\n",
    "* There are some finer oscillations of $x_1$ that are captured by $x_3$, leading to a still noticeable value of  $\\langle x_1 | x_3 \\rangle=14.7$. \n",
    "* The two sinusoids $x_2$ and $x_3$ are more or less uncorrelated, which is revealed by the value of $\\langle x_2 | x_3 \\rangle\\approx 0$. \n",
    "\n",
    "In other words, the above comparison reveals that the signal $x_1$ has a strong signal component of $2~\\mathrm {Hz}$ (frequency of $x_2$) and $6~\\mathrm {Hz}$ (frequency of $x_3$). Measuring correlations between an arbitrary signal and sinusoids of different frequencies is exactly the idea of performing a Fourier (or spectral) analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='dft'></a> \n",
    "## Definition of DFT\n",
    "\n",
    "Let $x\\in \\mathbb{C}^N$ be a vector of length $N\\in\\mathbb{N}$. The **discrete Fourier transform** (DFT) of $x$ is defined by:\n",
    "\n",
    "$$ X(k) := \\sum_{n=0}^{N-1} x(n) \\exp(-2 \\pi i k n / N) $$\n",
    "\n",
    "for $k \\in [0:N-1]$. The vector $X\\in\\mathbb{C}^N$ can be interpreted as a frequency representation of the time-domain signal $x$.   To obtain a geometric interpretation of the DFT, we define the vector $\\mathbf{e}_k \\in\\mathbb{C}^N$ with real part $\\mathbf{c}_k=\\mathrm{Re}(\\mathbf{e}_k)$ and imaginary part $\\mathbf{s}_k=\\mathrm{Im}(\\mathbf{e}_k)$ by\n",
    "\n",
    "$$\\mathbf{e}_k(n) :=  \\exp(2 \\pi i k n / N) = \\cos(2 \\pi i k n / N) + i \\sin(2 \\pi i k n / N)\n",
    "= \\mathbf{c}_k(n) + i \\mathbf{s}_k(n)$$\n",
    "\n",
    "for each $k \\in [0:N-1]$.\n",
    "\n",
    "\n",
    "This vector can be regarded as a [sampled version](PCP_08_signal.html) of the [exponential function](PCP_07_exp.html) of frequency $k/N$. Using inner products, the DFT can be expressed as\n",
    "\n",
    "$$ X(k) = \\sum_{n=0}^{N-1} x(n) \\overline{\\mathbf{e}_k}(n) = \\langle x | \\mathbf{e}_k \\rangle,$$\n",
    "\n",
    "thus measuring the similarity between the signal $x$ and the sampled exponential functions $\\mathbf{e}_k$. The absolute value $|X(k)|$ indicates the degree of similarity between the signal $x$ and $\\mathbf{e}_k$. In the case that $x\\in \\mathbb{R}^N$ is a real-valued vector (which is typically the case for audio signals), we obtain:\n",
    "\n",
    "$$ \n",
    "X(k) = \\langle x |\\mathrm{Re}(\\mathbf{e}_k) \\rangle - i\\langle x | \\mathrm{Im}(\\mathbf{e}_k) \\rangle\n",
    "= \\langle x |\\mathbf{c}_k \\rangle - i\\langle x | \\mathbf{s}_k \\rangle\n",
    "$$\n",
    "\n",
    "The following plot shows an example signal $x$ compared with functions $\\overline{\\mathbf{e}_k}$ for various frequency parameters $k$. The real and imaginary part of $\\overline{\\mathbf{e}_k}$ are shown in <font color='red'> red</font> and <font color='blue'> blue</font>, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_signal_e_k(ax, x, k, show_e=True, show_opt=False):\n",
    "    \"\"\"Plot signal and k-th DFT sinusoid\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        ax: Axis handle\n",
    "        x: Signal\n",
    "        k: Index of DFT\n",
    "        show_e: Shows cosine and sine (Default value = True)\n",
    "        show_opt: Shows cosine with optimal phase (Default value = False)\n",
    "    \"\"\"\n",
    "    N = len(x)\n",
    "    time_index = np.arange(N)\n",
    "    ax.plot(time_index, x, 'k', marker='.', markersize='10', linewidth=2.0, label='$x$')\n",
    "    plt.xlabel('Time (samples)')\n",
    "    e_k = np.exp(2 * np.pi * 1j * k * time_index / N)\n",
    "    c_k = np.real(e_k)\n",
    "    s_k = np.imag(e_k)\n",
    "    X_k = np.vdot(e_k, x)\n",
    "\n",
    "    plt.title(r'k = %d: Re($X(k)$) = %0.2f, Im($X(k)$) = %0.2f, $|X(k)|$=%0.2f' %\n",
    "              (k, X_k.real, X_k.imag, np.abs(X_k)))\n",
    "    if show_e is True:\n",
    "        ax.plot(time_index, c_k, 'r', marker='.', markersize='5',\n",
    "                 linewidth=1.0, linestyle=':', label='$\\mathrm{Re}(\\overline{\\mathbf{u}}_k)$')\n",
    "        ax.plot(time_index, s_k, 'b', marker='.', markersize='5',\n",
    "                 linewidth=1.0, linestyle=':', label='$\\mathrm{Im}(\\overline{\\mathbf{u}}_k)$')\n",
    "    if show_opt is True:\n",
    "        phase_k = - np.angle(X_k) / (2 * np.pi)\n",
    "        cos_k_opt = np.cos(2 * np.pi * (k * time_index / N - phase_k))\n",
    "        d_k = np.sum(x * cos_k_opt)\n",
    "        ax.plot(time_index, cos_k_opt, 'g', marker='.', markersize='5',\n",
    "                 linewidth=1.0, linestyle=':', label='$\\cos_{k, opt}$')\n",
    "    plt.grid()\n",
    "    plt.legend(loc='lower right')\n",
    "\n",
    "\n",
    "N = 64\n",
    "x, t = libpcp.signal.generate_example_signal(Fs=N, dur=1)\n",
    "\n",
    "plt.figure(figsize=(8, 15))\n",
    "for k in range(1, 8):\n",
    "    ax = plt.subplot(7, 1, k)\n",
    "    plot_signal_e_k(ax, x, k=k)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='phase'></a>\n",
    "## DFT Phase\n",
    "\n",
    "At first sight, the DFT may be a bit confusing: Why is a real-valued signal $x$ compared with a complex-valued sinusoid $\\mathbf{e}_k$? What does the resulting complex-valued Fourier coefficient\n",
    "\n",
    "$$\n",
    "c_k:= X(k) := \\langle x |\\mathrm{Re}(\\mathbf{e}_k) \\rangle - i\\langle x | \\mathrm{Im}(\\mathbf{e}_k) \\rangle. \n",
    "$$\n",
    "\n",
    "encode? To understand this, we represent the complex number $c_k$ in form of its [polar representation](PCP_06_complex.html#polar)\n",
    "\n",
    "$$\n",
    "c_k = |c_k| \\cdot \\mathrm{exp}(i \\gamma_k),    \n",
    "$$\n",
    "\n",
    "where $\\gamma_k$ is the [angle](PCP_06_complex.html) (given in radians).  Furthermore, let $\\mathbf{cos}_{k,\\varphi}:[0:N-1]\\to\\mathbb{R}$ be a sampled sinusoid with frequency parameter $k$ and phase $\\varphi\\in[0,1)$, defined by\n",
    "\n",
    "$$\n",
    "   \\mathbf{cos}_{k,\\varphi}(n) = \\mathrm{cos}\\big( 2\\pi (kn/N - \\varphi) \\big)\n",
    "$$\n",
    "\n",
    "for $n\\in[0,N-1]$. Defining $\\varphi_k := - \\frac{\\gamma_k}{2 \\pi}$, one obtains the following remarkable property of the Fourier coefficient $c_k$: \n",
    "\n",
    "\\begin{eqnarray}\n",
    "|c_k| &=& \\mathrm{max}_{\\varphi\\in[0,1)} \\langle x | \\mathbf{cos}_{k,\\varphi} \\rangle,\\\\\n",
    "\\varphi_k  &=& \\mathrm{argmax}_{\\varphi\\in[0,1)} \\langle x | \\mathbf{cos}_{k,\\varphi} \\rangle.\n",
    "\\end{eqnarray}\n",
    "\n",
    "In other words, the phase $\\varphi_k$ maximizes the correlation between $x$ and all possible sinusoids $\\mathbf{cos}_{k,\\varphi}$ with $\\varphi\\in[0,1)$. Furthermore, the magnitude $|c_k|$ yields this maximal value. Thus, computing a single correlation between $x$ and the complex-valued function $\\mathbf{e}_k$ (which real part coincides with $\\mathbf{cos}_{k,0}$, and its imaginary part with $\\mathbf{cos}_{k,0.25}$) solves an optimization problem. In the following code cell, we demonstrate this optimality property, where the $\\mathbf{cos}_{k,\\varphi}$ with optimal phase $\\varphi=\\varphi_k$ is shown in <font color='green'>green</font>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 15))\n",
    "for k in range(1, 8):\n",
    "    ax = plt.subplot(7, 1, k)\n",
    "    plot_signal_e_k(ax, x, k=k, show_e=False, show_opt=True)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='dftmatrix'></a>\n",
    "## DFT Matrix\n",
    "\n",
    "Being a linear operator $\\mathbb{C}^N \\to \\mathbb{C}^N$, the DFT can be expressed by some $N\\times N$-matrix. This leads to the famous DFT matrix $\\mathrm{DFT}_N \\in \\mathbb{C}^{N\\times N}$ matrix, which is given by\n",
    "\n",
    "$$\\mathrm{DFT}_N(n, k) = \\mathrm{exp}(-2 \\pi i k n / N)$$\n",
    "\n",
    "for $n\\in[0:N-1]$ and $k\\in[0:N-1]$. Let $\\rho_N:=\\exp(2 \\pi i / N)$ be the primitive $N^\\mathrm{th}$ [root of unity](PCP_07_exp.html#roots). Then \n",
    "\n",
    "$$\\sigma_N:= \\overline{\\rho_N} = \\mathrm{exp}(-2 \\pi i / N)$$\n",
    "\n",
    "also defines a primitive $N^\\mathrm{th}$ [root of unity](PCP_07_exp.html#roots). From the [properties of exponential functions](PCP_07_exp.html), one obtains that\n",
    "\n",
    "$$ \\sigma_N^{kn} = \\mathrm{exp}(-2 \\pi i / N)^{kn} = \\mathrm{exp}(-2 \\pi i k n / N)$$\n",
    "\n",
    "From this, one obtains:\n",
    "\n",
    "$$\n",
    "\\mathrm{DFT}_N =\n",
    "\\begin{pmatrix}\n",
    "    1 & 1 & 1 & \\dots  & 1 \\\\\n",
    "    1 & \\sigma_N & \\sigma_N^2 & \\dots  & \\sigma_N^{N-1} \\\\\n",
    "    1 & \\sigma_N^2 & \\sigma_N^4 & \\dots  & \\sigma_N^{2(N-1)} \\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    1 & \\sigma_N^{N-1} & \\sigma_N^{2(N-1)} & \\dots  & \\sigma_N^{(N-1)(N-1)} \\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "In the following visualization, the real and imaginary part of $\\mathrm{DFT}_N$ are shown, where the values are encoded by suitable colors. Note that the $k^\\mathrm{th}$ row of $\\mathrm{DFT}_N$ corresponds to the vector $\\mathbf{e}_k$ as defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_matrix_dft(N, K):\n",
    "    \"\"\"Generate a DFT (discete Fourier transfrom) matrix\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        N: Number of samples\n",
    "        K: Number of frequency bins\n",
    "\n",
    "    Returns:\n",
    "        dft: The DFT matrix\n",
    "    \"\"\"\n",
    "    dft = np.zeros((K, N), dtype=np.complex128)\n",
    "    time_index = np.arange(N)\n",
    "    for k in range(K):\n",
    "        dft[k, :] = np.exp(-2j * np.pi * k * time_index / N)\n",
    "    return dft\n",
    "\n",
    "N = 32\n",
    "dft_matrix = generate_matrix_dft(N, N)\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.title('$\\mathrm{Re}(\\mathrm{DFT}_N)$')\n",
    "plt.imshow(np.real(dft_matrix), origin='lower', cmap='seismic', aspect='equal')\n",
    "plt.xlabel('Time (sample, index $n$)')\n",
    "plt.ylabel('Frequency (index $k$)')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.title('$\\mathrm{Im}(\\mathrm{DFT}_N)$')\n",
    "plt.imshow(np.imag(dft_matrix), origin='lower', cmap='seismic', aspect='equal')\n",
    "plt.xlabel('Time (samples, index $n$)')\n",
    "plt.ylabel('Frequency (index $k$)')\n",
    "plt.colorbar()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write a function that computes the discrete Fourier transform $X = \\mathrm{DFT}_N \\cdot x$ of a signal $x\\in\\mathbb{C}^N$. We apply the function from above sampled at $N=64$ time points. The peaks of the magnitude Fourier transform $|X|$ correspond to the main frequency components the signal is composed of. Note that the magnitude Fourier transform is symmetrical around the center. Why? For the interpretation of the time and frequency axis, see also <a href='#exercise_freq_index'>Exercise 1: Interpretation of Frequency Indices</a></li>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(x):\n",
    "    \"\"\"Compute the discete Fourier transfrom (DFT)\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        x: Signal to be transformed\n",
    "\n",
    "    Returns:\n",
    "        X: Fourier transform of x\n",
    "    \"\"\"\n",
    "    x = x.astype(np.complex128)\n",
    "    N = len(x)\n",
    "    dft_mat = generate_matrix_dft(N, N)\n",
    "    return np.dot(dft_mat, x)\n",
    "\n",
    "N = 64\n",
    "x, t = libpcp.signal.generate_example_signal(Fs=N, dur=1)\n",
    "X = dft(x)\n",
    "\n",
    "def plot_signal_dft(t, x, X, ax_sec=False, ax_Hz=False, freq_half=False, figsize=(10, 2)):\n",
    "    \"\"\"Plotting function for signals and its magnitude DFT\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        t: Time axis (given in seconds)\n",
    "        x: Signal\n",
    "        X: DFT\n",
    "        ax_sec: Plots time axis in seconds (Default value = False)\n",
    "        ax_Hz: Plots frequency axis in Hertz (Default value = False)\n",
    "        freq_half: Plots only low half of frequency coefficients (Default value = False)\n",
    "        figsize: Size of figure (Default value = (10, 2))\n",
    "    \"\"\"\n",
    "    N = len(x)\n",
    "    if freq_half is True:\n",
    "        K = N // 2\n",
    "        X = X[:K]\n",
    "    else:\n",
    "        K = N\n",
    "\n",
    "    plt.figure(figsize=figsize)\n",
    "    ax = plt.subplot(1, 2, 1)\n",
    "    ax.set_title('$x$ with $N=%d$' % N)\n",
    "    if ax_sec is True:\n",
    "        ax.plot(t, x, 'k', marker='.', markersize='3', linewidth=0.5)\n",
    "        ax.set_xlabel('Time (seconds)')\n",
    "    else:\n",
    "        ax.plot(x, 'k', marker='.', markersize='3', linewidth=0.5)\n",
    "        ax.set_xlabel('Time (samples)')\n",
    "    ax.grid()\n",
    "\n",
    "    ax = plt.subplot(1, 2, 2)\n",
    "    ax.set_title('$|X|$')\n",
    "    if ax_Hz is True:\n",
    "        Fs = 1 / (t[1] - t[0])\n",
    "        ax_freq = Fs * np.arange(K) / N\n",
    "        ax.plot(ax_freq, np.abs(X), 'k', marker='.', markersize='3', linewidth=0.5)\n",
    "        ax.set_xlabel('Frequency (Hz)')\n",
    "\n",
    "    else:\n",
    "        ax.plot(np.abs(X), 'k', marker='.', markersize='3', linewidth=0.5)\n",
    "        ax.set_xlabel('Frequency (index)')\n",
    "    ax.grid()\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "\n",
    "plot_signal_dft(t, x, X)\n",
    "plot_signal_dft(t, x, X, ax_sec=True, ax_Hz=True)\n",
    "plot_signal_dft(t, x, X, ax_sec=True, ax_Hz=True, freq_half=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='fft'></a>\n",
    "## Fast Fourier Transform (FFT)\n",
    "\n",
    "Next, we discuss the famous fast Fourier transform (FFT), which is a fast algorithm to compute the DFT. The FFT algorithm was originally found by Gauss in about 1805 and  then rediscovered by Cooley and Tukey in 1965. The FFT algorithm is based on the observation that applying a DFT of even size $N=2M$ can be expressed in terms of  applying two DFTs of half the size $M$. It exploits the fact that there are algebraic relations between the entries $\\sigma_N^{kn} = \\mathrm{exp}(-2 \\pi i / N)^{kn}$ of DFT matrices. In particular, one has \n",
    "\n",
    "$$\\sigma_M = \\sigma_N^2$$\n",
    "\n",
    "In the FFT algorithm, one computes the DFT of the even-indexed and the uneven-indexed entries of $x$:\n",
    "\n",
    "\\begin{align}\n",
    "(A(0), \\dots, A(N/2-1)) &= \\mathrm{DFT}_{N/2} \\cdot (x(0), x(2), x(4), \\dots, x(N-2))\\\\\n",
    "(B(0), \\dots, B(N/2-1)) &= \\mathrm{DFT}_{N/2} \\cdot (x(1), x(3), x(5), \\dots, x(N-1))\n",
    "\\end{align}\n",
    "\n",
    "With these two DFTs of size $N/2$, one can compute the full DFT of size $N$ via:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "C(k) &=& \\sigma_N^k \\cdot B(k)\\\\\n",
    "X(k) &=& A(k) + C(k)\\\\\n",
    "X(N/2 + k) &=& A(k) - C(k)\\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "for $k \\in [0: N/2 - 1]$. The numbers $\\sigma_N^k$ are also called *twiddle factors*. If $N$ is a power of two, this idea can be applied recursively until one reaches the computation of $\\mathrm{DFT}_{1}$ (the case $N=1$), which is simply multiplication by one (i.e. just returning the signal of length $N=1$). For further details, we refer to Section 2.4.3 of <a href=\"http://www.music-processing.de\">[Müller, FMP, Springer 2015])</a> (see also Table 2.1). \n",
    "\n",
    "In the following code, we provide a function `fft` that implements the FFT algorithm. We test the function  `fft` by comparing its output with the one when applying the `dft` on a test signal `x`. For the comparison of result matrices, we use the NumPy functions [`np.array_equal`](https://numpy.org/doc/stable/reference/generated/numpy.array_equal.html) and [`np.allclose`](https://numpy.org/doc/stable/reference/generated/numpy.allclose.html#numpy.allclose)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft(x):\n",
    "    \"\"\"Compute the fast Fourier transform (FFT)\n",
    "\n",
    "    Notebook: PCP_09_dft.ipynb\n",
    "\n",
    "    Args:\n",
    "        x: Signal to be transformed\n",
    "\n",
    "    Returns:\n",
    "        X: Fourier transform of x\n",
    "    \"\"\"\n",
    "    x = x.astype(np.complex128)\n",
    "    N = len(x)\n",
    "    log2N = np.log2(N)\n",
    "    assert log2N == int(log2N), 'N must be a power of two!'\n",
    "    X = np.zeros(N, dtype=np.complex128)\n",
    "\n",
    "    if N == 1:\n",
    "        return x\n",
    "    else:\n",
    "        this_range = np.arange(N)\n",
    "        A = fft(x[this_range % 2 == 0])\n",
    "        B = fft(x[this_range % 2 == 1])\n",
    "        range_twiddle_k = np.arange(N // 2)\n",
    "        sigma = np.exp(-2j * np.pi * range_twiddle_k / N)\n",
    "        C = sigma * B\n",
    "        X[:N//2] = A + C\n",
    "        X[N//2:] = A - C\n",
    "        return X\n",
    "    \n",
    "N = 64\n",
    "x, t = libpcp.signal.generate_example_signal(Fs=N, dur=1)\n",
    "\n",
    "X_via_dft = dft(x)\n",
    "X_via_fft = fft(x)\n",
    "X_via_fft_numpy = np.fft.fft(x)\n",
    "\n",
    "is_equal = np.array_equal(X_via_dft, X_via_fft)\n",
    "is_equal_tol = np.allclose(X_via_dft, X_via_fft)\n",
    "is_equal_tol_np = np.allclose(X_via_dft, X_via_fft_numpy)\n",
    "\n",
    "print('Equality test for dft(x) and fft(x) using np.array_equal:    ', is_equal)\n",
    "print('Equality test for dft(x) and fft(x) using np.allclose:       ', is_equal_tol)\n",
    "print('Equality test for dft(x) and np.fft.fft(x) using np.allclose:', is_equal_tol_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<strong>Note:</strong> The test shows that our <code>dft</code> and <code>fft</code> implementations do not yield the same result (due to numerical issues). However, the results are numerically very close, which is verified by the test using <code>np.allclose</code>.  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FFT reduces the overall number of operations from the order of $N^2$ (needed when computing the usual matrix&ndash;vector product $\\mathrm{DFT}_N \\cdot x$) to the order of $N\\log_2N$. The savings are enormous. For example, using $N=2^{10}=1024$, the FFT requires roughly $N\\log_2N=10240$ instead of $N^2=1048576$ operations in the naive approach. Using the module `timeit`, which provides a simple way to time small bits of Python code, the following code compares the running time when using the naive approach and the FFT. Furthermore, we compare the running time with the highly optimized NumPy implementation <code>np.fft.fft</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timeit\n",
    "\n",
    "rep = 3\n",
    "for N in [256, 512, 1024, 2048, 4096]:\n",
    "    time_index = np.arange(N)\n",
    "    x = np.sin(2 * np.pi *  time_index / N )\n",
    "    t_DFT = 1000 * timeit.timeit(lambda: dft(x), number=rep)/rep\n",
    "    t_FFT = timeit.timeit(lambda: fft(x), number=rep*5)/(rep*5)\n",
    "    t_FFT_np = timeit.timeit(lambda: np.fft.fft(x), number=rep*100)/(rep*100)\n",
    "    print(f'Runtime (ms) for N = {N:4d} : DFT {t_DFT:10.2f},  FFT {t_FFT:.5f},  FFT_np {t_FFT_np:.8f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises and Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libpcp.dft\n",
    "show_result = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise_freq_index'></a>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<strong>Exercise 1: Interpretation of Frequency Indices</strong><br>\n",
    "Given a dimension $N\\in\\mathbb{N}$, the $\\mathrm{DFT}_N$ transform a vector $x\\in\\mathbb{C}^N$ into another vector $X\\in\\mathbb{C}^N$. Assuming that $x$ represents a time-domain signal sampled with a sampling rate $F_\\mathrm{s}$, one can associate the index $n\\in[0:N-1]$ of the sample $x(n)$ with the physical time point $t = n/F_\\mathrm{s}$ given in seconds. In case of the vector $X$, the index $k\\in[0:N-1]$ of the coefficient $X(k)$ can be associated to a physical frequency value \n",
    "\n",
    "$$\n",
    "  \\omega=\\frac{k \\cdot F_\\mathrm{s}}{N}.\n",
    "$$\n",
    "    \n",
    "Furthermore, using a real-valued signal $x\\in\\mathbb{R}^N$, the upper part of $X\\in\\mathbb{C}^N$ becomes redundant, and it suffices to consider the first $K$ coefficients with $K=N/2$.\n",
    "    \n",
    "<ul>\n",
    "<li>Find explanations why these properties apply.</li>\n",
    "<li>Find out how the function <code>plot_signal_dft</code> uses these properties to convert and visualize the time and frequency axes.</li>   \n",
    "<li>Using the signal <code>x, t = libpcp.signal.generate_example_signal(Fs=64, dur=2)</code>, plot the signal and its magnitude Fourier transform once using axes given in indices and once using axes given in physical units (seconds, Hertz). Discuss the results.</li>\n",
    "<li>Do the same for the signal <code>x, t = libpcp.signal.generate_example_signal(Fs=32, dur=2)</code>. What is going wrong and why?</li>  \n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#<solution>\n",
    "# Your Solution\n",
    "#</solution>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libpcp.dft.exercise_freq_index(show_result=show_result) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise_missing_time'></a>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<strong>Exercise 2: Missing Time Localization</strong><br>\n",
    "The Fourier transform yields frequency information that is averaged over the entire time axis. However, the information on when these frequencies occur is hidden in the transform. To demonstrate this phenomenon, construct the following two different signals defined on a common time axis $[0, T]$ with $T$ given in seconds (e.g., $T=6~\\mathrm{sec}$). \n",
    "\n",
    "<ul>\n",
    "<li>A superposition of two sinusoids $f_1+f_2$ defined over the entire time interval $[0, T]$, where the first sinusoid $f_1$ has a frequency $\\omega_1=1~\\mathrm{Hz}$ and an amplitude of $1$, while the second sinusoid $f_2$ has a frequency $\\omega_2=5~\\mathrm{Hz}$ and an amplitude of $0.5$.</li>\n",
    "<li>A concatenation of two sinusoids, where $f_1$ (specified as before) is now defined only on the subinterval $[0, T/2]$, and $f_2$ is defined on the subinterval $[T/2, T]$.\n",
    "</ul>\n",
    "    \n",
    "Sample the interval $[0,T]$ to obtain $N$ samples (use <code>np.linspace</code>), with $N\\in\\mathbb{N}$ being power of two (e.g., $N=256$). Define DT-signals of the superposition and the concatenation and compute the DFT for each of the signals. Plot the signals as well as the resulting magnitude Fourier transforms and discuss the result.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#<solution>\n",
    "# Your Solution\n",
    "#</solution>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libpcp.dft.exercise_missing_time(show_result=show_result) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise_chirp'></a>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<strong>Exercise 3: Chirp Signal</strong><br>\n",
    "The function $f(t)=\\sin\\left(\\pi t^2\\right)$ defines a <strong>chirp signal</strong> (also called <strong>sweep signal</strong>), in which the frequency increases with time. The <strong>instantaneous frequency $\\omega_t$</strong> of the chirp signal at time $t$ is the derivate of the sinusoid's argument divided by $2\\pi$, thus $\\omega_t = t$. \n",
    "<ul>\n",
    "<li>Let $[t_0,t_1]$ be a time interval (given in seconds) with $0\\leq t_0<t_1$ and $N\\in\\mathbb{N}$ be power of two. Implement a function <code>generate_chirp</code> that outputs a sampled chirp signal <code>x</code> over the interval $[t_0,t_1]$ with $N$ samples (use <code>np.linspace</code>).</li>\n",
    "<li>Compute the DFT of <code>x</code> for various input parameters $t_0$, $t_1$, and $N$. Plot the chirp signal as well as the resulting magnitude Fourier transform. Discuss the result.</li>\n",
    "</ul>    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#<solution>\n",
    "# Your Solution\n",
    "#</solution>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libpcp.dft.exercise_chirp(show_result=show_result) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise_inverse'></a>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<strong>Exercise 4: Inverse DFT</strong><br>\n",
    "The discrete Fourier transform given by the matrix $\\mathrm{DFT}_N \\in \\mathbb{C}^{N\\times N}$ is an invertible operation, given by the inverse DFT matrix $\\mathrm{DFT}_N^{-1}$. \n",
    "<ul>\n",
    "<li>There is an explicit relation between $\\mathrm{DFT}_N$ and its inverse $\\mathrm{DFT}_N^{-1}$. Which one? </li>\n",
    "<li>Write a function <code>generate_matrix_dft_inv</code> that explicitly generates $\\mathrm{DFT}_N^{-1}$. \n",
    "<li>Check your function by computing $\\mathrm{DFT}_N \\cdot \\mathrm{DFT}_N^{-1}$ and $\\mathrm{DFT}_N^{-1} \\cdot \\mathrm{DFT}_N$ (using <code>np.matmul</code>) and comparing these products with the identity matrix (using <code>np.eye</code> and <code>np.allclose</code>).</li>\n",
    "<li>Furthermore, compute the inverse DFT by using <code>np.linalg.inv</code>. Compare the result with your function using <code>np.allclose</code>.\n",
    "<li>Similar to <code>fft</code>, implement a fast inverse Fourier transform <code>fft_inv</code></li>    \n",
    "</ul>    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#<solution>\n",
    "# Your Solution\n",
    "#</solution>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libpcp.dft.exercise_inverse(show_result=show_result)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<a href=\"https://opensource.org/licenses/MIT\"><img src=\"data_layout/PCP_License.png\" width=100% style=\"float: right;\" alt=\"PCP License\"></a>\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}