Resampling by any factor

Here a 10 kHz tone recorded at 32 kHz is converted to 48 kHz: up by 3, through one low-pass, then down by 2. Watch the 42 kHz image: the filter sinks it before halving the rate can fold it to 6 kHz.

Up by 3, filter once, down by 2

Say I have a recording made at 32 kHz, and the rest of my project runs at 48 kHz. Downsampling and decimation (22.1) divides a rate by a whole number, and Upsampling and interpolation (22.2) multiplies it by one. But 48/32 is 3/2, and neither step gets there alone.

Think of a recipe for 2 people that you want to cook for 3. You can scale it up to 6 portions and then take half. Rates work the same way: up by $L=3$ to 96 kHz, then down by $M=2$ to 48 kHz. Changing the rate by a ratio $L/M$ of whole numbers like this is rational resampling, and 96 kHz is the intermediate rate.

Each step brings its own problem. Putting zeros between the samples makes images (22.2): copies of what the input holds below its Nyquist frequency, 16 kHz, appear above it. Halving the rate folds everything above the new Nyquist frequency, 24 kHz, back down (22.1).

Both cures are a low-pass at 96 kHz, so one low-pass can do both jobs. It must stop everything above 16 kHz and everything above 24 kHz, so its cutoff is the smaller one, 16 kHz. Its gain is 3, the $L$ of 22.2, because the zeros divided every amplitude by 3.

My test signal is a 10 kHz tone at 32 kHz. After the zeros, 96 kHz holds three lines: the tone at 10 kHz, and images at 32 − 10 = 22 kHz and 32 + 10 = 42 kHz. Each has a third of the tone’s amplitude, which is −9.5 dB.

The filter is a Kaiser low-pass from Window-method FIR design (19.1), designed for 60 dB. Its pass band runs to 14 kHz and its stop band starts at 18 kHz, and kaiserord asks for 89 taps. Kaiser’s formula is an estimate again: the stop band reaches −59.5 dB, a little short of 60.

The picture at the top of the page follows the tone through the three stages, one panel for each rate. As the filter’s gain reaches the tone, the tone rises to 0 dB, because the gain of 3 undoes the zeros’ third. One filter serves both steps.

Notice the 42 kHz image. If the filter had let it through, halving the rate would have folded it to 6 kHz. That is in the middle of the audio band, a tone that was never in the input.

One filter for any ratio

Let’s say the same thing in rad/sample. At the intermediate rate of 96 kHz, 16 kHz is $\Omega=2\pi\cdot16/96=\pi/3$, which is $\pi/L$. And 24 kHz is $\pi/2$, which is $\pi/M$.

That holds for any ratio. After the zeros the images start above $\pi/L$, and keeping every $M$-th sample folds whatever lies above $\pi/M$. So the one filter needs

When the rate goes up, $L>M$, the images set the cutoff. When it goes down, the folding sets it, and the filter is 22.1’s decimation filter.

The three steps fit in one sum. Filter the zero-stuffed signal $x_u$ with $h$. It is zero except at the positions $Lk$, where it holds $x[k]$, so the filtered value at position $m$ is $\sum_k x[k]\,h[m-Lk]$. Keeping every $M$-th value means $m=Mn$:

Look at what this sum leaves out. It multiplies no zeros, and it computes no value that would be thrown away. SciPy’s upfirdn(h, x, L, M) computes this sum, and resample_poly(x, L, M) designs a Kaiser low-pass for you and then does the same.

Write the ratio in lowest terms. 6/4 is also 3/2, but it would need an intermediate rate of 192 kHz and a filter twice as long for the same job.

A good filter and a poor one

Does the filter matter much? Let’s compare two. The poor one is linear interpolation: join the input samples with straight lines and read the lines at the new times. It is the “straight lines” pulse of Reconstruction (10.3).

At 96 kHz each gap between two input samples gets two new samples, a third and two thirds of the way across. On the line they are $\tfrac23x[0]+\tfrac13x[1]$ and $\tfrac13x[0]+\tfrac23x[1]$. That is the zero-stuffed signal filtered by the triangle (1, 2, 3, 2, 1)/3. Its middle tap keeps each input sample, and its other taps reach into the gaps.

What is its gain? The triangle is the box (1, 1, 1) convolved with itself, divided by 3. Centred on 0, the box has gain $e^{j\Omega}+1+e^{-j\Omega}=1+2\cos\Omega$. Convolving multiplies gains (Properties of the DTFT, 12.3), so the triangle, centred on 0, has

At $\Omega=0$ it is 3, the gain $L$. It is 0 where $\cos\Omega=-1/2$, at $\Omega=2\pi/3$, which is 32 kHz: the image of 0 Hz. But our images sit at 22 and 42 kHz, on the slopes either side of that zero, and there the gain is far from 0.

The good filter is the 89-tap Kaiser low-pass of the first picture. The next picture runs the same conversion with each filter and shows the 48 kHz output. Its levels are measured from the output tone, so the tone always sits at 0 dB.

Watch the image and the alias as the filter changes: the steps are the same, and only the filter changed. Choose each filter and press “Hear it”. With linear interpolation the 6 kHz alias is a second, lower whistle beside the tone. The 22 kHz image is there too, but it is above what most adults can hear.

Notice one more cost of linear interpolation that the panel hides. It turns the tone itself down, to 0.7435 of its amplitude, or −2.57 dB. That is the droop of 10.3’s straight lines. The Kaiser filter passes the tone at 1.0000.

Four weights that slide

Rational resampling needs a fixed ratio. Sometimes there is none. A USB microphone and a sound card each run on their own clock. Both say 48 kHz, but the true ratio of their rates is a little off 1 and drifts as the clocks warm up.

Then each output sample falls somewhere between two input samples, and the place changes from one output to the next. I write $\mu$ for that place: the position between two input samples, with $0\le\mu<1$. At $\mu=0$ the output sits on $x[0]$; at $\mu=0.5$ it is halfway to $x[1]$. This $\mu$ is local to this page.

Linear interpolation reads the straight line, $(1-\mu)\,x[0]+\mu\,x[1]$. The fractional delay of Special FIR filters (19.4) does much better, with 16 taps of a shifted sinc. But every one of those taps changes whenever $\mu$ changes.

A middle way is Lagrange interpolation. Pass the one cubic through the four nearest samples, $x[-1]$, $x[0]$, $x[1]$ and $x[2]$, and read it at $\mu$. Written as weights on the samples, the output $y[n]$ that falls at $\mu$ is

Each weight is a cubic that is 1 at its own sample and 0 at the other three. Put $\mu=0$, and only the weight of $x[0]$ is left, and it is 1. So the curve passes through every sample, as 10.3’s sinc sum did.

The two middle weights do most of the work. The outer two dip a little below 0, to −0.064 at most. At every $\mu$ the four weights add to 1, so a constant signal comes through unchanged.

The Farrow structure

Each weight is a cubic in $\mu$, so the output is a cubic in $\mu$ too. Multiply out the four weights and collect the powers of $\mu$. The output is a number times $\mu^3$, plus a number times $\mu^2$, plus a number times $\mu$, plus $x[0]$. Each number is a fixed mix of the four samples:

The mixes are fixed, so a converter computes them once for each gap. For each new $\mu$ it evaluates the cubic the cheap way.

Start with the $\mu^3$ number and multiply by $\mu$. Add the $\mu^2$ number and multiply by $\mu$. Add the $\mu$ number, multiply by $\mu$ and add $x[0]$. That is three multiplies by $\mu$ for each output.

Fixed filters on the samples, followed by a polynomial in $\mu$, make a Farrow structure. It suits a $\mu$ that changes at every output, because nothing has to be redesigned.

Now back to the two clocks. The converter keeps a running position, counted in input samples. For each output it moves the position on by the measured ratio of the input rate to the output rate.

The whole part of the position says which four samples to use, and the fraction is $\mu$. This is asynchronous rate conversion. The ratio is measured again and again, so the converter follows the drift.

The interpolator can be 19.4’s shifted sinc or Lagrange in a Farrow structure. Lagrange is accurate for slow signals. At $\mu=0.5$ it passes $0.1\pi$ with gain 0.9998, but $0.5\pi$ with only 0.884. So a converter can first raise the rate by a fixed factor with a good low-pass, and then use Lagrange between those closer samples.

44.1 kHz and 48 kHz

CDs use 44.1 kHz, while video and many sound cards use 48 kHz. The largest whole number that divides both 44 100 and 48 000 is 300. So 48 000/44 100 is 160/147, and the converter goes up by 160 and down by 147.

The intermediate rate is $44.1\cdot160=7056$ kHz, which is 7.056 MHz. 159 of every 160 samples there are zeros. The filter’s cutoff is $\min(\pi/160,\pi/147)=\pi/160$, which is 22.05 kHz, the input’s Nyquist frequency.

SciPy’s resample_poly(x, 160, 147) uses a 3201-tap filter. In Real-time processing (21.4) we sized jobs against a 100 MHz processor that does one multiply each clock cycle. Here is this job, done three ways:

Only the last route fits, and it is the sum $y[n]=\sum_k x[k]\,h[Mn-Lk]$ from above. Arranging a filter so that it works this way is the subject of Polyphase structures (22.4).

Worked example

Let’s redo the page’s numbers, by hand where we can.

1. The ratio for 44.1 to 48 kHz. Factor both rates: $44\,100=2^2\cdot3^2\cdot5^2\cdot7^2$ and $48\,000=2^7\cdot3\cdot5^3$. The common part is $2^2\cdot3\cdot5^2=300$. Then 48 000/300 = 160 and 44 100/300 = 147. Since $147=3\cdot7^2$ and $160=2^5\cdot5$ share no factor, $L/M=160/147$ is in lowest terms.

2. Linear interpolation’s levels. At 96 kHz, $\Omega=2\pi f/96$ kHz, and the gain is $(1+2\cos\Omega)^2/3$:

Measured from the tone, the 22 kHz image is $20\log_{10}(0.5301/2.2304)=-12.48$ dB. The 42 kHz image, which becomes the 6 kHz alias, is $20\log_{10}(0.2396/2.2304)=-19.38$ dB. The tone keeps $2.2304/3=0.7435$ of its amplitude.

3. The two filters side by side. The Kaiser values come from freqz of the 89 taps.

4. Lagrange weights at μ = 0.25. The weight of $x[-1]$ is $-\tfrac16(0.25)(-0.75)(-1.75)=-0.0547$. The weight of $x[0]$ is $\tfrac12(1.25)(-0.75)(-1.75)=0.8203$. The weight of $x[1]$ is $-\tfrac12(1.25)(0.25)(-1.75)=0.2734$. The weight of $x[2]$ is $\tfrac16(1.25)(0.25)(-0.75)=-0.0391$.

As fractions they are −7/128, 105/128, 35/128 and −5/128, which add to 128/128 = 1. At μ = 0.5 they are −1/16, 9/16, 9/16 and −1/16, which also add to 1.

5. Lagrange on a slow wave. Take $x[n]=\cos(0.1\pi n)$. The samples at −1, 0, 1 and 2 are 0.95106, 1, 0.95106 and 0.80902. At $\mu=0.5$ the output is $0.5625\,(1+0.95106)-0.0625\,(0.95106+0.80902)=0.98746$.

The true value is $\cos(0.05\pi)=0.98769$, so the error is −0.00022. Linear interpolation gives $(1+0.95106)/2=0.97553$, an error of −0.0122, about 54 times larger.

Where you’ll meet this

Audio editors and the sound mixers of operating systems convert between 44.1 and 48 kHz, for example when a CD-rate file plays on a 48 kHz sound card. SciPy’s resample_poly is a rational resampler with a polyphase filter.

Video frame-rate conversion uses the same ideas along time, with frames in place of samples. Software radio uses them to move between the converter’s rate and the rate a receiver needs. Fractional delays come back in Audio effects (29.2), where delay lines that move smoothly need them.

Reference card