Multirate Signal Processing in Python: Decimation, Interpolation, and Polyphase Filters

Introduction

The Nyquist-Shannon sampling theorem guarantees that a band-limited signal can be perfectly reconstructed if it is sampled above the Nyquist rate. In practical systems, however, we frequently need to change the sampling rate in the middle of a processing chain:

• Audio format conversion: CD (44.1 kHz) ↔ professional audio and video (48 kHz)
• Delta-sigma ADCs: decimating a 1-bit stream oversampled at several MHz down to tens of kHz
• Communication systems: matching symbol rates to DAC/ADC rates, pulse shaping
• Wavelet transforms and filter banks: the structure itself downsamples by 2 at every stage

The discipline that handles systems with multiple coexisting sampling rates is multirate signal processing. This article systematically covers the frequency-domain behavior of downsampling (decimation) and upsampling (interpolation), rational-rate resampling, and the two ideas that make it all computationally efficient — polyphase decomposition and the Noble identities — with mathematical derivations and Python implementations using scipy.signal.decimate, resample, resample_poly, and upfirdn.

For the big picture of discrete-time signal processing, see the Discrete Signal Processing Fundamentals Roadmap; for the continuous-time side of interpolation, see Signal Reconstruction and Interpolation.

Downsampling (Decimation)

Definition and Frequency-Domain Analysis

Keeping only every \(M\) -th sample of a signal \(x[n]\) (downsampling by \(M\) ) is defined as:

The sampling rate drops from \(f_s\) to \(f_s / M\) . In the DTFT domain, the spectrum is stretched by a factor of \(M\) and \(M\) shifted copies are superimposed:

The \(k = 0\) term is the desired (stretched) spectrum; the \(k \neq 0\) terms are aliasing components. If the original signal is band-limited to \(|\omega| < \pi / M\) (i.e., \(f < f_s / 2M\) in physical frequency), the copies do not overlap and no aliasing occurs.

Why an Anti-Aliasing Filter Is Mandatory

Real signals usually contain energy above \(\pi / M\) , so a low-pass filter with cutoff \(\pi / M\) must precede the downsampler. The combination “LPF + downsampler” is called decimation:

where \(h[n]\) is the impulse response of the anti-aliasing filter. For the fundamentals of filter design, see Low-Pass Filter Design and Comparison and FIR vs IIR Filters.

Python: Naive Downsampling vs scipy.signal.decimate

scipy.signal.decimate bundles the anti-aliasing filter and the downsampler. Comparing it with naive slicing x[::M] makes the aliasing visible.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import decimate

# --- Test signal: desired 50 Hz component + 380 Hz high-frequency component ---
fs = 1000                      # original sampling rate [Hz]
t = np.arange(0, 2.0, 1/fs)
x = np.sin(2*np.pi*50*t) + 0.8*np.sin(2*np.pi*380*t)

M = 4                          # decimation factor: 1000 Hz -> 250 Hz (new Nyquist 125 Hz)
fs_new = fs / M

# (1) Naive slicing without a filter -> 380 Hz aliases down to 120 Hz
x_naive = x[::M]

# (2) decimate: anti-aliasing filter + downsampling
x_dec = decimate(x, M, ftype='fir', zero_phase=True)

def spectrum_db(sig, fs):
    """Hann-windowed magnitude spectrum in dB (normalized to the maximum)."""
    win = np.hanning(len(sig))
    X = np.abs(np.fft.rfft(sig * win))
    f = np.fft.rfftfreq(len(sig), 1/fs)
    return f, 20*np.log10(X / X.max() + 1e-12)

fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=False)

f0, S0 = spectrum_db(x, fs)
axes[0].plot(f0, S0)
axes[0].set_title('Original (fs=1000 Hz): peaks at 50 Hz and 380 Hz')
axes[0].set_xlim(0, 500)

f1, S1 = spectrum_db(x_naive, fs_new)
axes[1].plot(f1, S1, color='crimson')
axes[1].set_title('Naive x[::4] (fs=250 Hz): 380 Hz aliases to 120 Hz!')
axes[1].set_xlim(0, 125)

f2, S2 = spectrum_db(x_dec, fs_new)
axes[2].plot(f2, S2, color='seagreen')
axes[2].set_title('scipy.signal.decimate (fs=250 Hz): alias suppressed')
axes[2].set_xlim(0, 125)

for ax in axes:
    ax.set_ylabel('Magnitude [dB]')
    ax.set_ylim(-100, 5)
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

With naive slicing, the 380 Hz component folds around the new sampling rate \(f_s' = 250\) Hz: \(380 - 250 = 130\) Hz, then reflects at the new Nyquist (125 Hz) to appear as a spurious peak at \(250 - 130 = 120\) Hz, corrupting the band of interest. With decimate, this component is removed by the filter before downsampling.

Note that decimate defaults to an 8th-order Chebyshev type I IIR filter. With zero_phase=True (the default), filtfilt-style zero-phase filtering avoids phase distortion. If you want linear phase and explicit control of the filter characteristics, use ftype='fir'. For large factors (\(M > 13\) or so), the standard practice is to apply decimate in multiple stages (e.g., \(M = 32 \to 8 \times 4\) ).

Upsampling (Interpolation)

Zero Insertion and Spectral Images

Upsampling by \(L\) starts by inserting \(L - 1\) zeros between samples:

The spectrum becomes

i.e., the original spectrum compressed by \(1/L\) , with \(L\) copies (spectral images) centered at \(\omega = 2\pi k / L\) .

The Interpolation Filter

To remove the images, zero insertion is followed by a low-pass filter with gain \(L\) and cutoff \(\pi / L\) :

The gain of \(L\) compensates for the energy diluted by the inserted zeros. The combination “zero insertion + LPF” is called interpolation. The ideal filter is exactly sinc interpolation, connecting directly to the theory covered in Signal Reconstruction and Interpolation.

Python: Zero Insertion, upfirdn, and resample_poly

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample, resample_poly, upfirdn, firwin

# --- Test signal: 500 Hz + 1500 Hz ---
fs = 8000
t = np.arange(0, 0.2, 1/fs)
x = np.sin(2*np.pi*500*t) + 0.5*np.sin(2*np.pi*1500*t)

L = 3                          # 8000 Hz -> 24000 Hz
fs_up = fs * L

# (1) Zero insertion only -> images remain
x_zero = np.zeros(len(x) * L)
x_zero[::L] = x

# (2) Zero insertion + interpolation filter (gain L, cutoff pi/L) via upfirdn
h = L * firwin(121, 1/L)       # firwin cutoff is normalized to Nyquist
x_upfirdn = upfirdn(h, x, up=L)

# (3) resample_poly: polyphase implementation with an auto-designed Kaiser FIR
x_poly = resample_poly(x, L, 1)

# (4) FFT-based resample
x_fft = resample(x, len(x) * L)

def spectrum_db(sig, fs):
    win = np.hanning(len(sig))
    X = np.abs(np.fft.rfft(sig * win))
    f = np.fft.rfftfreq(len(sig), 1/fs)
    return f, 20*np.log10(X / X.max() + 1e-12)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))

f0, S0 = spectrum_db(x_zero, fs_up)
axes[0].plot(f0, S0, color='crimson')
axes[0].set_title('Zero insertion only: images at 6.5/8.5/9.5 kHz etc.')

f1, S1 = spectrum_db(x_upfirdn, fs_up)
axes[1].plot(f1, S1, color='seagreen')
axes[1].set_title('upfirdn (zero insertion + LPF): images removed')

f2, S2 = spectrum_db(x_poly, fs_up)
axes[2].plot(f2, S2, color='navy')
axes[2].set_title('resample_poly(x, 3, 1): same result, auto-designed filter')

for ax in axes:
    ax.set_ylabel('Magnitude [dB]')
    ax.set_ylim(-100, 5)
    ax.set_xlim(0, fs_up / 2)
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

With zero insertion alone, images appear at \(8000 \pm 500\) and \(8000 \pm 1500\) Hz (and beyond) at the 24 kHz rate; after the interpolation filter, only the original 500 / 1500 Hz components remain. For interpreting these spectra, see FFT-based spectral analysis and Window Functions and PSD.

Rational-Rate Resampling (L/M Conversion)

To convert the sampling rate by a rational factor \(L / M\) , connect upsampling (\(\uparrow L\) ) → low-pass filter → downsampling (\(\downarrow M\) ) in that order. Upsampling must come first — downsampling first would destroy information.

The single filter in the middle serves both as the interpolation filter (cutoff \(\pi / L\) ) and the anti-aliasing filter (cutoff \(\pi / M\) ):

The classic example is CD-to-48 kHz conversion, where \(48000 / 44100 = 160 / 147\) .

import numpy as np
import matplotlib.pyplot as plt
from math import gcd
from scipy.signal import resample_poly

# --- 44.1 kHz -> 48 kHz conversion ---
fs_in, fs_out = 44100, 48000
g = gcd(fs_out, fs_in)
L, M = fs_out // g, fs_in // g
print(f"L/M = {L}/{M}")        # -> 160/147

# Convert a 1 kHz sine and check the accuracy
duration = 0.05
t_in = np.arange(0, duration, 1/fs_in)
x_in = np.sin(2*np.pi*1000*t_in)

x_out = resample_poly(x_in, L, M)
t_out = np.arange(len(x_out)) / fs_out

# Compare against the reference (the continuous sine sampled directly at 48 kHz)
x_ref = np.sin(2*np.pi*1000*t_out)
edge = 200                     # exclude filter transients at the edges
err = np.max(np.abs(x_out[edge:-edge] - x_ref[edge:-edge]))
print(f"Max error: {err:.2e}")  # -> around 1e-5 (set by the filter accuracy)

plt.figure(figsize=(10, 4))
plt.plot(t_in[:100]*1000, x_in[:100], 'o-', label='44.1 kHz input',
         markersize=4, alpha=0.7)
plt.plot(t_out[:109]*1000, x_out[:109], '.-', label='48 kHz output',
         markersize=4, alpha=0.7)
plt.xlabel('Time [ms]')
plt.ylabel('Amplitude')
plt.title('Rational resampling 44.1 kHz -> 48 kHz (L/M = 160/147)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

A single call resample_poly(x, 160, 147) performs high-quality rate conversion. Internally, the “upsample by 160, then downsample by 147” chain is implemented in a polyphase structure, costing roughly \(1/160\) of a naive implementation. The next section explains how.

Polyphase Decomposition

The Waste in the Naive Implementation

Implementing the decimation of Eq. (3) naively means “run the full FIR filter at the input rate \(f_s\) , then keep only one output out of \(M\) ” — you compute and throw away \(M-1\) of every \(M\) outputs. For an FIR filter with \(N\) taps, that costs \(N f_s\) multiply-accumulate (MAC) operations per second.

Type-1 Polyphase Decomposition

Polyphase decomposition eliminates this waste. Split the filter \(H(z) = \sum_n h[n] z^{-n}\) into \(M\) subfilters (polyphase components) by taking every \(M\) -th coefficient:

For example, with \(M = 3\) and \(h = [h_0, h_1, h_2, h_3, h_4, h_5]\) , the components are \(E_0 = [h_0, h_3]\) , \(E_1 = [h_1, h_4]\) , \(E_2 = [h_2, h_5]\) .

The Polyphase Decimator

Substituting Eq. (8) into Eq. (3) and rearranging, the decimator output becomes

In words: distribute the input into \(M\) low-rate streams \(x_p\) (a commutator), pass each stream through its short subfilter \(E_p\) , and sum. All filtering now runs at the output rate \(f_s / M\) , so the MAC count per second drops to

a reduction by a factor of \(M\) . Interpolators can be polyphased in exactly the same way, saving a factor of \(L\) (multiplications by inserted zeros are skipped). Rational \(L/M\) converters benefit from both savings at once.

Python: Verifying a Polyphase Decimator

We implement Eq. (10) directly and confirm it matches both the naive “filter, then downsample” approach and scipy.signal.upfirdn.

import numpy as np
from scipy.signal import firwin, upfirdn

def polyphase_decimate(x, h, M):
    """Decimation in polyphase form (direct implementation of Eq. (10))."""
    N_out = len(x) // M
    y = np.zeros(N_out)
    for p in range(M):
        e_p = h[p::M]                        # subfilter E_p
        # low-rate input stream x_p[m] = x[mM - p]
        idx = np.arange(N_out) * M - p
        x_p = np.zeros(N_out)
        valid = idx >= 0
        x_p[valid] = x[idx[valid]]
        y += np.convolve(x_p, e_p)[:N_out]   # convolution at the low rate
    return y

# --- Verification ---
rng = np.random.default_rng(42)
fs, M = 48000, 4
x = rng.standard_normal(fs)                  # 1 second of white noise
h = firwin(128, 1/M)                         # FIR with cutoff pi/M

y_naive = np.convolve(x, h)[::M][:len(x)//M] # naive: full-rate filter -> downsample
y_poly = polyphase_decimate(x, h, M)
y_scipy = upfirdn(h, x, up=1, down=M)[:len(x)//M]

print("polyphase == naive :", np.allclose(y_poly, y_naive))   # True
print("polyphase == scipy :", np.allclose(y_poly, y_scipy))   # True

# --- Complexity comparison (theoretical) ---
N = len(h)
print(f"Direct form     : {N * fs / 1e6:.1f} M MAC/s")
print(f"Polyphase form  : {N * fs / M / 1e6:.1f} M MAC/s (1/{M})")

All three implementations produce identical outputs. scipy.signal.upfirdn and resample_poly implement this polyphase structure in C and never materialize the zero-stuffed intermediate signal.

Noble Identities

The correctness of the polyphase structure rests on the Noble identities — two equivalences for exchanging the order of filters and rate changers.

Downsampling side: downsampling by \(M\) followed by \(H(z)\) is equivalent to \(H(z^M)\) followed by downsampling by \(M\) :

Upsampling side: \(H(z)\) followed by zero insertion by \(L\) is equivalent to zero insertion followed by \(H(z^L)\) :

Here \(H(z^M)\) denotes the filter whose impulse response has \(M - 1\) zeros inserted between consecutive coefficients. Applying Eq. (12) to each term \(z^{-p} E_p(z^M)\) of the polyphase decomposition (Eq. 8) moves \(E_p(z^M)\) behind the downsampler — to the low-rate side — yielding the efficient structure of Eq. (10). The golden rule of multirate design, “always run filters on the low-rate side,” is guaranteed by these identities. The same rewriting is used to streamline multistage structures in wavelet filter banks and delta-sigma decimators.

Choosing Between decimate, resample, and resample_poly

Here is a summary of SciPy’s three resampling functions.

Practical guidelines:

• When in doubt, use resample_poly: few artifacts even on aperiodic real data, and it handles rational ratios
• Use resample only when the signal is periodic or well band-limited; watch out for edge ringing
• decimate is a convenience function for integer-factor downsampling; apply it in stages for large \(M\)
• If you need to control stopband attenuation or transition width yourself, design the FIR with firwin etc. and pass it to upfirdn (see the Digital Filter Design Guide for the design flow)

Applications and Design Tips

Combining a bandpass filter with downsampling — bandpass decimation — is a standard technique for drastically reducing the computational load when analyzing narrowband signals.

Summary

• Downsampling stretches the spectrum by \(M\) and superimposes copies, so an anti-aliasing filter with cutoff \(\pi/M\) is mandatory beforehand (= decimation)
• Upsampling creates spectral images through zero insertion, which are removed by an interpolation filter with gain \(L\) and cutoff \(\pi/L\)
• Rational \(L/M\) conversion is built as “\(\uparrow L\) → LPF → \(\downarrow M\) ”, with a single filter of cutoff \(\min(\pi/L, \pi/M)\) serving both roles
• Polyphase decomposition splits the filter into \(M\) subfilters running on the low-rate side, cutting the computation to \(1/M\) (or \(1/L\) for interpolators)
• The Noble identities justify exchanging filters and rate changers, underpinning the design principle “filters belong on the low-rate side”
• In SciPy, resample_poly is the practical default; use decimate for integer-factor downsampling, resample for periodic signals, and upfirdn for custom filters

Related Articles

• Nyquist-Shannon Sampling Theorem and Aliasing - The Nyquist frequency and aliasing theory that multirate processing builds on.
• Signal Reconstruction and Interpolation: Sinc, Linear, and Spline - A detailed treatment of sinc interpolation, the ideal interpolation filter.
• DTFT vs DFT vs FFT - The frequency-domain representations underlying the spectral analysis in this article (Eqs. 2 and 5).
• Discrete Signal Processing Fundamentals Roadmap - An integrated hub covering sampling, interpolation, DFT, and the Z-transform.
• FIR vs IIR Filters - Directly relevant to choosing the implementation of anti-aliasing/interpolation filters.
• Low-Pass Filter Design and Comparison - The basics of designing the LPF that precedes decimation.
• Digital Filter Design Guide - The design flow for custom filters to use with upfirdn.
• Fast Fourier Transform (FFT): Theory and Python Implementation - Used to verify the spectra of resampled signals.
• Window Functions and Power Spectral Density (PSD) - The theory behind the window functions used in the spectral plots here.
• Bandpass Filter Design: Theory and Python Implementation - The band-selection filter to combine with bandpass decimation.
• Time-Frequency Analysis Guide - Choosing analysis methods (STFT, wavelets, etc.) for the downsampled signal.

References

• Crochiere, R. E., & Rabiner, L. R. (1983). Multirate Digital Signal Processing. Prentice Hall.
• Vaidyanathan, P. P. (1993). Multirate Systems and Filter Banks. Prentice Hall.
• Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
• harris, f. j. (2004). Multirate Signal Processing for Communication Systems. Prentice Hall.
• SciPy Signal Processing — decimate
• SciPy Signal Processing — resample_poly
• SciPy Signal Processing — upfirdn