path: root/usth/ICT2.9/practical/filters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RsxAag6zVgGa"
   },
   "source": [
    "# Filter Design"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S-FLJ4AAY14k"
   },
   "source": [
    "## The preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mvme4rlMVkNV"
   },
   "source": [
    "First, we load a sample sound as an example for later experiments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FT2qT4zdVMxC"
   },
   "outputs": [],
   "source": [
    "import wave\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "with wave.open('blue-ocean-floor.wav', 'r') as wav:\n",
    "    frame_rate, sample_width = wav.getframerate(), wav.getsampwidth()\n",
    "    ocean = (lambda a: a/ a.max())(np.frombuffer(\n",
    "        wav.readframes(wav.getnframes()), dtype=f'i{sample_width}'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5_47f9h6WcCM"
   },
   "source": [
    "Then for convenience purposes, we define a few `matplotlib` wrappers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0dkqk-_9WsOq"
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.fft import fft\n",
    "\n",
    "\n",
    "def plt_time(sequence):\n",
    "    \"\"\"Plot the signal in time domain.\"\"\"\n",
    "    length = len(sequence)\n",
    "    plt.plot(np.linspace(0, length/frame_rate, length), sequence)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plt_fft(sequence):\n",
    "    \"\"\"Plot the magnitude of the FT of the signal.\"\"\"\n",
    "    domain = np.linspace(0, pi*2, len(sequence))\n",
    "    plt.plot(domain, np.abs(fft(sequence)))\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plt_db(sequence):\n",
    "    \"\"\"Plot the FT of the signal in dB.\"\"\"\n",
    "    domain = np.linspace(0, pi*2, len(sequence))\n",
    "    FT = fft(sequence)\n",
    "    plt.plot(domain, 20 * np.log10(np.abs(FT/abs(FT).max())))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hzFZU7lxXlrk"
   },
   "source": [
    "The sample audio's waveform can then be visualized using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "-dMPhzl6Xw5R",
    "outputId": "593f98e0-f1d3-4c93-e385-aa8988020276"
   },
   "outputs": [],
   "source": [
    "plt_time(ocean)\n",
    "plt_fft(ocean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Euu6xt44YRG0"
   },
   "source": [
    "It is completely normal that the audio is mainly low frequency, given the `frame_rate` of 44100 Hz, while vocal is around 300 to 1000 Hz and we can usually hear only up to around 10000 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "g_psAgGiYuMN",
    "outputId": "9cda3d29-6fa2-444c-e227-162ff6a0ea30"
   },
   "outputs": [],
   "source": [
    "print(frame_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DiWkQgwXY-Am"
   },
   "source": [
    "## The sinc function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zO-UY4-6ZLAx"
   },
   "source": [
    "The ideal low pass filter has the impulse response of `h[n] = sin(w*n) / (pi*n)` which pass frequencies from 0 to `limit = w/pi * frame_rate` hertz.  For the ease of calculation, we define `r =  w/pi` and get `h[n] = sin(pi*n*r) / (pi*n*r/r) = sinc(n*r) * r`.  We then define our own `sinc` sampler as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "JU9Ra87canaJ"
   },
   "outputs": [],
   "source": [
    "def sinc(limit, length=101):\n",
    "    \"\"\"Return the impulse response of the low-pass filter passing\n",
    "    frequencies up to limit herts, sampled to the given length.\n",
    "    \"\"\"\n",
    "    absolute = length >> 1\n",
    "    n = np.arange(-absolute, absolute+1)\n",
    "    r = limit / frame_rate\n",
    "    return np.sinc(n*r) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VGEAxTCxbOw1"
   },
   "source": [
    "Just to be sure, we plot a filter passing up to 420 Hz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "GiIuofwLbsHG",
    "outputId": "d0847d62-cc59-4004-d7bd-4bf3678d55ca"
   },
   "outputs": [],
   "source": [
    "plt_fft(sinc(3000))\n",
    "plt_db(sinc(3000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WR0SOl2veOiO"
   },
   "source": [
    "In order to have something for the later filters to compare with, we apply it to the input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "yR-2V36OeXtU",
    "outputId": "dc55207b-2fb5-4ae7-fe4d-bd60a34512dd"
   },
   "outputs": [],
   "source": [
    "output_sinc = np.convolve(ocean, sinc(3000))\n",
    "plt_time(output_sinc)\n",
    "plt_fft(output_sinc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Z7PKrdjdb9Xu"
   },
   "source": [
    "## The FIR windows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "brssPyVncKdx"
   },
   "source": [
    "As seen from above, the low-pass filter sampled from `sinc` isn't ideal: there's quite some visible passband ripple and this could distort our output.  Luckily we can improve this by multiply the impulse response with a window, whose functions are provided by `scipy`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6nv8Ickhcwbt"
   },
   "outputs": [],
   "source": [
    "from scipy.signal import windows\n",
    "\n",
    "\n",
    "def fir(limit, window, length=101):\n",
    "    \"\"\"Return a low-pass filter using the specified windowing technique.\"\"\"\n",
    "    return getattr(windows, window)(length) * sinc(limit, length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1poaLOrZdJUl"
   },
   "source": [
    "We first try the Bartlett window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "Itc6VElKhW9g",
    "outputId": "ec91de0d-e007-4d83-84b4-50a63450d4ff"
   },
   "outputs": [],
   "source": [
    "bartlett = fir(3000, 'bartlett')\n",
    "plt_fft(bartlett)\n",
    "plt_db(bartlett)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A9qPm8wAiEDc"
   },
   "source": [
    "Considering only the plot in dB, this seems to be worse than the vanila sinc samples.  Next, we try the Hamming window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "R4Cp5Ya_iUyr",
    "outputId": "8736f38d-4965-4f3c-ab02-1a36902cdd25"
   },
   "outputs": [],
   "source": [
    "hamming = fir(3000, 'hamming')\n",
    "plt_fft(hamming)\n",
    "plt_db(hamming)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WuYjHCFOihPf"
   },
   "source": [
    "The Hamming window seems to provide a smoother filter (i.e. almost no ramples); other than that it's quite equivalent to the original filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "DflFR72-i582",
    "outputId": "b5f2cd8a-2055-4650-a08f-9c4917ed4133"
   },
   "outputs": [],
   "source": [
    "blackman = fir(3000, 'blackman')\n",
    "plt_fft(blackman)\n",
    "plt_db(blackman)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "sNmW3UPVjMNV"
   },
   "source": [
    "The Blackman seems to provide a slightly better result, with the gain dropping noticiably faster in the transition bandwidth."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kal1ldxbpxSh"
   },
   "source": [
    "There is very little point ploting the output signal however, since the original is already too complicated for human eyes.  Pretty much all we will see is very similar to that of the vanila sinc samples.  Instead, we write the output to files and listen to them instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "LuJy6wpdqUAi"
   },
   "outputs": [],
   "source": [
    "def write(signal, filename):\n",
    "    \"\"\"Write the given signal to the corresponding filename.\"\"\"\n",
    "    with wave.open(f'{filename}.wav', 'w') as wav:\n",
    "        wav.setnchannels(1)\n",
    "        wav.setsampwidth(sample_width)\n",
    "        wav.setframerate(frame_rate)\n",
    "        width = sample_width << 3\n",
    "        wav.writeframesraw(getattr(np, f'int{width}')(signal * 2**(width - 1)))\n",
    "\n",
    "\n",
    "write(np.convolve(ocean, sinc(3000)), 'output_sinc')\n",
    "write(np.convolve(ocean, bartlett), 'output_bartlett')\n",
    "write(np.convolve(ocean, hamming), 'output_hamming')\n",
    "write(np.convolve(ocean, blackman), 'output_blackman')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wIlx2i3Mqy2p"
   },
   "source": [
    "After a quick listen, `output_{sinc,bartlett}` seems to have better soundstage than the other two, which means Hamming and Blackman does better jobs filtering *away* high frequencies.  `output_hamming` seems to be the most *dull* to me, which means it's the best filtered one.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3k6HmI3XsTAK"
   },
   "source": [
    "## The IIR filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6tBjWTKNwBKq"
   },
   "source": [
    "In this section, we are going to examine a few IIR filters.  Before we start, let's define a thin wrapper for convenience:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ez0jOaiawNfc"
   },
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "\n",
    "\n",
    "def iir(limit, design, order=4):\n",
    "    \"\"\"Return the specified IIR filter.\"\"\"\n",
    "    if design == 'cheby1': return signal.cheby1(order, 1, limit/frame_rate)\n",
    "    return getattr(signal, design)(order, limit/frame_rate)\n",
    "\n",
    "\n",
    "def plt_iir(b, a):\n",
    "    \"\"\"Plot the given IIR filter in linear as well as log (dB) magnitude.\"\"\"\n",
    "    w, h = signal.freqz(b, a)\n",
    "    plt.plot(w, abs(h))\n",
    "    plt.show()\n",
    "    plt.plot(w, 20*np.log10(abs(h)))\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def apply_iir(data, b, a):\n",
    "    \"\"\"Apply the given IIR filter to data.\"\"\"\n",
    "    return signal.lfilter(b, a, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "iuTzduyKz7EG"
   },
   "source": [
    "We can then respectively plot the Chebyshev, Butterworth and Bessel filters in frequency domain:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "yYawU_7G0B_x",
    "outputId": "29933810-3ae2-4ae2-87da-7a858784009e"
   },
   "outputs": [],
   "source": [
    "cheby1 = iir(3000, 'cheby1')\n",
    "plt_iir(*cheby1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "r2qImpDh1OIs",
    "outputId": "17c3b88e-4e7a-432d-e6f2-8a3d0e7e0477"
   },
   "outputs": [],
   "source": [
    "butter = iir(3000, 'butter')\n",
    "plt_iir(*butter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "1ByMdegV1O9X",
    "outputId": "523fc3e9-09b2-49ee-e626-e667fe3d7abf"
   },
   "outputs": [],
   "source": [
    "bessel = iir(3000, 'bessel')\n",
    "plt_iir(*bessel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zkPBz9lt1eTY"
   },
   "source": [
    "We then apply these to the input a pay listen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "YA0Z9AZP1kI8"
   },
   "outputs": [],
   "source": [
    "write(apply_iir(ocean, *cheby1), 'output_chebyshev')\n",
    "write(apply_iir(ocean, *butter), 'output_butterworth')\n",
    "write(apply_iir(ocean, *bessel), 'output_bessel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "n5u3BNy-3Ys8"
   },
   "source": [
    "As shown in the frequency domain, these filters are not as effective as the ones we discovered in FIR, at least using the configurations above.  The outputs sound a lot clearer which indicates some higher frequency sounds excapes from the filters."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Untitled0.ipynb",
   "provenance": []
  },
  "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.8.3rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}