Polyphase Decimation and the Cost of Changing Sample Rate

Reducing a sample rate by an integer factor $M$ requires lowpass filtering before discarding samples, because otherwise the content above the new Nyquist frequency folds into the retained band. The straightforward implementation filters at the original rate and then keeps every $M$-th output, which computes $(M-1)$ filtered samples for every one that survives. That waste is not marginal: for a decimation by ten the direct form spends ninety percent of its arithmetic on values that are immediately discarded. The polyphase decomposition eliminates exactly that waste, and the reason it can is a structural commutativity between filtering and downsampling that holds for any filter.

The discard is what makes the direct form expensive

A length-$N$ FIR filter applied at the input rate costs $N$ multiply-accumulates per input sample. To decimate by $M$, the filter runs at the input rate and its output is subsampled, so the cost per retained output sample is $MN$ multiply-accumulates. The filter’s stopband must extend from the new Nyquist frequency $f_s/(2M)$ to the old Nyquist frequency $f_s/2$, a transition band that is usually a small fraction of the old Nyquist frequency. A narrow transition band in normalized terms means a long filter: the length scales inversely with the transition width, so a decimation by ten with a sharp transition can require several hundred taps, and $MN$ becomes thousands of operations per output.

The obvious optimization is to compute only the outputs that are kept. The filter’s output at index $n$ is the inner product of the input with the reversed coefficient vector shifted by $n$. If only every $M$-th output is needed, only those inner products are needed, and each still costs $N$ multiplies. Computing one output in $N$ operations rather than $M$ outputs in $MN$ gives a factor of $M$ saving before any restructuring, just from not evaluating the discarded values. That factor is real but it does not reduce the $N$ per output, which is the part that matters for a long filter.

The polyphase decomposition splits the coefficients

Write the filter’s transfer function as a sum over its coefficients grouped by residue modulo $M$:

$$H(z) = \sum_{r=0}^{M-1} z^{-r} E_r(z^M)$$

where $E_r(z^M)$ contains the coefficients at indices $r, r+M, r+2M, \dots$ Each $E_r$ is a shorter filter of length about $N/M$, operating on the input at rate $f_s$ but with its own delay phase. In the time domain this is a bank of $M$ sub-filters, each of which produces one of the $M$ interleaved phases of the full filter output.

The key step is the identity that downsampling and filtering commute when the filter is expressed in this form. Because each sub-filter $E_r$ is followed by a downsampler, the downsampler can be moved ahead of the sub-filter: instead of computing $E_r$ at every input sample and keeping every $M$-th result, keep every $M$-th input sample and run $E_r$ on that reduced stream. The number of operations per retained output becomes $M \cdot (N/M) = N$, distributed across $M$ filters that each run at $1/M$ the input rate.

The arithmetic per output sample is now $N$ multiplies plus $N-1$ additions, the same count as a single-rate filter of the same length would need for one output, and no samples are computed and discarded. The saving relative to naive direct-form decimation is exactly the factor $M$. The memory access pattern also improves, because each sub-filter reads a contiguous decimated stream rather than striding through the input.

import numpy as np

def polyphase_decimate(x, h, M):
    """Decimate by M using h, computing only retained outputs."""
    x = np.asarray(x, dtype=float)
    h = np.asarray(h, dtype=float)
    n_out = (x.size - h.size) // M + 1
    y = np.zeros(n_out)
    for r in range(M):
        # Taps h[r + j*M] meet input phase M-r, one index earlier than the output.
        u = x[r::M] if r == 0 else np.concatenate(([0.0], x[M - r::M]))
        c = np.convolve(u, h[r::M])
        y[:min(n_out, c.size)] += c[:n_out]
    return y

h = np.hamming(81) / np.hamming(81).sum()
x = np.random.default_rng(3).standard_normal(10_000)
y = polyphase_decimate(x, h, 5)
print(y.size, abs(y).max())   # 1984 outputs, matching a full-rate filter

Each phase convolution runs on a stream $M$ times shorter than the input and uses a filter $M$ times shorter, so the $M$ phases together cost the same as one full-rate filter per retained sample rather than $M$ of them. The alignment is what makes the phases add coherently. Filter tap $h[r + jM]$ multiplies input sample $x[kM - r - jM]$ for output index $k$, and for $r > 0$ that sample sits one index earlier in input phase $M - r$; the single zero prepended to that phase stands for the sample that lies before the start of the record. Standard implementations avoid the $M$ separate convolutions and instead use a single ring buffer of length $N$ with a commutator that selects which coefficients apply to each output.

An equivalent formulation, the commutator model, describes the decimator as a rotating switch that distributes input samples to $M$ sub-filters in turn and sums their outputs at the reduced rate. The two descriptions are algebraically identical; the commutator model is the one that maps naturally onto hardware, where the switch is a counter and the sub-filters are parallel multiply-accumulate units.

The filter must be designed in the right band

The decimation filter’s specification is set by two frequencies. Its passband must cover the band that is to be preserved, up to some fraction of the new Nyquist frequency. Its stopband must begin at the new Nyquist frequency, because content above that will fold. The gap between the passband edge and the new Nyquist is the transition band, and it is what fixes the length.

The subtlety is that the specification is stated in the input-rate frequency axis, where the new Nyquist is at $f_s/(2M)$. If the desired passband extends to $0.4 f_s/M$ and the stopband must start at $0.5 f_s/M$, the transition occupies $0.1 f_s/M$ out of a total axis of $f_s/2$, a normalized width of $0.2/M$ cycles per sample. The stopband level itself comes from the window choice, and the trade between attenuation and transition width follows FIR Design by Windowing and the Gibbs Ripple. For large $M$ that normalized width is tiny, and the required length grows proportionally to $M$. This is the reason a single-stage decimation by a large factor is impractical: the filter length needed is proportional to the factor, so the per-output cost $N$ grows with $M$ even after polyphase decomposition.

The standard remedy is a cascade of stages, each decimating by a modest factor with a filter whose transition band is a fixed fraction of that stage’s input rate. If the total factor is $M = M_1 M_2 \cdots$, the cost per output becomes a sum rather than a product, and the sum is minimized by putting most of the decimation in the early stages where the filters can be short. The optimal stage distribution for a given transition bandwidth is a classical result, and for a multistage design the total arithmetic is typically far below the single-stage cost for the same total factor.

Upsampling is the mirror image, and rational rates combine both

Interpolation by $L$ inserts $L-1$ zeros between input samples and then filters at the output rate. The zero insertion makes the naive filter wasteful in the opposite way: most of its input samples are zero, so most of its multiplies are with zeros. The polyphase decomposition for interpolation computes the filter as $L$ sub-filters, each of which processes only the nonzero input samples and contributes to one output phase. The cost per output sample is $N/L$, which is again the factor $L$ improvement over the naive form.

A rational rate change by $L/M$ combines both operations: interpolate by $L$, filter, and decimate by $M$, with the filter running at the least common multiple rate. The polyphase structure handles the combination directly, with $L$ phases for the interpolation and $M$ decimation phases, and the filter’s coefficients are indexed by the product. The design constraint is that the filter’s stopband must be below the smaller of the two Nyquist frequencies involved, which for a large $L$ or $M$ again forces a multistage cascade.

rational resampler, L/M:
  x[n] -> [up by L] -> [h, polyphase] -> [down by M] -> y[k]
  effective rate = L*fs/M

For rate changes that are not rational, or that must track a continuously varying clock, the resampler becomes an interpolating filter with a fractional delay. A Farrow structure implements that: the filter coefficients are polynomials in the fractional delay, so a new delay value requires evaluating the polynomials rather than reloading a coefficient set. A cubic Farrow structure matches the first few terms of a Taylor expansion of the ideal interpolator and is accurate to a small fraction of a sample over a moderate bandwidth, at a cost of a few multiplies per output that does not grow with the interpolation factor.

What the structure does to the signal

Decimation lowers the sample rate, so it lowers the Nyquist frequency and therefore the bandwidth available at the output. Content between the new and old Nyquist frequencies is removed by the filter, and content above the new Nyquist that leaks through the stopband folds into the retained band. The stopband attenuation therefore sets an alias floor: a filter with 60 dB stopband attenuation lets through half a millionth of the out-of-band power, which for a strong interferer may still exceed the in-band signal.

There is also a phase consideration. The decimation filter is causal and has a group delay of about half its length at the input rate, which becomes half its length divided by $M$ at the output rate. The delay is constant if the filter is linear phase, so it does not distort, but it must be accounted for when aligning decimated and undecimated streams, and it changes the effective latency of the whole chain. In a feedback loop the delay appears at the reduced rate multiplied by the decimation factor, which can be a significant fraction of the loop’s period.

Finally, the order of operations matters for numerical range. Filtering before decimation operates on the full-rate signal, whose amplitude sets the accumulator width; decimating first would lose the very content the filter is meant to remove. The correct sequence is always filter, then decimate, and the polyphase structure is a way of evaluating that sequence without computing the discarded intermediate values. Reversing the order is a common and silent error, because the result still looks like a signal; it simply contains aliased content that no later stage can distinguish from genuine in-band energy. The sampling argument that makes the anti-alias filter necessary in the first place is developed in Sampling Rate, Aliasing, and the Folding of the Spectrum.