From 9b82ec25a2915273b1a6dbd2e9fbb66b0a75587d Mon Sep 17 00:00:00 2001 From: Nguyễn Gia Phong Date: Fri, 12 Jun 2020 14:23:00 +0700 Subject: [usth/ICT2.9] Add final project --- .gitmodules | 4 + README.md | 3 +- usth/ICT2.9/anon.pdf | Bin 95868 -> 0 bytes usth/ICT2.9/anon.tex | 81 --- usth/ICT2.9/blue-ocean-floor.wav | Bin 740622 -> 0 bytes usth/ICT2.9/compose/anon.pdf | Bin 0 -> 95868 bytes usth/ICT2.9/compose/anon.tex | 81 +++ usth/ICT2.9/dsp.ipynb | 981 --------------------------- usth/ICT2.9/filters.ipynb | 613 ----------------- usth/ICT2.9/output_bartlett.wav | Bin 740222 -> 0 bytes usth/ICT2.9/output_bessel.wav | Bin 740022 -> 0 bytes usth/ICT2.9/output_blackman.wav | Bin 740222 -> 0 bytes usth/ICT2.9/output_butterworth.wav | Bin 740022 -> 0 bytes usth/ICT2.9/output_chebyshev.wav | Bin 740022 -> 0 bytes usth/ICT2.9/output_hamming.wav | Bin 740222 -> 0 bytes usth/ICT2.9/output_sinc.wav | Bin 740222 -> 0 bytes usth/ICT2.9/practical/blue-ocean-floor.wav | Bin 0 -> 740622 bytes usth/ICT2.9/practical/dsp.ipynb | 981 +++++++++++++++++++++++++++ usth/ICT2.9/practical/filters.ipynb | 613 +++++++++++++++++ usth/ICT2.9/practical/output_bartlett.wav | Bin 0 -> 740222 bytes usth/ICT2.9/practical/output_bessel.wav | Bin 0 -> 740022 bytes usth/ICT2.9/practical/output_blackman.wav | Bin 0 -> 740222 bytes usth/ICT2.9/practical/output_butterworth.wav | Bin 0 -> 740022 bytes usth/ICT2.9/practical/output_chebyshev.wav | Bin 0 -> 740022 bytes usth/ICT2.9/practical/output_hamming.wav | Bin 0 -> 740222 bytes usth/ICT2.9/practical/output_sinc.wav | Bin 0 -> 740222 bytes usth/ICT2.9/speakerid | 1 + 27 files changed, 1681 insertions(+), 1677 deletions(-) delete mode 100644 usth/ICT2.9/anon.pdf delete mode 100644 usth/ICT2.9/anon.tex delete mode 100644 usth/ICT2.9/blue-ocean-floor.wav create mode 100644 usth/ICT2.9/compose/anon.pdf create mode 100644 usth/ICT2.9/compose/anon.tex delete mode 100644 usth/ICT2.9/dsp.ipynb delete mode 100644 usth/ICT2.9/filters.ipynb delete mode 100644 usth/ICT2.9/output_bartlett.wav delete mode 100644 usth/ICT2.9/output_bessel.wav delete mode 100644 usth/ICT2.9/output_blackman.wav delete mode 100644 usth/ICT2.9/output_butterworth.wav delete mode 100644 usth/ICT2.9/output_chebyshev.wav delete mode 100644 usth/ICT2.9/output_hamming.wav delete mode 100644 usth/ICT2.9/output_sinc.wav create mode 100644 usth/ICT2.9/practical/blue-ocean-floor.wav create mode 100644 usth/ICT2.9/practical/dsp.ipynb create mode 100644 usth/ICT2.9/practical/filters.ipynb create mode 100644 usth/ICT2.9/practical/output_bartlett.wav create mode 100644 usth/ICT2.9/practical/output_bessel.wav create mode 100644 usth/ICT2.9/practical/output_blackman.wav create mode 100644 usth/ICT2.9/practical/output_butterworth.wav create mode 100644 usth/ICT2.9/practical/output_chebyshev.wav create mode 100644 usth/ICT2.9/practical/output_hamming.wav create mode 100644 usth/ICT2.9/practical/output_sinc.wav create mode 160000 usth/ICT2.9/speakerid diff --git a/.gitmodules b/.gitmodules index 0bc314b..affdd48 100644 --- a/.gitmodules +++ b/.gitmodules @@ -5,3 +5,7 @@ [submodule "usth/ICT2.8"] path = usth/ICT2.8 url = https://github.com/McSinyx/os2020.git + +[submodule "usth/ICT2.9/speakerid"] + path = usth/ICT2.9/speakerid + url = https://github.com/McSinyx/speakerid.git diff --git a/README.md b/README.md index dd9b524..839d24f 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ this README as well as commit messages are duolingo (anglais et vietnamien). | Thư mục | Nguồn đề bài | | ------------ | ------------------------------------------------------ | -| `\d{2}` | Đề thi, kiểm tra phân theo lớp | +| `2ndary` | Secondary school competitions | | `coci` | [Giải Tin học Croatia mở rộng][0] | | `codechef` | [Codechef][2] | | `codeforces` | [Codeforces][3] | @@ -18,7 +18,6 @@ this README as well as commit messages are duolingo (anglais et vietnamien). | `paip` | Paradigms of Artificial Intelligence Programming | | `sicp` | Structure and Interpretation of Computer Programs | | `thinkperl6` | Think Perl 6 | -| `tht` | Hội thi Tin học trẻ | | `toys` | Programs that don't deserve their own repo | | `usth` | L'Université des Sciences et des Technologies de Hanoï | diff --git a/usth/ICT2.9/anon.pdf b/usth/ICT2.9/anon.pdf deleted file mode 100644 index 3416256..0000000 Binary files a/usth/ICT2.9/anon.pdf and /dev/null differ diff --git a/usth/ICT2.9/anon.tex b/usth/ICT2.9/anon.tex deleted file mode 100644 index 4a3b62f..0000000 --- a/usth/ICT2.9/anon.tex +++ /dev/null @@ -1,81 +0,0 @@ -\documentclass[a4paper,12pt]{article} -\usepackage[english,vietnamese]{babel} -\usepackage{amsmath} -\usepackage{enumerate} -\usepackage{lmodern} - -\title{System Cascade Connection} -\author{{\selectlanguage{vietnamese}Nguyễn Gia Phong}} - -\begin{document} -\selectlanguage{english}\maketitle -Given two discrete-time systems $A$ and $B$ connected in cascade to form -a new system $C = x \mapsto B(A(x))$. - -\section{Linearity} -If $A$ and $B$ are linear, i.e. for all signals $x_i$ and scalars $a_i$, -\begin{align*} - A\left(n \mapsto \sum_i a_i x_i[n]\right) &= n \mapsto \sum_i a_i A(x_i)[n]\\ - B\left(n \mapsto \sum_i a_i x_i[n]\right) &= n \mapsto \sum_i a_i B(x_i)[n] -\end{align*} -then $C$ is also linear -\begin{align*} - C\left(n \mapsto \sum_i a_i x_i[n]\right) - &= B\left(A\left(n \mapsto \sum_i a_i x_i[n]\right)\right)\\ - &= B\left(n \mapsto \sum_i a_i A(x_i)[n]\right)\\ - &= n \mapsto \sum_i a_i B(A(x_i))[n]\\ - &= n \mapsto \sum_i a_i C(x_i)[n] -\end{align*} - -\section{Time Invariance} -If $A$ and $B$ are time invariant, i.e. for all signals $x$ and integers $k$, -\begin{align*} - A(n \mapsto x[n - k]) &= n \mapsto A(x)[n - k]\\ - B(n \mapsto x[n - k]) &= n \mapsto B(x)[n - k]\\ -\end{align*} -then $C$ is also time invariant -\begin{align*} - C(n \mapsto x[n - k]) - &= B(A(n \mapsto x[n - k]))\\ - &= B(n \mapsto A(x)[n - k])\\ - &= n \mapsto B(A(x))[n - k]\\ - &= n \mapsto C(x)[n - k] -\end{align*} - -\section{LTI Ordering} -If $A$ and $B$ are linear and time-invariant, there exists signals $g$ and $h$ -such that for all signals $x$, $A = x \mapsto x * g$ and $B = x \mapsto x * h$, -thus \[B(A(x)) = B(x * g) = x * g * h = x * h * g = A(x * h) = A(B(x))\] -or interchanging $A$ and $B$ order does not change $C$. - -\section{Causality} -If $A$ and $B$ are causal, i.e. for all signals $x$, $y$ and integers $k$, -\begin{multline*} - x[n] = y[n]\quad\forall n < k - \Longrightarrow\begin{cases} - A(x)[n] = A(y)[n]\quad\forall n < k\\ - B(x)[n] = B(y)[n]\quad\forall n < k - \end{cases}\\ - \Longrightarrow B(A(x))[n] = B(A(y))[n]\quad\forall n < k - \iff C(x)[n] = C(y)[n]\quad\forall n < k -\end{multline*} -then $C$ is also causal. - -\section{BIBO Stability} -If $A$ and $B$ are stable, i.e. there exists a signal $x$ -and scalars $a$, $b$ that -\begin{align*} - |x[n]| < a\quad\forall n \in \mathbb Z - &\Longrightarrow |A(x)[n]| < b\quad\forall n \in \mathbb Z\\ - |x[n]| < a\quad\forall n \in \mathbb Z - &\Longrightarrow |B(x)[n]| < b\quad\forall n \in \mathbb Z -\end{align*} -then $C$ is also stable, i.e. there exists a signal $x$ -and scalars $a$, $b$, $c$ that -\begin{align*} - |x[n]| < a\;\forall n \in \mathbb Z - &\Longrightarrow |A(x)[n]| < b\;\forall n \in \mathbb Z\\ - &\Longrightarrow |B(A(x))[n]| < c\;\forall n \in \mathbb Z - \iff |C(x)[n]| < c\;\forall n \in \mathbb Z -\end{align*} -\end{document} diff --git a/usth/ICT2.9/blue-ocean-floor.wav b/usth/ICT2.9/blue-ocean-floor.wav deleted file mode 100644 index 4208cd0..0000000 Binary files a/usth/ICT2.9/blue-ocean-floor.wav and /dev/null differ diff --git a/usth/ICT2.9/compose/anon.pdf b/usth/ICT2.9/compose/anon.pdf new file mode 100644 index 0000000..3416256 Binary files /dev/null and b/usth/ICT2.9/compose/anon.pdf differ diff --git a/usth/ICT2.9/compose/anon.tex b/usth/ICT2.9/compose/anon.tex new file mode 100644 index 0000000..4a3b62f --- /dev/null +++ b/usth/ICT2.9/compose/anon.tex @@ -0,0 +1,81 @@ +\documentclass[a4paper,12pt]{article} +\usepackage[english,vietnamese]{babel} +\usepackage{amsmath} +\usepackage{enumerate} +\usepackage{lmodern} + +\title{System Cascade Connection} +\author{{\selectlanguage{vietnamese}Nguyễn Gia Phong}} + +\begin{document} +\selectlanguage{english}\maketitle +Given two discrete-time systems $A$ and $B$ connected in cascade to form +a new system $C = x \mapsto B(A(x))$. + +\section{Linearity} +If $A$ and $B$ are linear, i.e. for all signals $x_i$ and scalars $a_i$, +\begin{align*} + A\left(n \mapsto \sum_i a_i x_i[n]\right) &= n \mapsto \sum_i a_i A(x_i)[n]\\ + B\left(n \mapsto \sum_i a_i x_i[n]\right) &= n \mapsto \sum_i a_i B(x_i)[n] +\end{align*} +then $C$ is also linear +\begin{align*} + C\left(n \mapsto \sum_i a_i x_i[n]\right) + &= B\left(A\left(n \mapsto \sum_i a_i x_i[n]\right)\right)\\ + &= B\left(n \mapsto \sum_i a_i A(x_i)[n]\right)\\ + &= n \mapsto \sum_i a_i B(A(x_i))[n]\\ + &= n \mapsto \sum_i a_i C(x_i)[n] +\end{align*} + +\section{Time Invariance} +If $A$ and $B$ are time invariant, i.e. for all signals $x$ and integers $k$, +\begin{align*} + A(n \mapsto x[n - k]) &= n \mapsto A(x)[n - k]\\ + B(n \mapsto x[n - k]) &= n \mapsto B(x)[n - k]\\ +\end{align*} +then $C$ is also time invariant +\begin{align*} + C(n \mapsto x[n - k]) + &= B(A(n \mapsto x[n - k]))\\ + &= B(n \mapsto A(x)[n - k])\\ + &= n \mapsto B(A(x))[n - k]\\ + &= n \mapsto C(x)[n - k] +\end{align*} + +\section{LTI Ordering} +If $A$ and $B$ are linear and time-invariant, there exists signals $g$ and $h$ +such that for all signals $x$, $A = x \mapsto x * g$ and $B = x \mapsto x * h$, +thus \[B(A(x)) = B(x * g) = x * g * h = x * h * g = A(x * h) = A(B(x))\] +or interchanging $A$ and $B$ order does not change $C$. + +\section{Causality} +If $A$ and $B$ are causal, i.e. for all signals $x$, $y$ and integers $k$, +\begin{multline*} + x[n] = y[n]\quad\forall n < k + \Longrightarrow\begin{cases} + A(x)[n] = A(y)[n]\quad\forall n < k\\ + B(x)[n] = B(y)[n]\quad\forall n < k + \end{cases}\\ + \Longrightarrow B(A(x))[n] = B(A(y))[n]\quad\forall n < k + \iff C(x)[n] = C(y)[n]\quad\forall n < k +\end{multline*} +then $C$ is also causal. + +\section{BIBO Stability} +If $A$ and $B$ are stable, i.e. there exists a signal $x$ +and scalars $a$, $b$ that +\begin{align*} + |x[n]| < a\quad\forall n \in \mathbb Z + &\Longrightarrow |A(x)[n]| < b\quad\forall n \in \mathbb Z\\ + |x[n]| < a\quad\forall n \in \mathbb Z + &\Longrightarrow |B(x)[n]| < b\quad\forall n \in \mathbb Z +\end{align*} +then $C$ is also stable, i.e. there exists a signal $x$ +and scalars $a$, $b$, $c$ that +\begin{align*} + |x[n]| < a\;\forall n \in \mathbb Z + &\Longrightarrow |A(x)[n]| < b\;\forall n \in \mathbb Z\\ + &\Longrightarrow |B(A(x))[n]| < c\;\forall n \in \mathbb Z + \iff |C(x)[n]| < c\;\forall n \in \mathbb Z +\end{align*} +\end{document} diff --git a/usth/ICT2.9/dsp.ipynb b/usth/ICT2.9/dsp.ipynb deleted file mode 100644 index 904ac67..0000000 --- a/usth/ICT2.9/dsp.ipynb +++ /dev/null @@ -1,981 +0,0 @@ -{ - "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 -} diff --git a/usth/ICT2.9/filters.ipynb b/usth/ICT2.9/filters.ipynb deleted file mode 100644 index f4ebdcb..0000000 --- a/usth/ICT2.9/filters.ipynb +++ /dev/null @@ -1,613 +0,0 @@ -{ - "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 -} diff --git a/usth/ICT2.9/output_bartlett.wav b/usth/ICT2.9/output_bartlett.wav deleted file mode 100644 index a3ae1ff..0000000 Binary files a/usth/ICT2.9/output_bartlett.wav and /dev/null differ diff --git a/usth/ICT2.9/output_bessel.wav b/usth/ICT2.9/output_bessel.wav deleted file mode 100644 index 27f8e33..0000000 Binary files a/usth/ICT2.9/output_bessel.wav and /dev/null differ diff --git a/usth/ICT2.9/output_blackman.wav b/usth/ICT2.9/output_blackman.wav deleted file mode 100644 index 43cc293..0000000 Binary files a/usth/ICT2.9/output_blackman.wav and /dev/null differ diff --git a/usth/ICT2.9/output_butterworth.wav b/usth/ICT2.9/output_butterworth.wav deleted file mode 100644 index 7880779..0000000 Binary files a/usth/ICT2.9/output_butterworth.wav and /dev/null differ diff --git a/usth/ICT2.9/output_chebyshev.wav b/usth/ICT2.9/output_chebyshev.wav deleted file mode 100644 index db55a6e..0000000 Binary files a/usth/ICT2.9/output_chebyshev.wav and /dev/null differ diff --git a/usth/ICT2.9/output_hamming.wav b/usth/ICT2.9/output_hamming.wav deleted file mode 100644 index 0596d61..0000000 Binary files a/usth/ICT2.9/output_hamming.wav and /dev/null differ diff --git a/usth/ICT2.9/output_sinc.wav b/usth/ICT2.9/output_sinc.wav deleted file mode 100644 index ee7a28e..0000000 Binary files a/usth/ICT2.9/output_sinc.wav and /dev/null differ diff --git a/usth/ICT2.9/practical/blue-ocean-floor.wav b/usth/ICT2.9/practical/blue-ocean-floor.wav new file mode 100644 index 0000000..4208cd0 Binary files /dev/null and b/usth/ICT2.9/practical/blue-ocean-floor.wav differ diff --git a/usth/ICT2.9/practical/dsp.ipynb b/usth/ICT2.9/practical/dsp.ipynb new file mode 100644 index 0000000..904ac67 --- /dev/null +++ b/usth/ICT2.9/practical/dsp.ipynb @@ -0,0 +1,981 @@ +{ + "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 +} diff --git a/usth/ICT2.9/practical/filters.ipynb b/usth/ICT2.9/practical/filters.ipynb new file mode 100644 index 0000000..f4ebdcb --- /dev/null +++ b/usth/ICT2.9/practical/filters.ipynb @@ -0,0 +1,613 @@ +{ + "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 +} diff --git a/usth/ICT2.9/practical/output_bartlett.wav b/usth/ICT2.9/practical/output_bartlett.wav new file mode 100644 index 0000000..a3ae1ff Binary files /dev/null and b/usth/ICT2.9/practical/output_bartlett.wav differ diff --git a/usth/ICT2.9/practical/output_bessel.wav b/usth/ICT2.9/practical/output_bessel.wav new file mode 100644 index 0000000..27f8e33 Binary files /dev/null and b/usth/ICT2.9/practical/output_bessel.wav differ diff --git a/usth/ICT2.9/practical/output_blackman.wav b/usth/ICT2.9/practical/output_blackman.wav new file mode 100644 index 0000000..43cc293 Binary files /dev/null and b/usth/ICT2.9/practical/output_blackman.wav differ diff --git a/usth/ICT2.9/practical/output_butterworth.wav b/usth/ICT2.9/practical/output_butterworth.wav new file mode 100644 index 0000000..7880779 Binary files /dev/null and b/usth/ICT2.9/practical/output_butterworth.wav differ diff --git a/usth/ICT2.9/practical/output_chebyshev.wav b/usth/ICT2.9/practical/output_chebyshev.wav new file mode 100644 index 0000000..db55a6e Binary files /dev/null and b/usth/ICT2.9/practical/output_chebyshev.wav differ diff --git a/usth/ICT2.9/practical/output_hamming.wav b/usth/ICT2.9/practical/output_hamming.wav new file mode 100644 index 0000000..0596d61 Binary files /dev/null and b/usth/ICT2.9/practical/output_hamming.wav differ diff --git a/usth/ICT2.9/practical/output_sinc.wav b/usth/ICT2.9/practical/output_sinc.wav new file mode 100644 index 0000000..ee7a28e Binary files /dev/null and b/usth/ICT2.9/practical/output_sinc.wav differ diff --git a/usth/ICT2.9/speakerid b/usth/ICT2.9/speakerid new file mode 160000 index 0000000..b79672f --- /dev/null +++ b/usth/ICT2.9/speakerid @@ -0,0 +1 @@ +Subproject commit b79672fba4a4419bc301184c1f5c1eac21572780 -- cgit 1.4.1