Downsampling and decimation

Here are two tones, at 0.10π and 0.35π rad/sample, and I keep only every $M$-th sample, for $M$ = 1 to 4. Watch the second tone’s line: from $M=3$ on it passes π and folds back.

Keep every M-th sample: the spectrum stretches and folds

In Shifting, reversing and scaling time (2.1) we sped a signal up by a whole number $M$: $y[n]=x[Mn]$. You keep every $M$-th sample and throw the rest away. On this page I call that downsampling by $M$, and I write the result as $x_d[n]=x[Mn]$, with d for downsampled.

Why throw samples away? A converter that samples fast, as in Anti-aliasing and practical converters (10.4), hands over far more samples than the band you keep needs. Fewer samples means less to store and less to compute. But dropping samples does something to the frequencies, and that is what this page is about.

Let’s take two tones:

Keep every $M$-th sample by putting $Mn$ in place of $n$:

Read each cosine as $\cos\big((M\Omega)\,n\big)$. Every frequency $\Omega$ is multiplied by $M$. The amplitudes do not change: a cosine of amplitude 1 still has amplitude 1 when you skip some of its samples.

What happens when $M\Omega$ passes $\pi$? Frequency in discrete time (12.1) showed that $\Omega$ and $\Omega+2\pi$ give the same samples. A cosine is even, so $\Omega$ and $-\Omega$ give the same samples too. So a frequency past $\pi$ is the same as one below it: subtract $2\pi$, then drop the sign.

For $M=3$ the second tone goes to $1.05\pi$. Subtract $2\pi$ to get $-0.95\pi$, and drop the sign: $0.95\pi$. This is the fold of Sampling & aliasing (10.1), in rad/sample instead of hertz.

Notice that there is no clock and no analog part here. Aliasing happens inside the computer, from nothing more than dropping samples.

Think of a film projector that skips frames. A turning wheel seems to turn faster, and past some speed it seems to turn backwards. That backwards turn is the fold. The picture at the top of the page drops samples from these two tones for $M$ = 1, 2, 3 and 4, and moves each line to where it lands.

Notice that the 0.35π tone folds as soon as $M=3$, because $3\cdot0.35\pi=1.05\pi$ is past $\pi$. The 0.10π tone moves up without folding, and neither line changes height.

After the clip, use the slider to set $M$ up to 8. Look at $M=6$: the second tone lands at 0.10π, exactly where the first tone was in $x[n]$.

Then try $M=8$. Both tones land on 0.80π, and the two lines stand side by side. They are now the same cosine, since $\cos(2.8\pi n)=\cos(0.8\pi n)$, so $x_d[n]=1.5\cos(0.8\pi n)$. Two tones became one, and nothing in the samples can split them again.

The spectrum after downsampling

For two tones the rule was: multiply each $\Omega$ by $M$, then fold. For a whole spectrum it reads

Here $X$ is the DTFT of $x$ (The DTFT, 12.2), and $X_d$ is the DTFT of $x_d$.

Read the $k=0$ term first. $X(e^{j\Omega/M})$ is $X$ stretched by $M$: whatever $X$ had at $\Omega_0$ now sits at $M\Omega_0$. The other $M-1$ terms are the same stretched spectrum, shifted by $2\pi k$. All $M$ terms are scaled by $1/M$.

Why the shifted copies? A spectrum of samples repeats every $2\pi$ (12.1), but the stretched one repeats only every $2\pi M$. The copies fill the gaps, so the result repeats every $2\pi$ again.

Where a copy overlaps the stretched original, the two add. That is aliasing: the overlapping copies of The sampling theorem (10.2).

And the $1/M$, when the lines in the instrument kept their heights? A tone’s spectrum is made of arrows (12.2). Stretching the $\Omega$ axis by $M$ makes each arrow’s area $M$ times larger, $\delta(\Omega/M)=M\,\delta(\Omega)$. That $M$ cancels the $1/M$, so a sinusoid keeps its amplitude.

When does nothing alias? If $X$ is zero for $\pi/M<\lvert\Omega\rvert\le\pi$, the stretched spectrum ends at $\pi$. The copies, $2\pi$ away, then do not overlap it. So nothing aliases when the signal has no content above $\pi/M$.

Where the formula comes from

Start from the DTFT of $x_d$:

Write $m=Mn$. The sum then runs only over the $m$ that are multiples of $M$, and $\Omega n=\Omega m/M$. To let it run over every $m$, I need a switch $c[m]$ that is 1 when $m$ is a multiple of $M$ and 0 otherwise. The $M$-th roots of unity of Complex numbers for signals (3.3) give one:

When $m$ is a multiple of $M$, every term is 1 and $c[m]=M/M=1$. Otherwise the $M$ arrows are equal steps round the circle, and they add to 0, as in The DFT (13.2). Put $c[m]$ in and swap the two sums:

The inner sum is the DTFT of $x$ at $(\Omega-2\pi k)/M$, and that gives the formula.

Filter first, then keep every M-th

The cure is the one of 10.4: remove what would fold before it can fold. Low-pass first, then keep every $M$-th sample. The pair is called decimation, and the low-pass is the decimation filter.

How strong must the filter be, and where? With $M=4$, let’s keep 0 to 0.2π, which becomes 0 to 0.8π after downsampling.

A component at $\Omega$ lands at $4\Omega$, and past $\pi$ it folds to $2\pi-4\Omega$. That is inside the kept band when $2\pi-4\Omega\le0.8\pi$, which is when $\Omega\ge0.3\pi$.

So content from 0.2π up to 0.3π lands at $2\pi-4\cdot0.3\pi=0.8\pi$ or above, outside the kept band, and the stop band may start at 0.3π. With a kept band from 0 to $\Omega_\text{pass}$, the same steps give

This is 10.4’s $f_\text{stop}=f_s-f_\text{pass}$ in the units of the new rate, which is $2\pi/M$ on the old $\Omega$ axis. The stop band does not have to start at $\pi/M$.

For the filter I use a Kaiser design from Window-method FIR design (19.1). For 60 dB, $\beta=0.1102\,(60-8.7)=5.6533$. The transition runs from 0.2π to 0.3π, so $\Delta\Omega=0.1\pi$, and Kaiser’s length formula asks for 74 taps.

I take 75, so that the delay is a whole number of samples, 37. The cutoff goes in the middle of the transition, 0.25π: in SciPy, firwin(75, 0.25, window=('kaiser', 5.6533)). I call it the 75-tap low-pass. It stays within ±0.00111 of 1 up to 0.2π, and from 0.3π it is at most −60.38 dB.

It is like blurring a photo slightly before you shrink it, so that fine stripes do not turn into false patterns. The instrument puts the 75-tap low-pass in front of the $M=4$ case. Its levels are in dB re 1: $20\log_{10}$ of the amplitude, so amplitude 1 is 0 dB and amplitude 0.5 is −6.0 dB.

Watch the alias at 0.60π: without the filter it is only 6 dB below the wanted tone, and after it, −74.2 dB, below the 60 dB line.

Notice where the filter’s stop band starts: at 0.3π, not at $\pi/4$. Its wide transition, from 0.2π to 0.3π, does no harm, because whatever lies there lands at 0.8π or above, outside the kept band.

Compute only what you keep

The filter makes one output for every input sample, and then three of every four are thrown away, so why compute them? Each output is a sum of 75 products, and I only need every 4th one. Computing only those costs 75 multiplies per output, not $4\times75=300$. Polyphase structures (22.4) builds this into the filter’s structure.

In SciPy, upfirdn(h, x, 1, 4) filters with h and keeps every 4th output, and it computes only those. decimate(x, 4, ftype='fir') does the same job with its own filter, an 81-tap Hamming design with cutoff 0.25π, and by default it removes the filter’s delay. Without ftype it uses an order-8 Chebyshev IIR filter (Analog prototype filters, 20.1), which this page leaves out.

One stage or two

A bigger drop costs more. Take 96 kHz down to 12 kHz, so $M=8$, keeping 0 to 5 kHz with a 60 dB stop band. The stop-edge rule, in hertz, puts each stage’s stop edge at its output rate minus 5 kHz.

In one stage the stop edge is $12-5=7$ kHz. The transition is only 2 kHz wide at a 96 kHz rate, and Kaiser’s formula asks for 176 taps: 176 multiplies per output.

Now do it in two stages, 4 then 2. The first goes from 96 to 24 kHz, with its stop edge at $24-5=19$ kHz. Its transition, 5 to 19 kHz, is 14 kHz wide, so it needs only 26 taps. Whatever lies between 5 and 12 kHz after it is the second stage’s job.

The second stage goes from 24 to 12 kHz, with its stop edge at 7 kHz. Its transition is 2 kHz wide again, but at a 24 kHz rate that is four times wider in $\Omega$, so it needs 45 taps. The first stage runs twice for each final output, so the cost is $2\times26+45=97$ multiplies per output, 55 % of one stage. This is multistage decimation.

The table adds a third design: three stages that each halve the rate.

The early stages are short because their transition bands are wide: Kaiser’s length grows as $1/\Delta\Omega$ (19.1). Three stages cost 117, more than two here: the two halving stages cost $4\times11+2\times14=72$ per output, against $2\times26=52$ for the first stage of the two-stage design.

The table counts every tap. A halving stage can use a half-band filter (Special FIR filters, 19.4), in which every other tap is zero, and skip those multiplies.

Worked example

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

1. Where the tones land. Multiply 0.10π and 0.35π by $M$, then fold into 0 to π. To fold, subtract whole turns of $2\pi$ until the frequency is within $\pi$ of 0, and drop the sign.

At $M=8$ both tones land on 0.80π.

2. The decimation filter. For 60 dB, $\beta=0.1102\cdot51.3=5.6533$. The transition is $\Delta\Omega=0.1\pi=0.3142$ rad/sample, and the length formula gives $N_h-1\ge(60-7.95)/(2.285\cdot0.3142)=72.51$. So $N_h-1=73$ and $N_h=74$, which I make 75, with a delay of $(75-1)/2=37$ samples.

SciPy measures the 75 taps: within ±0.00111 of 1 up to 0.2π, at most −60.38 dB from 0.3π, and −68.14 dB at 0.35π. The alias after filtering is $-6.02-68.14=-74.16$ dB re 1. As in 19.1, every 4th tap from the centre is 0: 18 of the 75.

3. Stages. Kaiser’s formula gives 176 taps for one stage, and 26 and 45 for two. So the cost is 176 against $2\times26+45=97$ multiplies per output. At 12 000 outputs a second that is 2.11 million against 1.16 million multiplies a second. On the 100 MHz processor of Real-time processing (21.4), that is 2.1 % against 1.2 %.

Kaiser’s formula is an estimate, as 19.1 found. Checked one tap at a time, the stop band first reaches −60 dB at 181 taps for one stage, and at 28 and 45 taps for two. That is 181 against 101, and two stages still cost 56 %.

Where you’ll meet this

A sigma-delta converter (Oversampling and noise shaping, 11.3) runs its 1-bit loop at megahertz rates. Its digital low-pass is a decimator, built in stages, that brings the rate down to 48 kHz or whatever the output needs.

Software-defined radios sample a wide band fast and decimate down to the rate one channel needs. Audio tools decimate when they turn a 96 kHz recording into 48 kHz. Data loggers sample a sensor fast and store fewer, filtered samples.

The reverse, raising the rate, is Upsampling and interpolation (22.2). Changing the rate by any ratio is Resampling by any factor (22.3). Computing only the kept outputs, built into the structure, is Polyphase structures (22.4), and banks of decimating filters are Filter banks (23.1).

Reference card