signal_processing_ee225d_f22.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d5367774",
   "metadata": {},
   "source": [
    "Basic Signal Processing review EE225D F22\n",
    "\n",
    "**Authors**: Gopala Anumanchipalli (Adapted from Anmol, Hoang & Jordan)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be645647",
   "metadata": {},
   "source": [
    "Load relevant libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fc5daee",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install -r requirements.txt\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import wavfile\n",
    "import IPython.display as ipd\n",
    "import scipy.signal as signal\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a184db48",
   "metadata": {},
   "source": [
    "Load and play audio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0304dcb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs, audio = wavfile.read(\"toms-diner.wav\")\n",
    "print(f\"Loaded {audio.size} samples at a sampling rate of {fs}Hz\")\n",
    "ipd.Audio(audio, rate=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "370e5396",
   "metadata": {},
   "source": [
    "Plot waveform and Spectrogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da27166",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.size(audio)\n",
    "plt.plot(np.arange(576000),audio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef363b66",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.specgram(audio, Fs=fs, cmap=\"rainbow\")\n",
    " \n",
    "# Set the title of the plot, xlabel and ylabel\n",
    "# and display using show() function\n",
    "plt.title('Spectrogram ')\n",
    "plt.ylim((0,16000))\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Hz\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b6fe1ba",
   "metadata": {},
   "source": [
    "# Table of Contents\n",
    "\n",
    "I. Time Domain Filtering\n",
    "\n",
    "II. DFT for Spectral Analysis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1365759",
   "metadata": {},
   "source": [
    "# I. Time Domain Filtering\n",
    "\n",
    "A discrete signal can be thought of as a function mapping integers to real values (i.e a function $f: \\mathbb{N}\\to\\mathbb{R})$). This is the so-called \"time-domain\" representation of the signal because the integers often represent time in some sense.\n",
    "\n",
    "A system is a process which takes a signal as an input and returns a signal as an output. Digital systems which we use to process signals are called filters. Systems can have several different properties. Two important ones are linearity and time-invariance.\n",
    "\n",
    "**Linearity**: A system $H$ is linear if given input signal $x$, input signal $y$, and scalars $\\alpha$ and $\\beta$, $H(\\alpha x + \\beta y) = \\alpha H(x) + \\beta H(y)$\n",
    "\n",
    "**Time-Invariance**: A system is time-invariant when shifting the input signal in time results in an equal shift in time for the output signal (i.e if $H$ transforms $x[n]$ into $y[n]$, then $H$ transforms $x[n-N]$ into $y[n-N]$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6637af2",
   "metadata": {},
   "source": [
    "## I.a Linear Filtering\n",
    "\n",
    "When a system is Linear and Time-Invariant, we can characterize systems by their impulse response. The impulse response of a system $H$ is given by $h[n] = H(\\delta[n])$ where $$\\delta[n] = \\begin{cases} 1 & \\text{ if, } n=0\\\\ 0 & \\text{ else.} \\end{cases}$$\n",
    "\n",
    "This is useful because it means we can compute the response of the system by doing a **convolution** of the input with the impulse response.\n",
    "\n",
    "$$(x * y)[n] = \\sum_{k=-\\infty}^{\\infty}x[k]y[n-k] $$\n",
    "\n",
    "For example, we can take a moving average by using the filter\n",
    "$$ h_{avg}[n] = \\begin{cases} \\frac{1}{5} & \\text{ if } 0 \\leq n \\leq 4\\\\ 0 & \\text{ else.}\\end{cases} $$\n",
    "\n",
    "We can also define a so-called \"edge detector\" filter in order to detect edges in the audio.\n",
    "$$ h_{edge}[n] = \\begin{cases} (-1)^n & \\text{ if } 0 \\leq n \\leq 1\\\\ 0 & \\text{ else.}\\end{cases} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3edbe33e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edge detector and moving average filters\n",
    "\n",
    "plt.stem(np.linspace(0, 100/44.1, 100), audio[:100])\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()\n",
    "\n",
    "hi_pass = np.array([(-1)**n for n in range(2)])/2.0\n",
    "lo_pass = np.array([1 for n in range(5)])/5.0\n",
    "\n",
    "plt.stem(hi_pass)\n",
    "plt.xlabel(\"samples\")\n",
    "plt.show()\n",
    "\n",
    "plt.stem(lo_pass)\n",
    "plt.xlabel(\"samples\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12ee703f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.stem(np.linspace(0, 100/44.1, 100), np.convolve(audio[:100], hi_pass, \"same\"))\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()\n",
    "\n",
    "plt.stem(np.linspace(0, 100/44.1, 100), np.convolve(audio[:100], lo_pass, \"same\"))\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4be99c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "hi_pass_song = np.convolve(audio, hi_pass)\n",
    "ipd.Audio(hi_pass_song, rate=fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bb4d44f",
   "metadata": {},
   "outputs": [],
   "source": [
    "lo_pass_song = np.convolve(audio, lo_pass)\n",
    "ipd.Audio(lo_pass_song, rate=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f9d2224",
   "metadata": {},
   "source": [
    "## I.b Autocorrelation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13bc42a8",
   "metadata": {},
   "source": [
    "Cross-correlation, in signal processing terms, is the process of convolving one signal with a flipped version of another. Cross-correlation produces a graph of correlation versus time, where correlation is the dot-product of the two signals at that particular point.\n",
    "- If you are looking for the timestamps of a particular noise within a longer signal with other sounds present, you may want to cross-correlate the two.\n",
    "- Cross-correlation is used in sonar to detect when the initial pulse (a known signal) is returned.\n",
    "\n",
    "Autocorrelation is the practice of cross-correlating a signal with itself. It is helpful for eliminating noise, as true frequencies will be preserved due to being periodic, while noise tends to be reduced.\n",
    "\n",
    "$$ r_x(n) = x[n] * x[-n] = \\sum_{k=-\\infty}^{\\infty} x[k] x[n-k] $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6f34d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO: Change this value and see how the noise amplitude affects the signal before and after autocorrelation!\n",
    "noise_amplitude = \n",
    "\n",
    "sample = np.sin(2 * np.pi * np.arange(50)/16)\n",
    "noise = 2*np.random.random(50) - 1\n",
    "noisy_sample = sample+noise_amplitude*noise\n",
    "autocorr = np.convolve(noisy_sample, np.flip(noisy_sample))\n",
    "plt.stem(np.linspace(0,49/44.1, 50), noisy_sample)\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()\n",
    "plt.stem(np.linspace(-49/44.1,50/44.1, 99), autocorr)\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()\n",
    "# ipd.Audio(audio, rate=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e896d34",
   "metadata": {},
   "source": [
    "Autocorrelation will always have a peak in the middle, which will grow larger relative to the rest of the signal the more noisy your signal is. This peak has a strength equal to the overall power of the signal, since it occurs at an offset of zero (meaning the signal is completely overlapping with itself, and the magnitude is $\\sum^{N}_{n=1} X[n]^2$. \n",
    "\n",
    "**Comprehension Question:**\n",
    "\n",
    "However, notice even when the signal is highly corrupted, you can still make out the base frequency in the autocorrelated signal. Why is this?\n",
    "\n",
    "**Answer:**\n",
    "\n",
    "Random noise tends to cancel when there is any offset, but pure frequencies still make it through. This is related to how convolution in the time domain equals multiplication in the frequency domain, and so pure frequencies will stand out above the noise when they are squared. We'll go over more about this later. Thus autocorrelation is often used to denoise signals. \n",
    "\n",
    "**Comprehension question:**\n",
    "\n",
    "Why is the signal contained within a triangular envelope?\n",
    "\n",
    "**Answer:**\n",
    "\n",
    "The length of the autocorrelation goes from -N to N, where N is the number of samples in the original signal. When some points are not overlapping, the non-overlapping points cannot contribute to the signal at all. The window this creates is effectively equivalent to convolving 2 boxes, which makes a triangular envelope."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f70007d3",
   "metadata": {},
   "source": [
    "## I.c Nonlinear Filtering\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0052399f",
   "metadata": {},
   "source": [
    "Sometimes you end up with a signal that has salt and pepper noise (random bits set to 0 or 1) due to corruption or problems with a sensor. Nonlinear filtering, such as median filtering, applies a non-linear filter so that extremely high peaks made by these errors can be filtered out without disproportionately affecting the surrounding signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418ff408",
   "metadata": {},
   "outputs": [],
   "source": [
    "salt_and_pepper = np.random.binomial(1, 0.01, size=audio.shape) * 32000 + np.random.binomial(1, 0.01, size=audio.shape) * -32000\n",
    "audio_corrupted = audio+salt_and_pepper\n",
    "plt.stem(np.linspace(0,999/44.1, 1000), audio_corrupted[:1000])\n",
    "plt.xlabel(\"time (ms)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b95cc7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ipd.Audio(audio_corrupted, rate=fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c3a001",
   "metadata": {},
   "outputs": [],
   "source": [
    "meanfilt = np.array([.2, .2, .2, .2, .2])\n",
    "audio_meanfilt = np.convolve(audio_corrupted, meanfilt)\n",
    "ipd.Audio(audio_meanfilt, rate=fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1008b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# median filter docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.medfilt.html\n",
    "# try changing the next cell to see how different mean and median filtering sound!\n",
    "audio_medfilt = signal.medfilt(audio_corrupted, kernel_size=5)\n",
    "ipd.Audio(audio_medfilt, rate=fs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6219acfe",
   "metadata": {},
   "source": [
    "# II. DFT\n",
    "\n",
    "Typically, when we look at signals, we look at them in the so-called time-domain. Each sample $x[k]$ represents the amplitude of the signal at time-step $k$. This tells us what the signal looks like. One question we might want to ask ourselves is _\"How fast is the signal changing?\"_\n",
    "\n",
    "For sinusoidal signals like $x[n] = \\cos(\\omega n)$ and $x[n] = \\sin(\\omega n)$, answering this question is easy because a larger $\\omega$ means the signal is changing faster ($\\omega$ is known as the angular frequency). For example, consider the plots below which each consist of 100 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "695eb69e",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.linspace(0, 100, 100)\n",
    "slow_cos = np.cos(2 * np.pi * n / 100)\n",
    "fast_cos = np.cos(2 * np.pi * 5 * n / 100)\n",
    "\n",
    "plt.figure(figsize=(15, 7))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.stem(n, slow_cos, use_line_collection=True)\n",
    "plt.title(\"$\\cos\\\\left(\\\\frac{2\\pi}{100} n\\\\right)$\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.title(\"$\\cos\\\\left(\\\\frac{10\\pi}{100} n\\\\right)$\")\n",
    "plt.stem(n, fast_cos, use_line_collection=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6b3299d",
   "metadata": {},
   "source": [
    "$\\cos\\left(\\frac{10\\pi}{100} t\\right)$ is clearly changing a lot faster. If we allow ourselves to consider complex signals, then we can generalized sinusoids using the complex exponential $e^{j\\omega}$. Just like real sinusoids, the $\\omega$ in the signal $x[n] = e^{j\\omega n}$ determines how fast the signal changes (i.e rotates around the unit circle). If we can somehow \"project\" our time-domain signal $x[n]$ onto a \"basis\" of complex exponential signals, then, then the coefficients $X[k]$ should tell us how much the signal changes.\n",
    "\n",
    "The Discrete Fourier Transform is the change of basis which we use for a finite, length-$N$ signal to understand how fast it is changing. The basis used in the DFT are the $N$th roots of unity (i.e the complex solutions to $\\omega=1$). More specifically, the $k$th basis vector is given by $\\phi_k[n] = e^{j\\frac{2\\pi}{N}kn}$. Using the complex inner product $\\langle \\vec{x}, \\vec{y} \\rangle = \\vec{y}^*\\vec{x}$, the DFT coefficients are given by\n",
    "\n",
    "$$X[k] = \\langle x, \\phi_k \\rangle = \\sum_{n=0}^{N-1}x[n]e^{-j\\frac{2\\pi}{N}kn}.$$\n",
    "\n",
    "From the DFT coefficients, we can recover the time-domain coefficients using the inverse DFT.\n",
    "\n",
    "$$x[n] = \\frac{1}{N} \\sum_{k=0}^{N-1}X[k]e^{j\\frac{2\\pi}{N}kn}.$$\n",
    "\n",
    "There are many ways to compute the DFT. The fastest method is the Fast Fourier Transform (FFT), which is an algorithm which computes the DFT. It is built into `numpy` as part of the `fft` submodule.\n",
    "\n",
    "If we look at the DFT coefficients of the two cosines we saw earlier, we can see that it is indeed doing exactly what we wanted it to: characterizing the frequency of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f20f1ae8",
   "metadata": {},
   "outputs": [],
   "source": [
    "slow_cos_fft = np.fft.fft(slow_cos)\n",
    "fast_cos_fft = np.fft.fft(fast_cos)\n",
    "\n",
    "plt.figure(figsize=(15, 7))\n",
    "plt.subplot(2, 2, 1)\n",
    "plt.stem(n, np.abs(slow_cos_fft), use_line_collection=True)\n",
    "plt.title(\"$|DFT\\{\\cos\\\\left(\\\\frac{2\\pi}{100} n\\\\right)\\}|$\")\n",
    "plt.subplot(2, 2, 2)\n",
    "plt.title(\"$|DFT\\{\\cos\\\\left(\\\\frac{10\\pi}{100} n\\\\right)\\}|$\")\n",
    "plt.stem(n, np.abs(fast_cos_fft), use_line_collection=True)\n",
    "plt.subplot(2, 2, 3)\n",
    "plt.stem(n, np.angle(slow_cos_fft), use_line_collection=True)\n",
    "plt.title(\"$\\\\arg \\\\left(DFT\\{\\cos\\\\left(\\\\frac{2\\pi}{100} n\\\\right)\\}\\\\right)$\")\n",
    "plt.subplot(2, 2, 4)\n",
    "plt.title(\"$\\\\arg \\\\left(DFT\\{\\cos\\\\left(\\\\frac{10\\pi}{100} n\\\\right)\\}\\\\right)$\")\n",
    "plt.stem(n, np.angle(fast_cos_fft), use_line_collection=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "908767f9",
   "metadata": {},
   "source": [
    "Since $\\cos\\left(\\frac{2\\pi}{100}n\\right) = \\frac{1}{2}\\left(e^{j\\frac{2\\pi}{100}n} + e^{-j\\frac{2\\pi}{100}n}\\right)$, we should expect peaks at $k = 1$ and $k =-1$ (note that because the roots of unity are periodic, $k=-1$ is the same basis vector as $k=99$). Likewise, since $\\cos\\left(\\frac{10\\pi}{100}n\\right) = \\frac{1}{2}\\left(e^{j\\frac{10\\pi}{100}n} + e^{-j\\frac{10\\pi}{100}n}\\right)$, we should expect peaks at $k=5$ and $k=-5$.\n",
    "\n",
    "There are a few things to note:\n",
    "1. The DFT coefficients are complex numbers, so we need both magnitude (top plots) and phase (bottom plots) to characterize the signal information\n",
    "2. For both $\\cos\\left(\\frac{2\\pi}{100}n\\right)$ and $\\cos\\left(\\frac{10\\pi}{100}n\\right)$, we should only expect 2 non-zero coefficients. However, we have apparently many non-zero coefficients. These are due to numerical instability in the FFT algorithm (if you print them out, these coefficients are on the order of $10^{-3}$ in magnitude and so are insignificant).\n",
    "3. The DFT basis is **not** orthonormal. This is why we must scale by $\\frac{1}{N}$ when applying the inverse DFT (`np.fft.ifft` in numpy). This is also why the peak magnitudes of the example signals above are 50 and not $\\frac{1}{2}$.\n",
    "4. DFT basis vectors are complex conjugates of each other (i.e $\\phi_k[n] = \\phi_{N-k}[n]^*$). This means for real signals, $X[k] = X^*[N-k]$.\n",
    "\n",
    "### Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19bc1309",
   "metadata": {},
   "source": [
    "To get a better feel for the DFT, compute and plot the magnitude of the DFT coefficients of our clip from Tom's Diner in decibels ($dB = 20\\log_{10}(\\cdot)$). Since our song is a real signal, do not plot the complex conjugate coefficients since they are redundant information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c19b7f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 7))\n",
    "\n",
    "#res = 20 * np.log10(np.abs(np.fft.fft(audio)))\n",
    "#plt.plot(res)\n",
    "freqs = np.linspace(0, fs / 2 , audio.size // 2)\n",
    "song_dft = 20 * np.log10(np.abs(np.fft.fft(audio)))\n",
    "plt.plot(freqs, song_dft[:audio.size // 2]) # Coefficents N/2 to N are complex coefficients\n",
    "plt.xlabel(\"Hz\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5642689d",
   "metadata": {},
   "source": [
    "**Comprehension Question**: Do you notice anything interesting about the chart above?\n",
    "\n",
    "**Answer**: Around index 150,000, there is a sharp decline in the magnitude of the DFT coefficients. It turns out that this DFT coefficient represents approximately 12.5 kHz (we'll see how to compute this later), which is close to the human hearing limit of about 20kHz."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc0b0b4f",
   "metadata": {},
   "source": [
    "**Comprehension Question**: What does the first coefficient $X[0]$ of the DFT represent in simple terms?\n",
    "\n",
    "**Answer**: It is the sum of the signal (we can see this from the formula by letting $k=0$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b76755a1",
   "metadata": {},
   "source": [
    "## II.a PSD\n",
    "\n",
    "In signal processing, due to noise, numerical stability, and other issues, we often care about the dominant frequencies in the signal (e.g when we are looking for formants in a vowel). This means we want to look at the magnitude of the DFT coefficients. However, sometimes peaks in the DFT are difficult to distinguish when looking at a magnitude plot. To better distinguish peaks, we can instead look at $|X[k]|^2$, the so-called **Power Spectral Density (PSD)**.\n",
    "\n",
    "The Power Spectral Density is the essentially the magnitude of the DFT of the auto-correlation of the signal $x$. This is because when $x[n]$ has DFT coefficients $X[k]$, then $x[-n]$ has DFT coefficients $X^*[k]$ and since auto-correlation is the convolution of $x[n] * x[-n]$, and convolution in the time-domain is multiplication in the frequency domain, $PSD = X[k] X^*[k] = |X[k]|^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bcdab90",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Remember that formants are the dominant frequencies in vowels. That means we can use the PSD to roughly find formants and distinguish vowels from each other.\n",
    "\n",
    "We have two mystery recordings taken from [this source](https://linguistics.ucla.edu/people/hayes/103/Charts/VChart/). They were sampled at 16000Hz. Try and distingiush them by their dominant frequencies, we will go through the following procedure.\n",
    "\n",
    "1. Split the recording into 25ms sections\n",
    "2. Find the PSD of each section\n",
    "3. Let the \"PSD\" of the recording be the mean value of the PSD of each section at each particular point. This will help pick out the frequencies that are dominant in any section of the recording\n",
    "4. Try and guess where the formants are. Can you tell by their relative positions which vowel each recording is?\n",
    "\n",
    "When plotting, you can use the `FREQS` variable for the x-axis of the plot, and remember, the complex conjugate coefficients of the DFT give you no extra information, so do not plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "630a1541",
   "metadata": {},
   "outputs": [],
   "source": [
    "FREQS = np.linspace(0, 8000, 200)\n",
    "_, vowel_1 = wavfile.read(\"mystery_vowel_1.wav\")\n",
    "_, vowel_2 = wavfile.read(\"mystery_vowel_2.wav\")\n",
    "\n",
    "# Cut each recording to an appropriate length\n",
    "vowel_1 = vowel_1[13:]\n",
    "vowel_2 = vowel_2[114:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1989222c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE #\n",
    "\n",
    "def compute_max_psd(recording):\n",
    "    reshaped = recording.reshape((-1, 400))\n",
    "    psd = np.log10(np.abs(np.fft.fft(reshaped)) ** 2)\n",
    "    max_psd = np.mean(psd, axis=0)\n",
    "    return max_psd[:200]\n",
    "\n",
    "\n",
    "vowel_1_psd = compute_max_psd(vowel_1)\n",
    "vowel_2_psd = compute_max_psd(vowel_2)\n",
    "\n",
    "plt.figure(figsize=(15, 7))\n",
    "\n",
    "# plt.subplot(2, 1, 1)\n",
    "plt.title(\"PSD\")\n",
    "plt.plot(FREQS, vowel_1_psd)\n",
    "plt.plot(FREQS, vowel_2_psd)\n",
    "plt.legend([\"Vowel 1\", \"Vowel 2\"])\n",
    "plt.xlabel(\"Hz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e4a26f9",
   "metadata": {},
   "source": [
    "**Answer**: There is a lot of energy in the lower end of the spectrum for Vowel 2. We can predict that there are two formants in that region since there are no other prominent peaks. That means F1 and F2 are very close together. In contrast, for Vowel 1, we see one peak region below 500Hz and another peak region after 2000Hz. This means Vowel 1 has a relatively low F1 and high F2, making it an \"i\", and Vowel 2 has a relatively high F1 and low F2, making it an \"a\"."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "604671ec",
   "metadata": {},
   "source": [
    "Listen to the audio and see if you were right!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc382851",
   "metadata": {},
   "outputs": [],
   "source": [
    "ipd.Audio(vowel_1, rate=16000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3902df50",
   "metadata": {},
   "outputs": [],
   "source": [
    "ipd.Audio(vowel_2, rate=16000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}