-
Notifications
You must be signed in to change notification settings - Fork 23
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
2 changed files
with
434 additions
and
635 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,328 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# FIR Resampling\n", | ||
| "<div align=\"right\"><a href=\"https://people.epfl.ch/paolo.prandoni\">Paolo Prandoni</a>, <a href=\"https://www.epfl.ch/labs/lcav/\">LCAV, EPFL</a></div>\n", | ||
| "\n", | ||
| "In this notebook we will implement a simple fractional resampler that uses local Lagrange interpolation. The resampler can be used to perform any rational sampling rate change (but beware of aliasing when downsampling!)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import numpy as np\n", | ||
| "import IPython\n", | ||
| "from scipy.io import wavfile" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Local Lagrange interpolation\n", | ||
| "\n", | ||
| "To perform resampling we will need to compute subsample values and, in order to compute subsample values without resorting to a full interpolation of the input signal, the idea is to fit a local polynomial interpolator arount the reference index $n$ and compute the value _of the polynomial_ at $n+\\tau$, where $|\\tau| < 1/2$. In the following figure, for instance, a quadratic interpolator is used:\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "To fit the polynomial (whose degree, for reasons of symmetry, should be even) we can use Lagrange polynomials. In this case the approximation is simply: \n", | ||
| "\n", | ||
| "$$\n", | ||
| " x_L(n; \\tau) = \\sum_{k=-N}^{N} x[n + k] L_k^{(N)}(\\tau)\n", | ||
| "$$\n", | ||
| "\n", | ||
| "where $L_k^{(N)}(t)$ is the $k$-th Lagrange polynomial of order $2N$." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Rational sampling rate change\n", | ||
| "\n", | ||
| "Given a nominal input sampling rate $F_i$ and an output sampling rate $F_o$, a fractional resampler generates $N_o$ output samples for every $N_i$ input samples where\n", | ||
| "\n", | ||
| "$$\n", | ||
| " \\frac{N_o}{N_i} = \\frac{F_o}{F_i}\n", | ||
| "$$\n", | ||
| "with $N_i, N_o$ coprime for maximum efficiency. So the first thing we need is a simple function that simplifies the ratio of sampling frequencies to its lowest terms. For this we will use Euclid's algorithm:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def simplify(A, B):\n", | ||
| " # Euclid's GCD algorithm\n", | ||
| " a = A\n", | ||
| " b = B\n", | ||
| " while a != b:\n", | ||
| " if a > b:\n", | ||
| " a = a - b\n", | ||
| " else:\n", | ||
| " b = b - a\n", | ||
| " return A // a, B // b" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We can test the function on the usual CD to DVD sampling rate change and indeed we obtain the familiar 160/147 ratio:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": { | ||
| "scrolled": true | ||
| }, | ||
| "outputs": [], | ||
| "source": [ | ||
| "print(simplify(48000, 44100))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Fractional resampling" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "A local fractional resampler works as follows:\n", | ||
| " * for each output index $m$, find the closest input index $n$, called the \"anchor\", so that the distance between input and output instants (in seconds) is less than one half of the input sampling period in magnitude; let $\\tau$ be this distance, normalized by the input sampling period; clearly $|\\tau| < 1/2$;\n", | ||
| " * perform a local Lagrange interpolation of the input around $n$;\n", | ||
| " * compute the value of the interpolation at $n +\\tau$.\n", | ||
| " \n", | ||
| "For a rational sampling rate change, the values of $\\tau$ repeat in a pattern every $N_o$ output point, so that the values of the Lagrange interpolation coefficients can be precomputed; for every block of $N_o$ output points we will use $N_i$ input points. " | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "The method is best understood graphically: in the figures below, the top line shows the samples $x[n]$ in the original sequence, while the bottom line shows the resampled values $y[n]$. The red dotted lines indicate the intervals for each _input_ sample where $\\tau$ less than one half in magnitude; the blue arrows show which input sample is used as an anchor point in the computation of an output sample.\n", | ||
| "\n", | ||
| "Let's first consider downsampling, with a ratio $N_o/N_i = 4/5$; the operation generates 4 output samples for every 5 input samples; the anchor points are determined like so:\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "As apparent from the figure, since the output rate is less than the input rate, every once in a while an input sample is skipped." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "In the following figure, the sgnal is upsampled with a ratio $N_o/N_i = 8/5$ generates 8 output samples for every 5 input samples; since the output rate is larger than the input rate, some of the input samples need to be reused.\n", | ||
| "\n", | ||
| "" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Mathematically, the anchor points are determined like so:\n", | ||
| "\n", | ||
| " * the output sample $y[m]$ occurs at time $t_m = m/F_o$\n", | ||
| " * the closest input sample will be $x[n]$ (occurring at time $t_n = n/F_i$) where $n$ is such that $|m/F_o - n/F_i| < 1/2$.\n", | ||
| " \n", | ||
| "The required value for $n$ can be found by setting \n", | ||
| "$$\n", | ||
| " \\left|\\mbox{frac}\\left(\\frac{m}{N_o} - \\frac{n}{N_i}\\right)\\right| < \\frac{1}{2}\n", | ||
| "$$\n", | ||
| "where `frac()` indicates the fractional part of a number. This yields\n", | ||
| "$$\n", | ||
| " n = \\mbox{round}\\left(m\\frac{N_i}{N_o}\\right);\n", | ||
| "$$\n", | ||
| "the fractional distance between $y[m]$ and $x[n]$ is given by the time difference normalized by the input's sampling period, that is\n", | ||
| "$$\n", | ||
| " \\tau = F_i\\left(\\frac{m}{F_o} - \\frac{n}{F_i}\\right) = m\\frac{N_i}{N_o} - \\mbox{round}\\left(m\\frac{N_i}{N_o}\\right).\n", | ||
| "$$\n", | ||
| "Note that $\\tau = 0$ every time $m$ is a multiple of $N_o$, which confirms the repetition pattern every $N_o$ output samples." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "### Setting up the filters\n", | ||
| "\n", | ||
| "In practice, we want the resampling to be driven by the input data, so we set up an array of $N_i$ filter _lists;_ for every new input sample we look at the next element in the list (modulo $N_i$) and produce as many output samples as there are filters in the list element.\n", | ||
| "\n", | ||
| "The following function sets up a set of $N_i$ quadratic interpolation filters for each anchor point:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def setup_filters(output_rate, input_rate):\n", | ||
| " No, Ni = simplify(output_rate, input_rate)\n", | ||
| " filterbank = [[] for _ in range(Ni)]\n", | ||
| " # while output index spans [0, No-1], the input spans [0, Ni-1]\n", | ||
| " for m in range(0, No):\n", | ||
| " anchor = int(m * Ni / No + 0.5) \n", | ||
| " tau = (m * Ni / No) - anchor\n", | ||
| " filterbank[anchor].append([\n", | ||
| " tau * (tau - 1) / 2, \n", | ||
| " (1 - tau) * (1 + tau), \n", | ||
| " tau * (tau + 1) / 2\n", | ||
| " ])\n", | ||
| " return filterbank" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We can test with the examples in the figures above:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "setup_filters(4, 5)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "setup_filters(8, 5)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "### The final interpolator\n", | ||
| "We are now ready to write the full interpolation function:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def resample(output_rate, input_rate, x):\n", | ||
| " No, Ni = simplify(output_rate, input_rate)\n", | ||
| " filterbank = setup_filters(No, Ni)\n", | ||
| " \n", | ||
| " y = np.zeros(No * len(x)) // Ni \n", | ||
| " m = 0\n", | ||
| " for n, x_n in enumerate(x[1:-1]):\n", | ||
| " for fb in filterbank[n % Ni]:\n", | ||
| " y[m] = x[n-1] * fb[0] + x[n] * fb[1] + x[n+1] * fb[2]\n", | ||
| " m += 1\n", | ||
| " return y" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We can now test the resampler on a simple sinusoid; we generate the sinusoid at 44.1 KHz and resample at 48KHz; the pitch should not change:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "x = np.cos(2 * np.pi * 440 / 44100 * np.arange(0, 44100))\n", | ||
| "IPython.display.Audio(x, rate=44100)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "y = resample(12000, 44100, x)\n", | ||
| "IPython.display.Audio(y, rate=12000)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We can now test the resampler on an audio file; note how aliasing appears when we downsample too much:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "Fi, x = wavfile.read('oao.wav')\n", | ||
| "IPython.display.Audio(x, rate=Fi)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "y = resample(48000, Fi, x)\n", | ||
| "IPython.display.Audio(y, rate=48000)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "y = resample(8000, Fi, x)\n", | ||
| "IPython.display.Audio(y, rate=8000)" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.7.4" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 4 | ||
| } |
Large diffs are not rendered by default.
Oops, something went wrong.