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

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FTq58_r1KXZd"
   },
   "source": [
    "# DSP Labwork"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nB6RL9TOKm9D"
   },
   "source": [
    "Given the following signals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "rqw23aUkKqS2"
   },
   "outputs": [],
   "source": [
    "input_1kHz_15kHz = [\n",
    "    +0.0000000000, +0.5924659585, -0.0947343455, +0.1913417162, +1.0000000000,\n",
    "    +0.4174197128, +0.3535533906, +1.2552931065, +0.8660254038, +0.4619397663,\n",
    "    +1.3194792169, +1.1827865776, +0.5000000000, +1.1827865776, +1.3194792169,\n",
    "    +0.4619397663, +0.8660254038, +1.2552931065, +0.3535533906, +0.4174197128,\n",
    "    +1.0000000000, +0.1913417162, -0.0947343455, +0.5924659585, -0.0000000000,\n",
    "    -0.5924659585, +0.0947343455, -0.1913417162, -1.0000000000, -0.4174197128,\n",
    "    -0.3535533906, -1.2552931065, -0.8660254038, -0.4619397663, -1.3194792169,\n",
    "    -1.1827865776, -0.5000000000, -1.1827865776, -1.3194792169, -0.4619397663,\n",
    "    -0.8660254038, -1.2552931065, -0.3535533906, -0.4174197128, -1.0000000000,\n",
    "    -0.1913417162, +0.0947343455, -0.5924659585, +0.0000000000, +0.5924659585,\n",
    "    -0.0947343455, +0.1913417162, +1.0000000000, +0.4174197128, +0.3535533906,\n",
    "    +1.2552931065, +0.8660254038, +0.4619397663, +1.3194792169, +1.1827865776,\n",
    "    +0.5000000000, +1.1827865776, +1.3194792169, +0.4619397663, +0.8660254038,\n",
    "    +1.2552931065, +0.3535533906, +0.4174197128, +1.0000000000, +0.1913417162,\n",
    "    -0.0947343455, +0.5924659585, +0.0000000000, -0.5924659585, +0.0947343455,\n",
    "    -0.1913417162, -1.0000000000, -0.4174197128, -0.3535533906, -1.2552931065,\n",
    "    -0.8660254038, -0.4619397663, -1.3194792169, -1.1827865776, -0.5000000000,\n",
    "    -1.1827865776, -1.3194792169, -0.4619397663, -0.8660254038, -1.2552931065,\n",
    "    -0.3535533906, -0.4174197128, -1.0000000000, -0.1913417162, +0.0947343455,\n",
    "    -0.5924659585, +0.0000000000, +0.5924659585, -0.0947343455, +0.1913417162,\n",
    "    +1.0000000000, +0.4174197128, +0.3535533906, +1.2552931065, +0.8660254038,\n",
    "    +0.4619397663, +1.3194792169, +1.1827865776, +0.5000000000, +1.1827865776,\n",
    "    +1.3194792169, +0.4619397663, +0.8660254038, +1.2552931065, +0.3535533906,\n",
    "    +0.4174197128, +1.0000000000, +0.1913417162, -0.0947343455, +0.5924659585,\n",
    "    +0.0000000000, -0.5924659585, +0.0947343455, -0.1913417162, -1.0000000000,\n",
    "    -0.4174197128, -0.3535533906, -1.2552931065, -0.8660254038, -0.4619397663,\n",
    "    -1.3194792169, -1.1827865776, -0.5000000000, -1.1827865776, -1.3194792169,\n",
    "    -0.4619397663, -0.8660254038, -1.2552931065, -0.3535533906, -0.4174197128,\n",
    "    -1.0000000000, -0.1913417162, +0.0947343455, -0.5924659585, -0.0000000000,\n",
    "    +0.5924659585, -0.0947343455, +0.1913417162, +1.0000000000, +0.4174197128,\n",
    "    +0.3535533906, +1.2552931065, +0.8660254038, +0.4619397663, +1.3194792169,\n",
    "    +1.1827865776, +0.5000000000, +1.1827865776, +1.3194792169, +0.4619397663,\n",
    "    +0.8660254038, +1.2552931065, +0.3535533906, +0.4174197128, +1.0000000000,\n",
    "    +0.1913417162, -0.0947343455, +0.5924659585, -0.0000000000, -0.5924659585,\n",
    "    +0.0947343455, -0.1913417162, -1.0000000000, -0.4174197128, -0.3535533906,\n",
    "    -1.2552931065, -0.8660254038, -0.4619397663, -1.3194792169, -1.1827865776,\n",
    "    -0.5000000000, -1.1827865776, -1.3194792169, -0.4619397663, -0.8660254038,\n",
    "    -1.2552931065, -0.3535533906, -0.4174197128, -1.0000000000, -0.1913417162,\n",
    "    +0.0947343455, -0.5924659585, +0.0000000000, +0.5924659585, -0.0947343455,\n",
    "    +0.1913417162, +1.0000000000, +0.4174197128, +0.3535533906, +1.2552931065,\n",
    "    +0.8660254038, +0.4619397663, +1.3194792169, +1.1827865776, +0.5000000000,\n",
    "    +1.1827865776, +1.3194792169, +0.4619397663, +0.8660254038, +1.2552931065,\n",
    "    +0.3535533906, +0.4174197128, +1.0000000000, +0.1913417162, -0.0947343455,\n",
    "    +0.5924659585, +0.0000000000, -0.5924659585, +0.0947343455, -0.1913417162,\n",
    "    -1.0000000000, -0.4174197128, -0.3535533906, -1.2552931065, -0.8660254038,\n",
    "    -0.4619397663, -1.3194792169, -1.1827865776, -0.5000000000, -1.1827865776,\n",
    "    -1.3194792169, -0.4619397663, -0.8660254038, -1.2552931065, -0.3535533906,\n",
    "    -0.4174197128, -1.0000000000, -0.1913417162, +0.0947343455, -0.5924659585,\n",
    "    -0.0000000000, +0.5924659585, -0.0947343455, +0.1913417162, +1.0000000000,\n",
    "    +0.4174197128, +0.3535533906, +1.2552931065, +0.8660254038, +0.4619397663,\n",
    "    +1.3194792169, +1.1827865776, +0.5000000000, +1.1827865776, +1.3194792169,\n",
    "    +0.4619397663, +0.8660254038, +1.2552931065, +0.3535533906, +0.4174197128,\n",
    "    +1.0000000000, +0.1913417162, -0.0947343455, +0.5924659585, +0.0000000000,\n",
    "    -0.5924659585, +0.0947343455, -0.1913417162, -1.0000000000, -0.4174197128,\n",
    "    -0.3535533906, -1.2552931065, -0.8660254038, -0.4619397663, -1.3194792169,\n",
    "    -1.1827865776, -0.5000000000, -1.1827865776, -1.3194792169, -0.4619397663,\n",
    "    -0.8660254038, -1.2552931065, -0.3535533906, -0.4174197128, -1.0000000000,\n",
    "    -0.1913417162, +0.0947343455, -0.5924659585, -0.0000000000, +0.5924659585,\n",
    "    -0.0947343455, +0.1913417162, +1.0000000000, +0.4174197128, +0.3535533906,\n",
    "    +1.2552931065, +0.8660254038, +0.4619397663, +1.3194792169, +1.1827865776,\n",
    "    +0.5000000000, +1.1827865776, +1.3194792169, +0.4619397663, +0.8660254038,\n",
    "    +1.2552931065, +0.3535533906, +0.4174197128, +1.0000000000, +0.1913417162,\n",
    "    -0.0947343455, +0.5924659585, +0.0000000000, -0.5924659585, +0.0947343455,\n",
    "    -0.1913417162, -1.0000000000, -0.4174197128, -0.3535533906, -1.2552931065]\n",
    "impulse_response = [\n",
    "    -0.0018225230, -0.0015879294, +0.0000000000, +0.0036977508, +0.0080754303,\n",
    "    +0.0085302217, -0.0000000000, -0.0173976984, -0.0341458607, -0.0333591565,\n",
    "    +0.0000000000, +0.0676308395, +0.1522061835, +0.2229246956, +0.2504960933,\n",
    "    +0.2229246956, +0.1522061835, +0.0676308395, +0.0000000000, -0.0333591565,\n",
    "    -0.0341458607, -0.0173976984, -0.0000000000, +0.0085302217, +0.0080754303,\n",
    "    +0.0036977508, +0.0000000000, -0.0015879294, -0.0018225230]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Su7BkK-hNrRr"
   },
   "source": [
    "## Signals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7ywQ1k6SMjB1"
   },
   "source": [
    "Then the input signal can be plotted as following\n",
    "(the x-axis is wrong of course, we will tackle this later on).\n",
    "At first glance, it might be obvious (not for me though)\n",
    "that the input is a superposition of two sine waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "NVOkSo_RMcNW",
    "outputId": "5d8993db-acf2-4591-a501-8ece487993b4",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def quickplt(sequence):\n",
    "    \"\"\"Plot the signal as-is.\"\"\"\n",
    "    plt.plot(sequence)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "quickplt(input_1kHz_15kHz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lk2QtF8INxxn"
   },
   "source": [
    "To confirm this suspection, we apply FFT on the signals\n",
    "and plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "id": "PwoTNVzAOBb9",
    "outputId": "3200e5ca-4b15-4ae9-df41-7f13bb31b849",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "\n",
    "import numpy as np\n",
    "from numpy.fft import fft\n",
    "\n",
    "\n",
    "def plt_polar(sequence):\n",
    "    \"\"\"Plot the complex signal in polar coordinate, from 0 to pi*2.\"\"\"\n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)\n",
    "    domain = np.linspace(0, pi*2, len(sequence))\n",
    "    ax1.plot(domain, np.abs(sequence))\n",
    "    ax1.set_title('magnitude')\n",
    "    ax2.plot(domain, np.angle(sequence))\n",
    "    ax2.set_title('phase')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "inpft = fft(input_1kHz_15kHz)\n",
    "plt_polar(inpft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IVXjR_4vchPw"
   },
   "source": [
    "Since the low frequencies are lurking around `k*pi*2` and the high frequencies are around `k*pi*2 + pi`, we can make a good guess that the higher peak is of 1 kHz and the lower one is of 15 kHz.  The sample rate would be at exactly`k*pi*2 + pi` and can be calculated as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "ePfjHH_Adju4",
    "outputId": "cd35dd14-e771-42d3-b6f3-68b1ce6bc96d"
   },
   "outputs": [],
   "source": [
    "sample_rate = (len(inpft)/2) / np.argmax(inpft) * 1000\n",
    "print(sample_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rnHm4xSld9h2"
   },
   "source": [
    "We can then replot the signal with the correct scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "t1zuXQQWeH8o",
    "outputId": "236316d3-89ff-4e0f-88e0-0f7e529543f5",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plt_time(sequence):\n",
    "    \"\"\"Plot the signal in time domain.\"\"\"\n",
    "    length = len(sequence)\n",
    "    plt.plot(np.linspace(0, length/sample_rate, length), sequence)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "plt_time(input_1kHz_15kHz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WscxWf4PP0dN"
   },
   "source": [
    "The plot of these in cartesian coordinates doesn't give me any further understanding, however:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "id": "hyV1DJzaQLad",
    "outputId": "68b6097b-0d80-4fd2-babb-0dd3efbfa38e"
   },
   "outputs": [],
   "source": [
    "def plt_rect(sequence):\n",
    "    \"\"\"Plot the complex signal in rectangular coordinate, from 0 to pi*2.\"\"\" \n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)\n",
    "    domain = np.linspace(0, pi*2, len(sequence))\n",
    "    ax1.plot(domain, np.real(sequence))\n",
    "    ax1.set_title('real')\n",
    "    ax2.plot(domain, np.imag(sequence))\n",
    "    ax2.set_title('imaginary')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "plt_rect(inpft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-IiQK8jyUWaJ"
   },
   "source": [
    "## Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "p9V5bnZ4-E2K"
   },
   "source": [
    "### Low-pass Filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Vgs77T9yUhAU"
   },
   "source": [
    "In this section, we also try to do the same thing\n",
    "for the impulse response, which seems to be a sinc function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 793
    },
    "colab_type": "code",
    "id": "fh5OzwtRUv54",
    "outputId": "727ad099-50c3-4ebc-bb50-7fd096d55ed3"
   },
   "outputs": [],
   "source": [
    "quickplt(impulse_response)\n",
    "lfft = fft(impulse_response)\n",
    "plt_polar(lfft)\n",
    "plt_rect(lfft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3VGWOQv-VREF"
   },
   "source": [
    "As shown in the graphs, the system is a low-pass filter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "q70T5CckfPwf"
   },
   "source": [
    "Applying the system to the input we get what is undeniably the sinusoidal signal of frequency of 1 kHz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "153JptX0fPYL",
    "outputId": "8b688ebd-10bf-408b-9d40-6e0ed1c633d3"
   },
   "outputs": [],
   "source": [
    "output = np.convolve(input_1kHz_15kHz, impulse_response)\n",
    "plt_time(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "q0m_kWbQfy-1"
   },
   "source": [
    "Alternative to convolution in time domain, we can multiply the FTs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "EFZ0F9gTivdG",
    "outputId": "536337f3-821e-4f63-ea2d-c62c966003a7"
   },
   "outputs": [],
   "source": [
    "from numpy.fft import ifft\n",
    "from scipy import interpolate\n",
    "\n",
    "f = interpolate.interp1d(np.linspace(0, pi*2, len(lfft)), lfft, kind='zero')\n",
    "plt_time(np.real(ifft(inpft*f(np.linspace(0, pi*2, len(inpft))))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eoTRi4pfj3mD"
   },
   "source": [
    "It is noticeable that the wave is now distorted in shape.  Funny enough, using other interpolation methods, the result is much worse (using the convoluted one as the reference), for example the linear one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "4JbRKH-zkGZx",
    "outputId": "715bdc0a-f48f-4a26-b20e-095516f5bf61"
   },
   "outputs": [],
   "source": [
    "f = interpolate.interp1d(np.linspace(0, pi*2, len(lfft)), lfft, kind='linear')\n",
    "plt_time(np.real(ifft(inpft*f(np.linspace(0, pi*2, len(inpft))))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YP4VVonnyIfq"
   },
   "source": [
    "Notice that the low-pass filter filtered out the 15 kHz wave:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 545
    },
    "colab_type": "code",
    "id": "c7_k8qUAyLTi",
    "outputId": "73ac8110-c9b8-4ec6-cbf8-d3907fd72d8b"
   },
   "outputs": [],
   "source": [
    "plt_polar(fft(output))\n",
    "plt_rect(fft(output))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ywBqbDYo-O8c"
   },
   "source": [
    "### High-pass filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WTJ256AU-UeI"
   },
   "source": [
    "To turn the given low-pass filter to a high-pass one,\n",
    "we subtract it from the impulse signal (which is equivalent to\n",
    "subtracting it from 1 in the frequency domain thanks to linearity):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 793
    },
    "colab_type": "code",
    "id": "TGtdsdLi-jmX",
    "outputId": "f9720e14-eb43-4368-b514-07cd7d1e8a76"
   },
   "outputs": [],
   "source": [
    "high_pass = ((lambda m: [0]*m + [1] + [0]*m)(np.argmax(impulse_response))\n",
    "             - np.array(impulse_response))\n",
    "quickplt(high_pass)\n",
    "plt_polar(fft(high_pass))\n",
    "plt_rect(fft(high_pass))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ScsGWwmIAP6E"
   },
   "source": [
    "We then apply it to the input and get the high frequency signal of 15 kHz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 793
    },
    "colab_type": "code",
    "id": "OVMsNKwgAk6b",
    "outputId": "6810e7a3-8db3-4853-e19d-30551a4e4ece"
   },
   "outputs": [],
   "source": [
    "outputhf = np.convolve(input_1kHz_15kHz, high_pass)\n",
    "plt_time(outputhf)\n",
    "plt_polar(fft(outputhf))\n",
    "plt_rect(fft(outputhf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JUMVl3wrBAFX"
   },
   "source": [
    "Other methods to create a HF filter have been tried, however the result is nowhere as good:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nAg7C28hBNKd"
   },
   "source": [
    "1. Shifting the low-pass filter by pi in frequency domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 545
    },
    "colab_type": "code",
    "id": "Cib04NkFBIVD",
    "outputId": "9166d52e-c8e4-49f7-e2aa-ddf1978cb885"
   },
   "outputs": [],
   "source": [
    "high_pass_bad = ifft(np.roll(lfft, len(impulse_response)>>1))\n",
    "outputhf_bad = np.convolve(input_1kHz_15kHz, high_pass_bad)\n",
    "plt_rect(outputhf_bad)\n",
    "plt_polar(fft(outputhf_bad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "CbW0fDIQCTfP"
   },
   "source": [
    "2. Multiply the low-pass in time domain with (-1)^n (which effectively also shift it by pi):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 529
    },
    "colab_type": "code",
    "id": "3UYAyK6eC-Yg",
    "outputId": "7d0c08c1-589c-4326-bc25-17a29c457ed9"
   },
   "outputs": [],
   "source": [
    "high_pass_worse = np.fromiter(((-1)**n for n in range(len(impulse_response))),\n",
    "                              dtype=float) * impulse_response\n",
    "outputhf_worse = np.convolve(input_1kHz_15kHz, high_pass_worse)\n",
    "plt_time(outputhf_worse)\n",
    "plt_polar(fft(outputhf_worse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "l1qP70KvDvIO"
   },
   "source": [
    "While both of these produce a high frequency signal for most of the interval, at the start and end the *volume* is significantly higher and there are many different frequencies instead of just 15 kHz.  This seems disagrees with the theory at first, but the theory is only supposed to apply to infinite length impulse response; manipulate a rather small in size sample cannot give perfect results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "p2H5i0CX5OSY"
   },
   "source": [
    "## Electrocardiogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zf8MKQ_p5b5t"
   },
   "source": [
    "Given the following ECG signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 793
    },
    "colab_type": "code",
    "id": "ulTRe5DT5rNV",
    "outputId": "25dae93f-6ffe-4e13-dcc1-de77564d87b8"
   },
   "outputs": [],
   "source": [
    "ecg = [\n",
    "    0, 0.0010593, 0.0021186, 0.003178, 0.0042373, 0.0052966, 0.0063559,\n",
    "    0.0074153, 0.0084746, 0.045198, 0.081921, 0.11864, 0.15537, 0.19209,\n",
    "    0.22881, 0.26554, 0.30226, 0.33898, 0.30226, 0.26554, 0.22881, 0.19209,\n",
    "    0.15537, 0.11864, 0.081921, 0.045198, 0.0084746, 0.0077684, 0.0070621,\n",
    "    0.0063559, 0.0056497, 0.0049435, 0.0042373, 0.0035311, 0.0028249,\n",
    "    0.0021186, 0.0014124, 0.00070621, 0, -0.096045, -0.19209, -0.28814,\n",
    "    -0.073446, 0.14124, 0.35593, 0.57062, 0.78531, 1, 0.73729, 0.47458,\n",
    "    0.21186, -0.050847, -0.31356, -0.57627, -0.83898, -0.55932, -0.27966, 0,\n",
    "    0.00073692, 0.0014738, 0.0022108, 0.0029477, 0.0036846, 0.0044215,\n",
    "    0.0051584, 0.0058954, 0.0066323, 0.0073692, 0.0081061, 0.008843, 0.00958,\n",
    "    0.010317, 0.011054, 0.011791, 0.012528, 0.013265, 0.014001, 0.014738,\n",
    "    0.015475, 0.016212, 0.016949, 0.03484, 0.052731, 0.070621, 0.088512,\n",
    "    0.1064, 0.12429, 0.14218, 0.16008, 0.17797, 0.16186, 0.14576, 0.12966,\n",
    "    0.11356, 0.097458, 0.081356, 0.065254, 0.049153, 0.033051, 0.016949,\n",
    "    0.013559, 0.010169, 0.0067797, 0.0033898, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
    "    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0010593, 0.0021186, 0.003178,\n",
    "    0.0042373, 0.0052966, 0.0063559, 0.0074153, 0.0084746, 0.045198, 0.081921,\n",
    "    0.11864, 0.15537, 0.19209, 0.22881, 0.26554, 0.30226, 0.33898, 0.30226,\n",
    "    0.26554, 0.22881, 0.19209, 0.15537, 0.11864, 0.081921, 0.045198, 0.0084746,\n",
    "    0.0077684, 0.0070621, 0.0063559, 0.0056497, 0.0049435, 0.0042373,\n",
    "    0.0035311, 0.0028249, 0.0021186, 0.0014124, 0.00070621, 0, -0.096045,\n",
    "    -0.19209, -0.28814, -0.073446, 0.14124, 0.35593, 0.57062, 0.78531, 1,\n",
    "    0.73729, 0.47458, 0.21186, -0.050847, -0.31356, -0.57627, -0.83898,\n",
    "    -0.55932, -0.27966, 0, 0.00073692, 0.0014738, 0.0022108, 0.0029477,\n",
    "    0.0036846, 0.0044215, 0.0051584, 0.0058954, 0.0066323, 0.0073692,\n",
    "    0.0081061, 0.008843, 0.00958, 0.010317, 0.011054, 0.011791, 0.012528,\n",
    "    0.013265, 0.014001, 0.014738, 0.015475, 0.016212, 0.016949, 0.03484,\n",
    "    0.052731, 0.070621, 0.088512, 0.1064, 0.12429, 0.14218, 0.16008, 0.17797,\n",
    "    0.16186, 0.14576, 0.12966, 0.11356, 0.097458, 0.081356, 0.065254, 0.049153,\n",
    "    0.033051, 0.016949, 0.013559, 0.010169, 0.0067797, 0.0033898, 0, 0, 0, 0,\n",
    "    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0010593,\n",
    "    0.0021186, 0.003178, 0.0042373, 0.0052966, 0.0063559, 0.0074153, 0.0084746,\n",
    "    0.045198, 0.081921, 0.11864, 0.15537, 0.19209, 0.22881, 0.26554, 0.30226,\n",
    "    0.33898, 0.30226, 0.26554, 0.22881, 0.19209, 0.15537, 0.11864, 0.081921,\n",
    "    0.045198, 0.0084746, 0.0077684, 0.0070621, 0.0063559, 0.0056497, 0.0049435,\n",
    "    0.0042373, 0.0035311, 0.0028249, 0.0021186, 0.0014124, 0.00070621, 0,\n",
    "    -0.096045, -0.19209, -0.28814, -0.073446, 0.14124, 0.35593, 0.57062,\n",
    "    0.78531, 1, 0.73729, 0.47458, 0.21186, -0.050847, -0.31356, -0.57627,\n",
    "    -0.83898, -0.55932, -0.27966, 0, 0.00073692, 0.0014738, 0.0022108,\n",
    "    0.0029477, 0.0036846, 0.0044215, 0.0051584, 0.0058954, 0.0066323,\n",
    "    0.0073692, 0.0081061, 0.008843, 0.00958, 0.010317, 0.011054, 0.011791,\n",
    "    0.012528, 0.013265, 0.014001, 0.014738, 0.015475, 0.016212, 0.016949,\n",
    "    0.03484, 0.052731, 0.070621, 0.088512, 0.1064, 0.12429, 0.14218, 0.16008,\n",
    "    0.17797, 0.16186, 0.14576, 0.12966, 0.11356, 0.097458, 0.081356, 0.065254,\n",
    "    0.049153, 0.033051, 0.016949, 0.013559, 0.010169, 0.0067797, 0.0033898, 0,\n",
    "    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
    "    0.0010593, 0.0021186, 0.003178, 0.0042373, 0.0052966, 0.0063559, 0.0074153,\n",
    "    0.0084746, 0.045198, 0.081921, 0.11864, 0.15537, 0.19209, 0.22881, 0.26554,\n",
    "    0.30226, 0.33898, 0.30226, 0.26554, 0.22881, 0.19209, 0.15537, 0.11864,\n",
    "    0.081921, 0.045198, 0.0084746, 0.0077684, 0.0070621, 0.0063559, 0.0056497,\n",
    "    0.0049435, 0.0042373, 0.0035311, 0.0028249, 0.0021186, 0.0014124,\n",
    "    0.00070621, 0, -0.096045, -0.19209, -0.28814, -0.073446, 0.14124, 0.35593,\n",
    "    0.57062, 0.78531, 1, 0.73729, 0.47458, 0.21186, -0.050847, -0.31356,\n",
    "    -0.57627, -0.83898, -0.55932, -0.27966, 0, 0.00073692, 0.0014738,\n",
    "    0.0022108, 0.0029477, 0.0036846, 0.0044215, 0.0051584, 0.0058954,\n",
    "    0.0066323, 0.0073692, 0.0081061, 0.008843, 0.00958, 0.010317, 0.011054,\n",
    "    0.011791, 0.012528, 0.013265, 0.014001, 0.014738, 0.015475, 0.016212,\n",
    "    0.016949, 0.03484, 0.052731, 0.070621, 0.088512, 0.1064, 0.12429, 0.14218,\n",
    "    0.16008, 0.17797, 0.16186, 0.14576, 0.12966, 0.11356, 0.097458, 0.081356,\n",
    "    0.065254, 0.049153, 0.033051, 0.016949, 0.013559, 0.010169, 0.0067797,\n",
    "    0.0033898, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
    "    0, 0, 0, 0, 0.0010593, 0.0021186, 0.003178, 0.0042373, 0.0052966,\n",
    "    0.0063559, 0.0074153, 0.0084746, 0.045198, 0.081921, 0.11864, 0.15537,\n",
    "    0.19209, 0.22881, 0.26554, 0.30226, 0.33898, 0.30226, 0.26554, 0.22881,\n",
    "    0.19209, 0.15537, 0.11864, 0.081921, 0.045198, 0.0084746, 0.0077684,\n",
    "    0.0070621, 0.0063559, 0.0056497, 0.0049435, 0.0042373, 0.0035311,\n",
    "    0.0028249, 0.0021186, 0.0014124, 0.00070621, 0, -0.096045, -0.19209,\n",
    "    -0.28814, -0.073446, 0.14124, 0.35593, 0.57062, 0.78531, 1, 0.73729,\n",
    "    0.47458, 0.21186, -0.050847, -0.31356, -0.57627, -0.83898, -0.55932,\n",
    "    -0.27966, 0, 0.00073692, 0.0014738, 0.0022108, 0.0029477, 0.0036846,\n",
    "    0.0044215, 0.0051584, 0.0058954, 0.0066323, 0.0073692, 0.0081061, 0.008843,\n",
    "    0.00958, 0.010317, 0.011054, 0.011791, 0.012528, 0.013265, 0.014001,\n",
    "    0.014738, 0.015475, 0.016212, 0.016949, 0.03484, 0.052731, 0.070621,\n",
    "    0.088512, 0.1064, 0.12429, 0.14218, 0.16008, 0.17797, 0.16186, 0.14576,\n",
    "    0.12966, 0.11356, 0.097458, 0.081356, 0.065254, 0.049153, 0.033051,\n",
    "    0.016949, 0.013559, 0.010169, 0.0067797, 0.0033898, 0, 0, 0, 0, 0, 0, 0, 0,\n",
    "    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
    "quickplt(ecg)\n",
    "plt_polar(fft(ecg))\n",
    "plt_rect(fft(ecg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Aw8a-DUZFS6m"
   },
   "source": [
    "The plots of the ECG gives us the initial intuition that it's a rather low frequency signal (in fact we have the heartrate of somewhere between 50 to 500 Hz) with multiple face each priod.  This is confirmed by the low-passed response, which looks surprisingly similar to the original in time domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 793
    },
    "colab_type": "code",
    "id": "uTj33b3YF2Ij",
    "outputId": "f1677908-82a9-47c3-cb40-f2cbbdf7068d"
   },
   "outputs": [],
   "source": [
    "ecglf = np.convolve(ecg, impulse_response)\n",
    "quickplt(ecglf)\n",
    "plt_polar(fft(ecglf))\n",
    "plt_rect(fft(ecglf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EmUprtsbGVdT"
   },
   "source": [
    "Just for fun, we add some high frequency noise to the signal (because heartbeats are repetitive and boring!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "UEQt8JUpGjdB",
    "outputId": "8166fb68-294d-4e35-b2db-bcfeb1d3535f"
   },
   "outputs": [],
   "source": [
    "real, imag = np.random.random((2, len(ecg)-len(high_pass)+1))\n",
    "whitenoise = ifft(real + imag*1j)\n",
    "noise_hf = abs(np.convolve(whitenoise, high_pass))\n",
    "quickplt(noise_hf)\n",
    "noisy_ecg = ecg + noise_hf\n",
    "quickplt(noisy_ecg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1BHX-N2qG_Cx"
   },
   "source": [
    "To smoothen the noisy signal back to normal, the low pass filter shoud be able to does the job:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "7TpUggulHG0b",
    "outputId": "bbdd7926-b39d-4344-fe4e-9af6531f9aa7"
   },
   "outputs": [],
   "source": [
    "recovered_ecg = np.convolve(noisy_ecg, impulse_response)\n",
    "quickplt(recovered_ecg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0G7413ANVjDq"
   },
   "source": [
    "## Bonus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pwceV47wVoNB"
   },
   "source": [
    "In this section, we'll try to have some fun playing the signals\n",
    "using [palace](https://pypi.org/project/palace/).  While the main purpose\n",
    "of palace is positional audio rendering and environmental effect,\n",
    "it also provide a handy decoder base class, which can be easily derived:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 124
    },
    "colab_type": "code",
    "id": "Q5rVNQDeKMbQ",
    "outputId": "b5f12c91-391b-427b-c73c-83cc2c3eea6a",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!pip install palace\n",
    "from palace import BaseDecoder, Buffer, Context, Device\n",
    "\n",
    "\n",
    "class Dec(BaseDecoder):\n",
    "    \"\"\"Generator of elementary signals.\"\"\"\n",
    "    def __init__(self, content):\n",
    "        self.content, self.size = content.copy(), len(content)\n",
    "\n",
    "    @BaseDecoder.frequency.getter\n",
    "    def frequency(self) -> int: return int(sample_rate)\n",
    "\n",
    "    @BaseDecoder.channel_config.getter\n",
    "    def channel_config(self) -> str:\n",
    "        return 'Mono'\n",
    "\n",
    "    @BaseDecoder.sample_type.getter\n",
    "    def sample_type(self) -> str:\n",
    "        return '32-bit float'\n",
    "\n",
    "    @BaseDecoder.length.getter\n",
    "    def length(self) -> int: return self.size\n",
    "\n",
    "    def seek(self, pos: int) -> bool: return False\n",
    "\n",
    "    @BaseDecoder.loop_points.getter\n",
    "    def loop_points(self): return 0, 0\n",
    "\n",
    "    def read(self, count: int) -> bytes:\n",
    "        if count > len(self.content):\n",
    "            try:\n",
    "                return np.float32(self.content).tobytes()\n",
    "            finally:\n",
    "                self.content = []\n",
    "        data, self.content = self.content[:count], self.content[count:]\n",
    "        return np.float32(data).tobytes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cP3jOX1IlQvc"
   },
   "source": [
    "The input and output signals can then be played by running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 271
    },
    "colab_type": "code",
    "id": "vbgVYiGelOft",
    "outputId": "91bc61f4-baaf-4ecb-e4c2-315dc661e33a"
   },
   "outputs": [],
   "source": [
    "from time import sleep\n",
    "\n",
    "with Device() as d, Context(d) as c:\n",
    "    with Buffer.from_decoder(Dec(input_1kHz_15kHz), 'input') as b, b.play() as s:\n",
    "        sleep(1)\n",
    "    with Buffer.from_decoder(Dec(output), 'lf') as b, b.play() as s:\n",
    "        sleep(1)\n",
    "    with Buffer.from_decoder(Dec(outputhf), 'hf') as b, b.play() as s:\n",
    "        sleep(1)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "dsp.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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
}