03_fourier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dual-tone Multifrequency Dialing\n",
    "![](data/dtmf.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import fftpack as fft\n",
    "from sklearn.linear_model import Lasso\n",
    "from IPython.display import Audio\n",
    "plt.style.use('ggplot')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FS = 4000 # Sampling frequency in Hertz\n",
    "\n",
    "f1, f2 = 697, 1336 # Pick any two touchtone frequencies\n",
    "\n",
    "duration = 1 / 2\n",
    "n_samples = 500\n",
    "\n",
    "t = np.linspace(0, duration, int(duration * FS))\n",
    "s1 = np.sin(2 * np.pi * f1 * t)\n",
    "s2 = np.sin(2 * np.pi * f2 * t)\n",
    "S = s1 + s2\n",
    "\n",
    "fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(10, 8))\n",
    "ax[0].plot(t[:n_samples], s1[:n_samples])\n",
    "ax[0].set_title('Frequency = {} Hz'.format(f1))\n",
    "ax[1].plot(t[:n_samples], s2[:n_samples])\n",
    "ax[1].set_title('Frequency = {} Hz'.format(f2))\n",
    "ax[2].plot(t[:n_samples], S[:n_samples])\n",
    "ax[2].set_title('Touchtone signal')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(np.tile(S, 10), rate=8000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = fft.rfft(S)\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(np.abs(F))\n",
    "plt.xlabel('Frequency [Hz]')\n",
    "plt.ylabel('Spectral Power')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrete Fourier Transform\n",
    "## Convert a sequence of complex numbers* into another sequence of complex numbers*\n",
    "## $$X_{k} = \\sum_{n=0}^{N-1} x_{n}e^{-\\frac{i2\\pi}{N}kn}$$\n",
    "\n",
    "### Question: How do we convert this in equation to matrix form?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = K = S.shape[0]\n",
    "D = np.ones((N, K), dtype=complex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D *= np.exp(-1j * 2 * np.pi / N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we have is a matrix $\\mathbf{D}$, such that\n",
    "### $$\\mathbf{D} \\in \\mathbb{C}^{N \\times K}$$\n",
    "### where\n",
    "### $$d_{n, k} = e^{-i \\frac{2\\pi}{N}}$$\n",
    "### $$\\forall n \\in [0, N - 1]$$\n",
    "### $$\\forall k \\in [0, K - 1]$$\n",
    "\n",
    "## What we want is a $\\mathbf{D}$ such that\n",
    "### $$d_{n, k} = e^{-i \\frac{2\\pi}{N} \\times kn}$$\n",
    "\n",
    "## Hint: $x^{yz} = (x^{y})^{z}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmat = np.ones(D.shape, dtype=int) * np.arange(K)\n",
    "kmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nmat = np.ones(D.shape, dtype=int) * np.arange(N).reshape(-1, 1)\n",
    "nmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = D ** (kmat * nmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))\n",
    "ax[0].imshow(np.real(D))\n",
    "ax[1].imshow(np.imag(D))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moment of truth!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_HAT = np.dot(S, D)\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(np.abs(F_HAT))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question: What went wrong?\n",
    "### If $N$ is even:\n",
    "### $$\\hat{F} = [F_{0}, \\Re(F_{1}), \\Im(F_{1}), \\Re(F_{2}), \\Im(F_{2}), \\dots, \\Im(F_{N/2})]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FREAL, FIMAG = np.real(F_HAT), np.imag(F_HAT)\n",
    "F_HAT_NEW = np.c_[FREAL[:1000], FIMAG[:1000]].ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(np.abs(F_HAT_NEW))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Pick any key from the touchtone keypad, construct the signal, and find it's Spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter code here"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}