Upsampling and interpolation

To raise a sample rate, I start by putting zeros between the samples. Here a tone, $\cos(0.40\pi n)$, gets $L-1$ zeros after each sample, for $L$ from 1 to 4. Watch its spectrum: the tone’s line slides down, and new lines appear in the space it left.

Insert zeros: the spectrum squeezes and images appear

In Downsampling and decimation (22.1) we lowered the sample rate by keeping every $M$-th sample. Now I want to go the other way and raise it. Shifting, reversing and scaling time (2.1) already met the problem. Slowing down by 2 asks for samples halfway between the measured ones, at instants that were never measured.

So let’s leave the gaps empty first and fill them later. To raise the rate by a whole number $L$, I put $L-1$ zeros after each sample. This is upsampling by $L$, also called zero insertion or zero-stuffing. I write the result $x_u[n]$:

Nothing is lost. Every original sample is still there, now $L$ places from the next one.

What do the zeros do to a tone? Take $x[n]=\cos(0.40\pi n)$ and $L=2$. At even $n$, $x_u[n]=x[n/2]=\cos(0.40\pi\cdot n/2)=\cos(0.20\pi n)$, and at odd $n$ it is 0. The factor $\tfrac12\big(1+(-1)^n\big)$ is 1 at even $n$ and 0 at odd $n$, so

The second line takes one step. $(-1)^n$ is $e^{j\pi n}$, and by Properties of the DTFT (12.3), multiplying by it slides every line of a spectrum by $\pi$. The tone’s lines at $\pm0.20\pi$ move to $1.20\pi$ and $0.80\pi$. On the circle of Frequency in discrete time (12.1), $1.20\pi$ is $-0.80\pi$, so the product is $\cos(0.80\pi n)$.

So the zeros did two things. The tone moved to half its frequency, $0.20\pi$. And a second tone appeared at $\pi-0.20\pi=0.80\pi$, which was not in $x$ at all. I call it an image.

The tone and its image each have amplitude $\tfrac12$.

Think of a picket fence with a blank picket painted in after every real one. The old rhythm is still there, at half the pace, and a new, false rhythm sits beside it. The picture at the top of the page does this to the tone for $L$ from 1 to 4.

Watch two things. Every line’s height falls as $1/L$: 1, 0.500, 0.333, 0.250. And the tone slides down to $0.40\pi/L$, while new image lines grow in the space the tone left. After the clip, the slider goes on to $L=5$ and 6.

Squeezed by L

Why does every frequency divide by $L$? Write the DTFT of $x_u$ (The DTFT, 12.2). Only the terms with $n=Lm$ are not zero, and there $x_u[Lm]=x[m]$:

So the new spectrum is the old one, read at $\Omega L$. What sat at $\Omega$ now sits at $\Omega/L$: the spectrum is squeezed by $L$. The old spectrum repeats every $2\pi$ (12.1), so the new one repeats every $2\pi/L$. That fits $L$ copies into one turn, and the extra copies are the images.

Now follow a tone at $\Omega_1$. The old spectrum has its line at $\Omega_1$ and at every copy, $\Omega_1+2\pi k$. The new one has a line wherever $\Omega L$ hits one of those, so at

Larger $k$ only repeat these a whole turn further on. Fold each one into 0 to $\pi$, as in 22.1. For $\Omega_1=0.40\pi$ and $L=4$ that gives $0.10\pi$, $0.60\pi$, $1.10\pi$ and $1.60\pi$, which fold to $0.10\pi$, $0.60\pi$, $0.90\pi$ and $0.40\pi$.

Why is each line $1/L$ high? For $L=2$ the factor $\tfrac12\big(1+(-1)^n\big)$ did it. For any $L$ the same job is done by

When $L$ divides $n$, every arrow is 1 and the factor is 1. Otherwise the $L$ arrows are spread evenly round the circle and cancel, like the rows of The DFT as a matrix (13.5). So for the tone, $x_u[n]$ is $\cos(\Omega_1n/L)$ times this factor.

Each arrow $e^{j2\pi kn/L}$ slides the squeezed tone by $2\pi k/L$ (12.3) and carries the weight $1/L$. That makes $L$ lines, each a cosine of amplitude $1/L$.

Remove the images, and multiply by L

The zeros moved the tone to the right place, $\Omega_1/L$, and added images. So the repair is a low-pass that keeps the squeezed spectrum and removes the rest. The old band, 0 to $\pi$, now fills 0 to $\pi/L$, so the cutoff is $\pi/L$.

The tone also came out $L$ times too small. So the filter gets a gain of $L$ as well. A low-pass at $\pi/L$ with gain $L$ is the interpolation filter.

For $L=4$ the cutoff is $\pi/4=0.25\pi$, and 22.1 already has that filter. It is the 75-tap low-pass, a Kaiser design (Window-method FIR design, 19.1) with $\beta=5.6533$ and cutoff $0.25\pi$. It passes 0 to $0.2\pi$ within $\pm0.00111$, stays at least 60.38 dB down from $0.3\pi$, and delays by 37 samples. Here I multiply its taps by 4.

Will it do? The tone sits at $0.10\pi$, in the pass band. The images sit at $0.40\pi$, $0.60\pi$ and $0.90\pi$, all past $0.3\pi$, in the stop band.

In fact it serves any $x$ with nothing above $0.8\pi$. That content squeezes to at most $0.2\pi$, and its nearest image starts at $(2\pi-0.8\pi)/4=0.3\pi$.

Think of a flip-book with a blank page after every drawing. The filter draws the in-between pictures.

The second panel gives levels in dB re 1: a line of amplitude $A$ sits at $20\log_{10}A$. So the four lines of amplitude 0.250 sit at −12.0 dB, and the filter’s gain of 4 is +12.04 dB.

Watch the images sink and the gaps fill. The images fell from −12.0 dB to −72.4, −74.3 and −79.6 dB. The tone, which was 0.250, came out at 1.000. And in time, the output samples lie on the sine to within 0.0006.

Why the gain is L, and why the gaps fill with the right values

The gain is $L$ because zero insertion divided the tone’s amplitude by $L$. Only one sample in $L$ carries the signal, and the filter has to bring the level back.

Why does the output land on the sine? The filter’s taps are a windowed sinc, centred on tap 37. Taps 4, 8, 12, … away from the centre are 0, because the ideal low-pass at $0.25\pi$ has $\sin(0.25\pi m)=0$ when $m$ is a multiple of 4. The centre tap of $4h$ is 0.99969.

Take an output at the place of an original sample. The input’s other non-zero samples are 4, 8, 12, … places away, so they meet only zero taps. The output there is 0.99969 times the original sample, so the original samples are kept almost exactly.

In between, every non-zero tap adds in a scaled pulse from each nearby original sample. The ideal taps times 4 are $\sin(\pi m/4)/(\pi m/4)=\mathrm{sinc}(m/4)$: a sinc pulse that is 0 at every old sample instant. That is the sum of Reconstruction (10.3), a sinc pulse at every sample, added, now done in discrete time and cut to 75 taps by the window.

The same 75 taps decimated by 4 in 22.1. Here, times 4, they interpolate by 4. In SciPy, upfirdn(4*h, x, 4) inserts the zeros and filters in one call. resample_poly(x, 4, 1) does the same with a Kaiser filter of its own, and it removes the delay for you.

Notice that three inputs in every four are zeros. For each output, only 18 or 19 of the 75 taps meet a non-zero input, and the rest multiply zeros. Polyphase structures (22.4) skips those multiplications.

Worked example

Let’s redo the page’s numbers, by hand where we can and with SciPy where we can’t.

1. The images of $0.40\pi$. The lines sit at $(0.40\pi+2\pi k)/L$, folded into 0 to $\pi$. For $L=3$: $k=0$ gives $0.40\pi/3=0.1333\pi$, $k=1$ gives $2.40\pi/3=0.8000\pi$, and $k=2$ gives $4.40\pi/3=1.4667\pi$, which folds to $2\pi-1.4667\pi=0.5333\pi$.

The same rule gives these lines, each of amplitude $1/L$:

2. After the filter. The tone’s amplitude becomes $0.25\cdot4\,\lvert H(e^{j0.10\pi})\rvert=1.0000$. Each image becomes $0.25\cdot4\,\lvert H\rvert=\lvert H\rvert$ at its frequency: −72.4497 dB at $0.40\pi$, −74.3090 dB at $0.60\pi$ and −79.6460 dB at $0.90\pi$.

3. The taps against the sinc. Next to the centre, $4h$ is 0.89836, 0.63167 and 0.29499. The ideal $\mathrm{sinc}(m/4)$ for $m=1,2,3$ is 0.90032, 0.63662 and 0.30011. The window and SciPy’s scaling pull each tap down slightly.

4. One in-between value. Take an output one place after an original sample where the sine is at its peak. The sine there is $\cos(0.10\pi)=0.95106$. The filter gives 0.95096, which is 0.0001 low. The original sample beside it comes out at 0.99969 in place of 1.

Where you’ll meet this

Digital-to-analog converters upsample by 4 to 256 before the analog filter. The images then sit far above the audio band, so a gentle analog filter removes them. That is the easier filter of Anti-aliasing and practical converters (10.4), on the output side, and sigma-delta converters (Oversampling and noise shaping, 11.3) take it furthest.

Audio effects that would otherwise alias, such as distortion, run at a raised rate (Audio effects, 29.2). Zooming an image is interpolation along its rows and then its columns.

A ratio that is not a whole number is Resampling by any factor (22.3). In Filter banks (23.1), each band is upsampled and filtered like this, and the bands are added to rebuild the signal.

Reference card