7.5 Resampling

7.5 Resampling#

It is often useful to change the sample rate of audio after it has already been sampled, perhaps to shrink a file for transmission, or to combine two recordings made at different rates. This operation is called resampling. Its type signature maps one sample vector to another, generally of a different length:

\[\mathbf{x} = [x[0], \ldots, x[N-1]] \;\longrightarrow\; \mathbf{y} = [y[0], \ldots, y[M-1]].\]

Changing the sample rate#

Suppose we want to convert audio from a rate \(f_s^1\) to a new rate \(f_s^2\), keeping its duration in seconds unchanged. Since duration is \(N / f_s^1 = M / f_s^2\), the new length must be

\[M = N \cdot \frac{f_s^2}{f_s^1}.\]

To fill in the new samples, we read the original signal at the corresponding fractional positions. The \(m\)-th output sample comes from position \(p = m \cdot f_s^1 / f_s^2\) in the original, which generally falls between two original samples:

\[y[m] = \text{Interpolate}\!\left(\mathbf{x}, \; p = m \cdot \frac{f_s^1}{f_s^2}\right).\]

We already met this idea in Chapter 3, where wavetable synthesis read a table at fractional positions. The simplest choice is linear interpolation between the two neighboring samples:

\[y[m] = (1 - \alpha)\, x[\lfloor p \rfloor] + \alpha \, x[\lfloor p \rfloor + 1], \qquad \alpha = p - \lfloor p \rfloor.\]

A standalone linear resampler, along with the aliasing and quantization helpers from this chapter, is in code/sampling.py.

The widget below resamples a 1 Hz sine from \(f_s^1 = 8\) Hz to a new rate \(f_s^2\) of your choosing. Each red cross is read from a fractional position between two blue samples, on the straight line that joins them.

# hide
# no-output
from IPython.utils.capture import capture_output
with capture_output():
    %pip install -q plotly anywidget

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import ipywidgets as widgets
import icm_plotly
from icm_plotly import RED, BLUE, GOLD, IRON, TEAL, STEEL

Drag the new rate \(f_s^2\). The blue dots are the original samples at \(f_s^1 = 8\) Hz, and the red crosses are the resampled points. Each cross sits on the dashed line joining its two blue neighbors, which is exactly the value linear interpolation reads. The gray curve is the true signal.

# hide
# autorun
FS1 = 8.0                           # the original rate, as in the figure above
DUR = 2.0                           # two full cycles of a 1 Hz sine
N = int(DUR * FS1)
X = np.sin(2 * np.pi * np.arange(N) / FS1)
T_SMOOTH = np.linspace(0.0, DUR, 600)
FS2_0 = 12.0                        # starting parameter

def resample(fs2, N=N, X=X, FS1=FS1):
    # the chapter's linear interpolation; the two cycles repeat, so the
    # neighbor after the last sample is the first one again
    M = int(round(N * fs2 / FS1))
    p = np.arange(M) * FS1 / fs2
    i = np.floor(p).astype(int)
    a = p - i
    y = (1 - a) * X[i] + a * X[(i + 1) % N]
    return np.arange(M) / fs2, y

def figure():
    fig = go.Figure()
    tm, ym = resample(FS2_0)
    fig.add_scatter(x=T_SMOOTH, y=np.sin(2 * np.pi * T_SMOOTH), mode="lines",
                    line=dict(color=STEEL, width=1.8))
    fig.add_scatter(x=np.arange(N + 1) / FS1, y=np.append(X, X[0]), mode="lines",
                    line=dict(color=BLUE, width=1.2, dash="dash"))
    fig.add_scatter(x=np.arange(N) / FS1, y=X, mode="markers",
                    marker=dict(color=BLUE, size=9))
    fig.add_scatter(x=tm, y=ym, mode="markers",
                    marker=dict(color=RED, size=10, symbol="x-thin",
                                line=dict(color=RED, width=2.4)))
    fig.update_xaxes(range=[-0.05, DUR + 0.05], title_text="Time (s)", fixedrange=True)
    fig.update_yaxes(range=[-1.2, 1.2], title_text="Amplitude", fixedrange=True)
    return fig

def controls(fig):
    fs2 = widgets.FloatSlider(description="New rate $f_s^2$ (Hz)", min=3, max=24,
                              value=FS2_0, step=1, readout_format=".0f")
    readout = widgets.HTML()

    # the defaults snapshot the helpers; the page's notebooks share one kernel
    def update(fs2, resample=resample, N=N, FS1=FS1, readout=readout):
        tm, ym = resample(fs2)
        with fig.batch_update():
            fig.data[3].x, fig.data[3].y = tm, ym
        readout.value = (f"<span style='font-size:0.9em'><i>M</i> = <i>N</i> &middot; "
                         f"<i>f</i><sub>s</sub><sup>2</sup> / <i>f</i><sub>s</sub><sup>1</sup> = "
                         f"{N} &middot; {fs2:.0f} / {FS1:.0f} = {len(tm)} samples"
                         f"</span>")

    widgets.interactive_output(update, {"fs2": fs2})
    return widgets.VBox([fs2, readout])

icm_plotly.show(figure, controls)

There is one critical caveat. When we lower the sample rate (\(f_s^2 < f_s^1\)), we shrink the Nyquist frequency, and any content above the new Nyquist \(f_s^2/2\) will alias, just as in the analog case. So before downsampling, we must first filter out everything above \(f_s^2/2\), an anti-aliasing step we will be equipped to implement after studying filters in Chapter 9. In practice, high-quality resamplers combine this filtering with a more sophisticated interpolation than the linear scheme above. Pyquist’s Audio.resample handles both:

import pyquist as pq

audio = pq.Audio.from_file("drums.wav")   # 44.1 kHz
half = audio.resample(22050)              # bandlimited, anti-aliased
low = audio.resample(8000)

Listen to a recording resampled to progressively lower rates. As the sample rate drops, the Nyquist frequency falls below the signal’s high-frequency content, and that content is (properly) removed, so the sound grows progressively duller:

Original (44.1 kHz)

Resampled to 22.05 kHz

Resampled to 8 kHz

A recording resampled to lower rates. The 8 kHz version has a Nyquist frequency of only 4 kHz, so everything above that is gone and the sound is noticeably muffled. 666866 by MrJmix, License: Attribution 4.0.

7.5 Resampling

Contents

7.5 Resampling#

Changing the sample rate#