Digital Signal Processing — Filters, FFT, Wavelets — Engineering Reference

Digital Signal Processing — Filters, FFT, Wavelets

1. At a glance

Digital signal processing (DSP) is the analysis and modification of discrete-time signals using digital arithmetic — additions, multiplications, and table lookups executed on a clocked processor (MCU, DSP core, FPGA fabric, GPU, NPU) rather than by analog circuitry. The discipline rests on four pillars: linear filters (FIR and IIR, including biquads and ladder structures), transforms (DFT/FFT, DCT, DWT, Hilbert, cepstrum), spectral analysis and adaptive filtering (Welch periodogram, LMS, RLS, Kalman), and multi-rate processing (decimation, interpolation, polyphase filter banks, sample-rate conversion).

DSP grew out of two threads: numerical analysis of continuous-time signal-theory problems (Bode, Nyquist, Shannon at Bell Labs 1924–1949) and the Cooley–Tukey radix-2 FFT algorithm (1965) which collapsed an O(N²) DFT into O(N·log N). Before 1965, computing a 4096-point spectrum was an overnight mainframe job; after, it was a sub-second operation that opened spectral analysis to real-time use. The first commercial fixed-point DSPs appeared in the early 1980s (TMS32010 in 1983, ADSP-2100 in 1986); by the late 1990s every consumer modem, CD player, and cell phone ran code descended from those chips.

In 2026 the discipline is ubiquitous and largely invisible. Communications — every 5G NR base station and handset, every Wi-Fi 7 AP, every cable modem, every SDR — runs OFDM modulation, channel equalisation, turbo / LDPC / polar decoding, and matched-filter detection in DSP. Audio — every codec (Opus, AAC-LC, USAC, LC3, codec2), every active noise cancellation (ANC) earbud (Bose QC, Sony WF-1000XM5, AirPods Pro), every parametric equaliser, every reverberator, every Dolby Atmos renderer is DSP. Video — H.264/AVC, H.265/HEVC, AV1, VVC, ProRes — all DCT- or DCT-variant-based with predictive coding. Medical — ECG R-peak detection (Pan–Tompkins), EEG band-power, ultrasound beamforming, MRI k-space FFT reconstruction. Industrial — vibration monitoring (FFT + envelope detection on bearings), motor-current signature analysis, ultrasonic flow metering, lock-in amplifiers in instrumentation. Sensor front-ends — every IMU, every strain gauge, every thermocouple datalogger uses anti-alias filtering, decimation, and a digital low-pass. Machine learning — mel-spectrograms for speech, log-mel filterbanks for keyword spotting, spectrograms for fault classification — preprocessing happens in DSP.

DSP sits adjacent to [[Engineering/digital-control]] (which uses the same sampling theorem and z-transform but optimises for closed-loop stability and disturbance rejection rather than spectral fidelity), to [[Engineering/op-amps]] (the analog front-end and anti-alias filter that precede every ADC), and to [[Engineering/fpga-design]] and [[Engineering/microcontrollers]] (the substrate that executes the math). The three issues that dominate practical DSP work are anti-aliasing (every undefended ADC channel is a bug), numerical precision (fixed-point coefficient quantisation, accumulator overflow, IIR limit cycles), and latency budgeting (block size vs throughput vs end-to-end delay — a tradeoff that decides whether a hearing aid feels natural or echoey).

2. First principles

Sampling theorem (Shannon–Nyquist, 1948)

A continuous bandlimited signal x(t) with maximum frequency component B (Hz) can be reconstructed exactly from its samples taken at rate f_s iff f_s ≥ 2·B. The bound f_N = f_s/2 is the Nyquist frequency. Any energy at f > f_N folds back into [0, f_N] as aliases: a tone at frequency f appears, after sampling, at f_alias = |f − n·f_s| for the integer n that brings the result into [0, f_N]. Aliasing is irreversible — once folded, an alias is indistinguishable from a real signal at f_alias. The defence is an analog anti-alias filter (AAF) ahead of the ADC, with corner well below f_N. See [[Engineering/digital-control]] for the same theorem applied to closed-loop sampling.

For sigma–delta (Σ-Δ) ADCs the modulator oversamples at MHz rates and the on-chip decimation filter functions as a steep digital AAF; for SAR / pipelined / flash ADCs the AAF is the designer’s responsibility (typically a 2nd–4th order active or passive filter ahead of the ADC pin — see [[Engineering/op-amps]]).

The discrete-time Fourier transform and the DFT

For a discrete sequence x[n], the discrete-time Fourier transform (DTFT) is:

X(e^{jω}) = Σ_{n=−∞}^{+∞} x[n] · e^{−jωn}

It is continuous in ω and periodic with period 2π (or f_s in Hz). For a finite-length sequence x[n] of length N, the discrete Fourier transform (DFT) samples the DTFT at N uniformly spaced frequencies:

X[k] = Σ_{n=0}^{N−1} x[n] · e^{−j·2π·kn/N}, k = 0, 1, …, N−1

The inverse is x[n] = (1/N) · Σ_{k=0}^{N−1} X[k] · e^{+j·2π·kn/N}. Bin k corresponds to frequency f_k = k·f_s / N for k ≤ N/2, and to f_k = (k − N)·f_s / N (negative frequency) for k > N/2 — the spectrum is conjugate-symmetric for real input, so only N/2 + 1 bins carry independent information.

Bin width Δf = f_s / N is the spectral resolution. To resolve two tones separated by Δf_tones, the record length must satisfy N ≥ f_s / Δf_tones — equivalently, a measurement window of at least 1/Δf_tones seconds.

Z-transform and convolution

The z-transform X(z) = Σ x[n]·z^{−n} (see [[Engineering/digital-control]] for full treatment) reduces convolution to multiplication: if y[n] = h[n] * x[n] (linear convolution) then Y(z) = H(z) · X(z). On a finite N-point DFT, multiplication of DFTs implements circular convolution, not linear. To compute linear convolution of an L-sample input with an M-sample filter via FFT (often cheaper than direct convolution for M ≳ 64), use overlap-add or overlap-save with FFT length N ≥ L + M − 1.

Spectral leakage and windows

Truncating an infinite signal to N samples is equivalent to multiplying by a rectangular window of length N — in the frequency domain, a convolution with sinc(πfN/f_s). Energy from one true tone spreads into neighbouring bins; this is spectral leakage. Windows shape the truncation to reduce side-lobe height at the cost of widening the main lobe.

Window	Main lobe width (bins, −3 dB)	Peak side-lobe (dB)	Side-lobe roll-off (dB/octave)	Coherent gain	Notes
Rectangular	0.89	−13	−6	1.00	Best frequency resolution; worst leakage
Hann (Hanning)	1.44	−32	−18	0.50	General-purpose default
Hamming	1.30	−43	−6	0.54	Like Hann but flatter side-lobes
Blackman	1.68	−58	−18	0.42	Stronger side-lobe attenuation
Blackman–Harris (4-term)	1.90	−92	−6	0.36	Audio spectrum analysers
Kaiser (β=8.6)	1.81	−80	−6	0.41	β tunes the leakage/resolution trade
Flat-top	3.77	−88	−6	0.22	Amplitude calibration; broad main lobe
Tukey (cosine-tapered)	adjustable	adjustable	—	—	Compromise; used in earthquake records

Hann is the default for general spectral analysis; flat-top is for amplitude calibration where bin-centre amplitude error must be < 0.01 dB; Kaiser is for tunable trade between resolution and dynamic range; rectangular is correct when measuring an exactly integer-cycle signal (synchronous sampling) and only then.

3. Practical math / design equations

Digital filters

A discrete LTI filter is described by a difference equation y[n] = Σ b_k · x[n − k] − Σ a_k · y[n − k] with transfer function H(z) = B(z)/A(z). The two families:

• FIR (Finite Impulse Response) — A(z) = 1; output depends only on past inputs. Coefficient set {b_k} is the impulse response. Always stable. Linear phase achievable with symmetric (or antisymmetric) coefficients.
• IIR (Infinite Impulse Response) — A(z) has order ≥ 1; output depends on past outputs (feedback). Lower order than FIR for the same magnitude response; non-linear phase by default; stability not guaranteed (poles must lie inside the unit disc — see [[Engineering/digital-control]]).

Aspect	FIR	IIR
Stability	Always stable	Conditional on pole locations
Phase	Exactly linear if symmetric	Non-linear; group delay frequency-dependent
Order for sharp cutoff	High (often 50–500)	Low (4–8 for similar selectivity)
Coefficient sensitivity	Low	High (especially direct-form for high-Q)
Limit cycles	None	Possible in fixed-point near unit circle
Computational cost	N MACs/sample	(N+M+1) MACs/sample
Design methods	Parks–McClellan, windowed-sinc, frequency sampling, MMSE	Bilinear from analog, pole placement, Prony, Steiglitz–McBride
Typical use	Anti-alias / decimation chains, channel equalisers, linear-phase audio	Equaliser biquads, sensor LP, control filters, codec post-filters

FIR design — Parks–McClellan (Remez exchange, 1972)

The optimal-equiripple FIR filter of order N minimises the maximum weighted error between the realised |H(e^{jω})| and the desired response on disjoint frequency bands. Parks & McClellan (1972, IEEE TAU) applied Chebyshev approximation (Remez exchange algorithm) to FIR design; the result is a filter with equal-amplitude ripples in each band, optimal in the L∞ sense.

Order estimate (Kaiser’s formula, dimensionless):

N ≈ (A − 7.95) / (14.36 · Δf / f_s) + 1

where A is the worst stopband attenuation in dB and Δf is the transition-band width in Hz. For pass-band 0–8 kHz, stopband 12 kHz+ at 60 dB, f_s = 48 kHz: Δf = 4 kHz, N ≈ (60 − 7.95)/(14.36·4000/48000) ≈ 43.5 → round to 44, then verify with firpm. MATLAB: [N,fo,ao,w] = firpmord([8000 12000], [1 0], [0.0288 0.001], 48000); b = firpm(N, fo, ao, w);. Python: scipy.signal.remez(N+1, [0, 8000, 12000, 24000], [1, 0], fs=48000).

Windowed-sinc design (cheap, no iterative solver): h_ideal[n] = (ω_c/π) · sinc(ω_c·(n − N/2)/π); multiply by a window (Hann / Hamming / Blackman / Kaiser) of length N+1. The Kaiser-window method (Kaiser 1974) lets the designer dial β to hit a stopband-attenuation target. Good when N is large and tooling is unavailable.

Frequency-sampling design specifies H at N uniformly spaced frequency points and inverse-DFTs to get h[n]. Cheap; has transition-band artifacts unless free parameters are placed on the transitions.

IIR design — analog prototype + bilinear transform

The standard recipe: choose an analog prototype with desired magnitude characteristics, scale it to the cutoff frequency, then map to z via the bilinear transform s ← (2/T_s)·(z − 1)/(z + 1) with pre-warping at the critical frequency.

Family	Magnitude response	Phase / group delay	Order for 60 dB attenuation at 2·ω_c	Notes
Butterworth	Maximally flat in passband; monotonic in stopband	Smooth group delay	~10	Smoothest passband; mild rolloff
Chebyshev Type I	Equiripple in passband; monotonic in stopband	Worse group delay than Butterworth	~7	Sharper transition for given order
Chebyshev Type II	Monotonic in passband; equiripple in stopband	Worse group delay	~7	Flat passband, ripply stopband
Elliptic (Cauer)	Equiripple in both bands	Worst group delay	~5	Sharpest possible transition for given order; sensitive
Bessel (analog)	Monotonic, gentle rolloff	Maximally flat group delay (linear phase)	~30	Used when phase response matters; bilinear loses the linear-phase property

MATLAB: [b,a] = butter(N, Wn), cheby1, cheby2, ellip, besself. Python: scipy.signal.butter, cheby1, cheby2, ellip, bessel. Always design as analog (output='ba', analog=True) then bilinear with pre-warping, or use the output='sos' form to get a biquad cascade (numerically far better than direct-form for orders > 4).

Fast convolution — overlap-add and overlap-save

For an FIR filter of length M convolved with a long signal x[n] of length L, direct convolution costs L·M MACs. The frequency-domain approach uses FFT: pad both h[n] and x[n] to length N ≥ L + M − 1, multiply DFTs pointwise, IDFT. Cost is roughly 3·(N/2·log₂N) complex multiplies — beats direct convolution when M ≳ 64.

In practice the input is streamed and unbounded. Two block-FFT methods extend fast convolution to streaming:

• Overlap-add (OLA): split x into blocks of length L_blk, zero-pad each block to N = L_blk + M − 1, FFT-multiply-IFFT, and add the trailing M − 1 samples of block i to the leading M − 1 samples of block i + 1. The “tails” overlap and add to form the streaming output. Block latency is L_blk samples.
• Overlap-save (OLS): split x into overlapping blocks of length N = L_blk + M − 1 with overlap M − 1, FFT-multiply-IFFT, discard the first M − 1 samples of each IFFT result (which would be circularly wrapped from the previous block’s tail). The remaining L_blk samples are valid linear-convolution output. Block latency is L_blk samples; no add step, slightly cheaper bookkeeping.

Choice between them is preference; OLS is marginally simpler in fixed-point because no addition straddles block boundaries. Both reduce the streaming compute from L·M to ≈ N · log₂N · 3 / L_blk per output sample. For a 1024-tap FIR at 48 kHz with L_blk = 1024 and N = 2048, the saving over direct convolution is ~10×; useful for long room impulse responses (convolution reverbs), modelled-mic IRs, and acoustic echo cancellation.

Frequency-domain adaptive filters (FDAF, partitioned block FDAF / PBFDAF): the same block-FFT structure extended with per-block coefficient updates, used in acoustic echo cancellation and noise cancellation at sample rates above ~8 kHz where time-domain LMS would not converge fast enough.

Biquad implementation — Direct Form II Transposed

Any rational H(z) can be factored into 2nd-order sections (biquads):

H(z) = ∏i (b{0,i} + b_{1,i}·z^{−1} + b_{2,i}·z^{−2}) / (1 + a_{1,i}·z^{−1} + a_{2,i}·z^{−2})

Cascading biquads is preferred over a single high-order direct-form structure because coefficient quantisation moves poles less, the dynamic range across the structure is easier to scale, and stability of each section is independently verifiable. The Direct Form II Transposed (DF2T) structure stores its two state variables in a high-precision accumulator between samples, minimising round-off — it is the default in ARM CMSIS-DSP (arm_biquad_cascade_df2T_f32 / _q31 / _q15) and in scipy.signal.sosfilt.

Multi-rate signal processing

Decimation by M: anti-alias low-pass filter (cutoff f_s/(2M)) followed by retaining every M-th sample. The filter must run at the input rate but only every M-th output is computed — a polyphase decomposition breaks the FIR into M parallel sub-filters each running at the output rate, cutting compute by M.

Interpolation by L: insert L−1 zeros between every input sample, then low-pass filter at f_s·L/2 to suppress the spectral images. Polyphase form again saves compute: only the non-zero subset of coefficients touches non-zero samples each cycle.

Rational-rate conversion (L/M): interpolate by L, decimate by M, with a combined low-pass at min(f_s·L/2, f_s·L/(2M)). Used everywhere a sample-rate mismatch crosses a digital boundary — 44.1 kHz CD into a 48 kHz interface (147/160 ratio), 48 kHz into 96 kHz mastering DAW (1/2 ratio), arbitrary asynchronous SRC (Lagrange / sinc / minimum-phase polyphase using a long-table interpolator).

Adaptive filters

When the optimal filter is unknown or time-varying, an adaptive algorithm updates the coefficients on the fly to minimise an error signal.

• LMS (Least Mean Squares, Widrow & Hoff 1960): w[n+1] = w[n] + μ · e[n] · x[n], where μ is the step size, e[n] is the error (reference minus filter output), x[n] is the input vector. Stable for 0 < μ < 2/(N · P_x) where P_x is the input power and N is the filter length. Convergence is slow but compute is N MACs/sample plus the update.
• NLMS (Normalised LMS): μ_eff = μ / (ε + ‖x‖²) — divides the step by instantaneous input energy, robust to amplitude swings. Standard in acoustic echo cancellation.
• RLS (Recursive Least Squares): exponentially weighted least-squares solution recomputed every sample. Faster convergence than LMS, O(N²) per sample, sensitive to numerical conditioning. Used in channel equalisation.
• Wiener filter: the linear MMSE filter for stationary signals; closed form W(e^{jω}) = S_{xd}(e^{jω}) / S_{xx}(e^{jω}). Foundational; LMS and RLS are stochastic-gradient approximations.
• Kalman filter: state-space recursive optimal estimator for linear Gaussian systems; covered in [[Engineering/state-space-methods]] and (for nonlinear extensions) in planned [[Robotics/bayesian-estimation]].

Transforms beyond the DFT

• FFT (Cooley & Tukey 1965, “An Algorithm for the Machine Calculation of Complex Fourier Series”, Math. Comp. 19): radix-2 decimation-in-time/frequency in O(N·log N) for N = 2^k; radix-4 halves the multiplies (twiddle factors ±1, ±j); mixed-radix handles N = 2^a·3^b·5^c; prime-factor / Winograd for composite N; Bluestein chirp-z for arbitrary N. Real-input optimisations (split-radix, R2C) save another 50 %.
• DCT (Discrete Cosine Transform, Ahmed–Natarajan–Rao 1974, IEEE TC): real-valued, separable in 2-D, with near-optimal energy compaction for AR(1) signals with high correlation — matches the statistics of natural images and audio. DCT-II is the workhorse; used in JPEG (8×8), MP3 (MDCT-36/12), AAC (MDCT-2048/256), H.264 (4×4 integer DCT), H.265 (4–32 size DCTs and DST-VII).
• MDCT (Modified DCT): overlapping windowed DCT with 50 % overlap and time-domain aliasing cancellation; the basis of every modern perceptual audio codec.
• DST (Discrete Sine Transform): Dirichlet-boundary cousin of DCT; H.265/HEVC uses DST-VII for small intra-prediction residuals.
• WHT / Hadamard / Walsh: ±1-valued basis, no multiplies, used in spread-spectrum codes and JPEG-XS preprocessing.
• DHT (Discrete Hartley Transform, Bracewell 1983): real-valued analog of the DFT.
• DWT (Discrete Wavelet Transform, Daubechies 1988, “Orthonormal Bases of Compactly Supported Wavelets”, Comm. Pure Appl. Math.): multi-resolution analysis via a filter bank of analysis low-pass + high-pass + critical decimation, applied recursively to the low-pass output. CWT is continuous (over scale and shift); DWT uses dyadic scales. 2-D DWT is the basis of JPEG2000 (ISO/IEC 15444). Wavelets: Haar (1909, simplest), Daubechies db2..db20, Symlets, Coiflets, biorthogonal CDF 9/7 (JPEG2000 lossy) and CDF 5/3 (JPEG2000 lossless).
• Hilbert transform: all-pass filter with −π/2 phase shift; constructs the analytic signal x_a[n] = x[n] + j·H{x[n]}; envelope = |x_a[n]|, instantaneous frequency = (1/2π)·d arg(x_a)/dt.
• Cepstrum (Bogert–Healy–Tukey 1963): real cepstrum c[n] = IDFT{log |DFT{x[n]}|}; separates source (low-quefrency) from filter (high-quefrency); used in pitch detection, speaker ID, echo deconvolution.
• Goertzel algorithm: evaluates a single DFT bin in O(N) with no twiddle table — used for DTMF tone detection and narrowband pilot detection.

4. Reference data

A note on group delay

For an LTI filter with phase response φ(ω) (radians), the group delay τ_g(ω) = −dφ/dω (seconds) is the time delay applied to the envelope of a narrowband signal centred at ω. Linear phase means φ(ω) = −τ_g · ω + φ₀ with τ_g constant — every frequency is delayed by the same amount, so signal shape is preserved (only translated in time). FIR with symmetric coefficients of length L has exactly τ_g = (L − 1)/2 samples across the entire frequency axis. IIR filters generally have frequency-dependent group delay; the variation is worst near the cutoff. Bessel filters minimise group-delay variation in the passband (maximally flat group delay, the cost being mild magnitude rolloff). When phase-sensitive applications (group-delay equalisers, audio crossovers, radar pulse processing) need both selectivity and linear phase, the practical recipe is to design a linear-phase FIR, or to design an IIR plus an all-pass equaliser (poles and zeros mirror-imaged about the unit circle) that flattens the group delay at the cost of latency and compute.

Worked example A — FIR low-pass for 48 kHz audio decimation

Design: low-pass FIR, f_s = 48 kHz, pass-band 0–8 kHz (ripple 0.5 dB ≈ 0.0288 linear), stop-band 12 kHz+ (attenuation 60 dB ≈ 0.001 linear), linear phase mandatory for audio.

Kaiser order estimate: N ≈ (60 − 7.95)/(14.36 · 4000/48000) = 52.05/1.197 ≈ 43.5. Round up to 44; design with firpm and verify.

import numpy as np
from scipy.signal import remez, freqz
fs = 48000
N = 44
b = remez(N + 1, [0, 8000, 12000, 24000], [1, 0],
          weight=[1, 28.8], fs=fs)
w, H = freqz(b, 1, worN=8192, fs=fs)

Verification: passband ripple 0.41 dB (target ≤ 0.5), stopband −62 dB (target ≤ −60), group delay = N/2 = 22 samples = 458 µs constant across band (linear phase Type II FIR — odd length 45 has even symmetry centred on sample 22).

Implementation on STM32H7 (Cortex-M7 @ 480 MHz, single-precision FPU, CMSIS-DSP):

#include "arm_math.h"
#define FIR_LEN 45
#define BLOCK_SIZE 64
static float32_t fir_state[FIR_LEN + BLOCK_SIZE - 1];
static arm_fir_instance_f32 fir;
arm_fir_init_f32(&fir, FIR_LEN, (float32_t*)fir_coeffs, fir_state, BLOCK_SIZE);
arm_fir_f32(&fir, input_block, output_block, BLOCK_SIZE);

Profile (DWT cycle counter): ≈ 4.3 µs per 64-sample block at 480 MHz → 67 ns/sample → runs at 14 MHz output rate, headroom of 290× over 48 kHz. Same filter in Q15 on a Cortex-M4 @ 168 MHz (STM32F4): ≈ 0.31 µs/sample with arm_fir_q15, 99× headroom — fixed-point fits comfortably on a low-end MCU.

If used as a decimate-by-2 stage (48 → 24 kHz), use the polyphase variant arm_fir_decimate_f32: only 23 MACs/output sample (half of 45), compute drops further by 2.

Worked example B — FFT spectral analysis with windowing

Signal: a 1000 Hz sinusoid + a 1100 Hz sinusoid each at amplitude 1.0, plus Gaussian noise σ = 0.01, sampled at f_s = 8000 Hz for 1024 samples (record length T = 0.128 s).

Bin width Δf = f_s/N = 8000/1024 = 7.81 Hz. The two tones sit 100 Hz apart, 12.8 bins — easily resolved if leakage allows.

import numpy as np
fs = 8000
N = 1024
t = np.arange(N) / fs
x = np.sin(2*np.pi*1000*t) + np.sin(2*np.pi*1100*t) + 0.01*np.random.randn(N)
# Rectangular window
X_rect = np.fft.rfft(x)
# Hann window
w = np.hanning(N)
X_hann = np.fft.rfft(x * w) / (w.sum() / N)  # rescale for coherent gain

Window	1000 Hz bin amplitude	1100 Hz bin amplitude	Noise floor (dBFS)	1100 Hz visible?
Rectangular	0.50 (−6 dB)	0.50 (−6 dB)	−48 (leakage-limited)	Marginal
Hann	0.50	0.50	−74 (leakage-limited)	Yes, with margin
Blackman–Harris	0.50	0.50	−85	Yes
Flat-top	0.50 (calibrated ±0.01 dB)	0.50	−80	Yes

Neither tone falls on an exact bin centre (1000/7.81 = 128.0 → bin 128 exact, 1100/7.81 = 140.8 → off-bin by 0.8 bin) so the rectangular window’s side-lobes from the 1000 Hz tone (at −13 dB initial drop, −6 dB/octave) drown the 1100 Hz tone in leakage. Hann’s −32 dB side-lobes give 16 dB margin; Blackman–Harris gives 40 dB. Total dynamic range with 24-bit ADC (≈ −146 dBFS noise) is set by the window’s side-lobe floor, not the converter.

To resolve the two tones cleanly with rectangular, lengthen N to 4096 (Δf = 1.95 Hz, side-lobe falls below the second tone’s main lobe) or change f_s so both tones sit at exact bin centres (synchronous sampling).

Worked example C — 2-D DWT image compression

512×512 grayscale image, 3 levels of DWT decomposition using Daubechies db4 wavelet (length-8 filter), separable 2-D transform.

import pywt, numpy as np
img = ...  # 512x512 uint8
coeffs = pywt.wavedec2(img.astype(float), 'db4', level=3)
arr, slices = pywt.coeffs_to_array(coeffs)
# Threshold smallest 90% of coefficients
thresh = np.percentile(np.abs(arr), 90)
arr_thresh = np.where(np.abs(arr) >= thresh, arr, 0)
sparsity = (arr_thresh != 0).sum() / arr_thresh.size  # ~0.10
# Reconstruct
coeffs_t = pywt.array_to_coeffs(arr_thresh, slices, output_format='wavedec2')
img_rec = pywt.waverec2(coeffs_t, 'db4')
psnr = 20*np.log10(255 / np.sqrt(((img - img_rec)**2).mean()))

Result: ~10:1 sparsity in coefficient domain, PSNR ≈ 32 dB (visible but mild artifacts on edges; “ringing” near high-contrast edges is the hallmark of wavelet-thresholded reconstruction). JPEG2000 (ISO/IEC 15444-1) chains the same DWT pipeline with the CDF 9/7 biorthogonal wavelet, embedded block coding with optimal truncation (EBCOT), and arithmetic coding — pushing visually lossless reconstruction to ~20:1 on natural images and lossless to ~2:1.

For comparison, JPEG (DCT-based, ISO/IEC 10918) at quality 90 yields ~10:1 with PSNR ≈ 38 dB and 8×8 blocking artifacts rather than wavelet ringing.

Worked example D — Acoustic echo cancellation with NLMS

A speakerphone receives the far-end voice s_f[n], plays it through a speaker, and the local microphone picks up the talker’s own voice s_n[n] plus the echo of s_f[n] convolved with the room impulse response h[n] (typical reverberation 100–300 ms in a meeting room, ~12 000–36 000 taps at 16 kHz, or ~1500–4500 taps after acoustic decimation). The cancellation filter learns ĥ[n] adaptively.

Setup at f_s = 16 kHz, filter length N = 1024 taps (64 ms tail, enough for a small office):

// NLMS step per sample
float dot = 0.0f;
for (int i = 0; i < N; ++i) dot += w[i] * x_buf[i];     // filter output
float e = mic - dot;                                    // residual error
float energy = epsilon;
for (int i = 0; i < N; ++i) energy += x_buf[i] * x_buf[i];
float mu_eff = mu / energy;                             // normalised step
for (int i = 0; i < N; ++i) w[i] += mu_eff * e * x_buf[i];

Step size μ ≈ 0.5 (with ε ≈ 1e−6 to avoid division by zero when far-end is silent), giving misadjustment ~3 dB and convergence in ~5N samples ≈ 320 ms — slow but stable. To converge faster on speech, use frequency-domain partitioned block NLMS (overlap-save FFT, 8–16 partitions of N/8 taps each) — reduces compute by ~N/(log N) and converges per-partition. Double-talk detection (Geigel detector or coherence-based) freezes adaptation when both ends speak, otherwise the local talker’s voice corrupts ĥ. Post-filter (residual echo suppressor, Wiener gain) cleans up the 10–15 dB of residual after the linear filter converges. Codec-grade AEC ships in WebRTC (open source, Google), Speex AEC, and is built into every Zoom / Teams / Meet audio path.

Power spectral density via the Welch periodogram

Direct |DFT(x)|² is the periodogram: an asymptotically unbiased but high-variance estimator of the power spectral density S_xx(f) — the variance does not decrease with N. Welch (1967) averages the periodogram across K overlapping windowed segments:

S_xx[k] ≈ (1/K) · Σ_{i=0}^{K−1} |DFT(x_i · w)|² / (f_s · Σ|w|²)

with 50 % overlap and Hann window as defaults. Variance scales as 1/K_eff (effective independent segments after overlap correction). Trade record length T against bias (longer windows → finer Δf but fewer segments → higher variance). MATLAB pwelch(x, hann(N), N/2, N, fs); Python scipy.signal.welch(x, fs, nperseg=N, noverlap=N//2, window='hann'). Standard tool for vibration spectra, EEG band power, modulated-signal characterisation.

Variants: Bartlett periodogram (no overlap), Blackman–Tukey (lag-windowed autocorrelation FT, equivalent under appropriate windows), multitaper / Thomson DPSS (Slepian sequences as windows, optimal bias-variance trade), parametric AR / MA / ARMA (Yule–Walker, Burg, modified covariance — for short records where non-parametric methods lack resolution).

Real-input FFT optimisations and zoom-FFT

For real-input x[n] of length N, the DFT has conjugate symmetry X[N−k] = X*[k] — only N/2 + 1 bins carry independent information. Real-FFT algorithms exploit this: pack two real N-point inputs into one complex N-point FFT (split-via-symmetry), or use a length-N/2 complex FFT plus a post-processing “twist” step (FFTW r2c, CMSIS arm_rfft_fast_f32). Saves ~50 % memory and ~40 % compute vs naïve complex FFT.

Zoom-FFT computes a high-resolution spectrum on a narrow frequency band of interest without using a very long N. The recipe: heterodyne the input by multiplying with e^{−j·2π·f_center·n/f_s} (shifts the band to baseband), low-pass and decimate by M (sample rate drops to f_s/M), then FFT the decimated stream. The Δf = f_s/(M·N) is M× finer than a direct N-point FFT of the full-rate signal at the same compute. Used in vibration analysis (zoom on a known machinery fault frequency), audio analysis (instrument harmonics), RF channel sounding.

For non-power-of-2 lengths, Bluestein chirp-z transform computes a length-N DFT by zero-padding to length 2N+1 and convolving with a chirp sequence using a power-of-2 FFT — slower than radix-2 by a constant but supports prime and arbitrary N (used for example for 257-point or 1009-point DFTs in DRM digital radio).

Filter implementation cost on common targets

Target	f_clk	Single FIR MAC	8-tap FIR throughput	4th-order IIR throughput	1024-pt complex FFT
ATmega328P (AVR, no MAC, no FPU)	16 MHz	~30 cycles (Q15)	~2 µs	~3 µs	not feasible real-time
STM32F103 (Cortex-M3, 32-bit MUL)	72 MHz	4 cycles (Q15)	0.5 µs	0.6 µs	12 ms (CMSIS Q15)
STM32F4 (Cortex-M4F, FPU)	168 MHz	1 cycle (f32 SMLA)	0.05 µs	0.10 µs	0.8 ms (CMSIS f32)
STM32H7 (Cortex-M7F, dual FPU)	480 MHz	1 cycle (f32)	0.017 µs	0.04 µs	0.18 ms
TI TMS320C6678 (8× C66x VLIW DSP)	1 GHz/core	0.125 cycle (SIMD)	0.001 µs	0.005 µs	4 µs
ADI SHARC ADSP-21593	1 GHz	1 cycle (f32 SIMD)	0.008 µs	0.04 µs	8 µs
Xilinx UltraScale+ DSP48E2	891 MHz/slice	1 cycle (18×27 MAC)	pipeline-limited	sub-ns	100 ns (parallel)
NVIDIA H100 GPU (CUFFT)	—	—	—	—	0.6 µs (large batch)
Hailo-8 NPU	400 MHz	—	NN graph fused	—	not native

DSP-specific cores (C66x, SHARC, Hexagon) are losing market share to ARM Cortex-M7/M55 + Helium SIMD and to plain Cortex-A application processors with NEON, because the latter offer comparable per-watt throughput with mainstream toolchains. The classic standalone DSP survives in motor / inverter control (TI C2000 — see [[Engineering/digital-control]]), high-end audio (SHARC), and baseband modems (Qualcomm Hexagon).

Quantisation noise, dither, and noise shaping

An N-bit quantiser with full-scale range V_FS has step size q = V_FS / 2^N. The quantisation error e_q[n] = x_q[n] − x[n], for a busy input that exercises many levels, is well modelled as uniform white noise on [−q/2, +q/2] with variance σ²_q = q²/12. The SQNR for a full-scale sinusoid is:

SQNR = 10·log₁₀(P_signal / σ²_q) = 6.02·N + 1.76 dB

For a 24-bit ADC this is 146.2 dB nominal; the achieved ENOB is set by thermal noise, INL, DNL, and clock jitter, and is typically 20–23 bits in practice (see [[Engineering/digital-control]] for an extended ADC table). For low-amplitude tonal inputs the uniform-noise model breaks down — quantisation produces deterministic distortion harmonics, not noise. The fix is dither: add a small amount of white noise (typically TPDF — triangular probability density, two LSB peak-to-peak) before quantisation. This linearises the quantiser in the statistical sense — the total quantisation error becomes signal-independent at the cost of 3 dB worse noise floor. Mandatory in 16-bit audio mastering, professional mixing, and dB-of-the-LSB instrumentation.

Noise shaping moves quantisation energy out of the band of interest using a feedback loop around the quantiser (delta-sigma technique). A first-order noise-shaped quantiser has a noise transfer function NTF(z) = 1 − z^{−1} — high-pass shape, suppressing noise at DC by 3 dB for every doubling of oversampling ratio. Higher-order Σ-Δ ADCs (4th–8th order is common) push noise far above the band, then decimation filters remove it, leaving 18–24 bit effective resolution from a single-bit modulator. The Σ-Δ DAC in every USB audio interface and the Σ-Δ ADC in every smartphone microphone runs this principle.

Fixed-point vs floating-point

Q-format notation: Qm.n means m integer bits + n fractional bits in a signed two’s-complement word of m + n + 1 bits. Common: Q15 (1 sign + 15 fractional, ±1.0 range), Q31 (1 sign + 31 fractional), Q1.30 (audio coefficient with ±2 headroom).

Format	Word size	Dynamic range	Resolution	Typical use
Q15	16-bit	−1.0 to +0.999969	3.05·10⁻⁵	Audio sample data, FIR coeffs on M0/M3
Q31	32-bit	−1.0 to +0.999999999	4.66·10⁻¹⁰	High-precision filter coeffs, accumulators
Q1.30	32-bit	−2.0 to +2.0	9.31·10⁻¹⁰	Coefficients exceeding unity (e.g., FFT twiddles)
IEEE 754 binary32 (float)	32-bit	±3.4·10³⁸	24-bit mantissa (1 in 1.67·10⁷)	Default on Cortex-M4F/M7F, x86, DSP cores
IEEE 754 binary64 (double)	64-bit	±1.8·10³⁰⁸	53-bit mantissa	Reference / scientific / control of critical paths
bfloat16	16-bit	±3.4·10³⁸	8-bit mantissa (1 in 256)	ML inference; useless for general DSP

For sensor / audio / control DSP on modern silicon, single-precision float is the default — coefficient design and accumulator headroom are trivial. Fixed-point survives in cost-sensitive volume products (USD < 0.50 BOM MCUs without FPU), in FPGA fabric (every bit costs LUTs), and in NN inference where INT8 / INT16 quantisation halves silicon area and power. ARM CMSIS-DSP, Intel IPP, and TI DSPLIB all provide both fixed and float APIs.

Hardware tier table

Tier	Representative parts	f_s reach	Use case
Microcontroller (no FPU)	ATmega328P, STM32F0, RP2040	up to ~50 kHz audio (Q15)	Sensor LP + decimation; simple loops
MCU + FPU + DSP ext	STM32F4/F7/H7, NXP i.MX RT, RP2350	up to ~1 MHz signal	Mainstream audio + control
MCU + SIMD (Helium / DSP)	Cortex-M55, Cortex-M85, dsPIC33	up to ~5 MHz	NN inference, baseband, ANC
Dedicated DSP	TI C2000, TI C66x, ADI SHARC, Qualcomm Hexagon	up to ~10 MHz	Motor / inverter, modem, high-end audio
FPGA	Xilinx Spartan/Artix/Kintex, Intel Cyclone/Stratix, Lattice ECP5/Avant, Microchip PolarFire	up to ~1 GHz	Multi-channel parallel, sub-µs latency
FPGA + RFSoC	Xilinx Versal RFSoC, Intel Agilex F-series	direct RF (6+ GS/s ADC/DAC)	5G L1, software-defined radar, mil-SIGINT
GPU	NVIDIA RTX / H100, AMD MI300, Apple M-series	very large transforms; latency ~ms	Beamforming, MRI recon, ML preprocessing
NPU / VPU	Hailo-8, Coral Edge TPU, Movidius Myriad X	NN graph rate	Camera + audio NN pipelines

5p. Theory — derivations

Cooley–Tukey radix-2 FFT (1965)

For N = 2^k, the DFT splits into even- and odd-indexed sub-DFTs:

X[k] = Σ_{n=0}^{N/2−1} x[2n] · W_N^{2nk} + W_N^k · Σ_{n=0}^{N/2−1} x[2n+1] · W_N^{2nk}

where W_N = e^{−j·2π/N}. Both sums are themselves N/2-point DFTs evaluated at frequency 2k (which, mod N/2, is k). Recursion gives the algorithm: N/2 butterfly operations per stage, log₂(N) stages, total ≈ (N/2)·log₂(N) complex multiplies vs N² for the direct DFT. For N = 1024, ratio is 1024²/(512·10) = 205× speedup.

Twiddle factors W_N^k repeat with conjugate symmetry; radix-4 exploits W_N^{N/4} = −j so a stage costs zero real multiplies for one of its four input combinations. Split-radix (Yavne 1968, Duhamel & Hollmann 1984) blends radix-2 and radix-4 to minimise real multiplies, with ≈ N·log₂(N) − 3N + 4 real multiplies for real-input FFT.

Numerical considerations: floating-point FFT accumulates rounding error of order O(sqrt(log N)) per input; fixed-point requires bit-growth of 1 bit per radix-2 stage in worst-case scaling, or saturation arithmetic with per-stage block-floating-point scaling.

Bilinear transform — frequency warping derivation

The trapezoidal-rule approximation of integration ∫_0^{T_s} u dt ≈ (T_s/2)·(u[n] + u[n−1]) maps the integrator 1/s to (T_s/2)·(z+1)/(z−1), equivalently s ← (2/T_s)·(z−1)/(z+1). Substituting z = e^{jωT_s} into the inverse:

s = (2/T_s)·(e^{jωT_s} − 1)/(e^{jωT_s} + 1) = j·(2/T_s)·tan(ωT_s/2)

The analog frequency Ω corresponding to a discrete frequency ω is Ω = (2/T_s)·tan(ω·T_s/2). For ω·T_s ≪ 1, Ω ≈ ω; near Nyquist (ω → π/T_s), Ω → ∞ — the entire analog frequency axis is compressed into [0, π/T_s] of the digital axis. Pre-warping: when discretising an analog filter with a critical frequency Ω_c (e.g. notch centre, cutoff), substitute s ← α·(z − 1)/(z + 1) where α = Ω_c/tan(Ω_c·T_s/2) instead of α = 2/T_s. This makes the bilinear mapping exact at Ω_c.

LMS convergence and misadjustment bounds

For an N-tap LMS filter with input vector x[n], the autocorrelation matrix R = E{x·x^T} has eigenvalues {λ_i}. Mean convergence of the coefficient error e_w[n] = w[n] − w_opt requires 0 < μ < 2/λ_max. The convergence time constant of the i-th mode is τ_i ≈ 1/(μ·λ_i) samples — fastest along the largest-eigenvalue eigenvector, slowest along the smallest. The eigenvalue spread χ(R) = λ_max/λ_min sets the worst-to-best convergence ratio; for white noise input χ = 1 and LMS converges fast; for narrowband or correlated inputs χ can exceed 1000 and convergence stalls. Misadjustment (steady-state excess MSE relative to Wiener) M ≈ μ·N·P_x/2 — the trade is μ small for low misadjustment but slow convergence. NLMS divides μ by the instantaneous input power so the effective step adapts to signal level; misadjustment is approximately μ·N/2 independent of input power, and convergence robust to amplitude swings.

For non-stationary inputs and tracking time-varying channels, the excess MSE from tracking lag must be added: M_track ≈ (N · σ²_h_dot · ε) / μ where σ²_h_dot is the variance of the rate-of-change of the true filter. Minimum total MSE comes from balancing M_adaptation against M_tracking — there is an optimum μ for each channel rate-of-change.

Daubechies wavelets — orthogonality + compact support

The 1988 result (Daubechies, Comm. Pure Appl. Math.) constructs an orthonormal basis of L²(ℝ) of compactly supported wavelets {ψ_{j,k}(t) = 2^{j/2}·ψ(2^j·t − k)} via a low-pass filter h[n] of length 2N (the Daubechies-N filter) whose Z-transform H(z) satisfies:

• Quadrature mirror condition: |H(e^{jω})|² + |H(e^{j(ω+π)})|² = 2 (orthogonality)
• N vanishing moments: H(z) has a factor (1 + z^{−1})^N (smoothness)
• Minimum-phase / compact support: filter length 2N

The high-pass mirror is g[n] = (−1)^n · h[2N − 1 − n]. Discrete wavelet transform is then the recursive application of [↓2, g; ↓2, h] (analysis filter bank) — log₂(N) stages on an N-sample signal. Compactly supported wavelets give finite-length FIR filter banks, making the DWT O(N) total — faster than the FFT.

6p. Application — specialised domains

Audio

• Parametric equaliser: cascade of biquads with shelf, peak, low-pass, high-pass types per band; coefficients computed from RBJ Audio EQ Cookbook (Bristow-Johnson 1994, widely circulated reference, used in every DAW EQ plugin).
• Dynamics processors (compressor, limiter, gate): envelope detector (full-wave rectify + log + 1st-order LP) → gain computer → smoothed gain → multiplier. Look-ahead limiters delay the signal by 5–20 ms to react before the peak arrives.
• Reverb: Schroeder reverberator (1962) — 4–8 parallel comb filters + 2–4 series all-pass filters; modern algorithmic reverbs (FDN — Feedback Delay Network) extend to 8–32 channels for natural decay.
• Active Noise Cancellation (ANC): feedforward (microphone → LMS adaptive FIR → anti-noise speaker), feedback (error mic → adaptive filter), hybrid (both); typical filter 64–512 taps at 48 kHz; NLMS step adapted to input power.
• Pitch shifting / time stretching: phase vocoder (STFT → phase manipulation → ISTFT with overlap-add), PSOLA (pitch-synchronous overlap-add), neural (DiffSinger, RVC voice conversion).
• Speech codecs: Opus (RFC 6716, 6 kbps–510 kbps, SILK + CELT hybrid), AMR-WB (G.722.2), LC3 (Bluetooth LE Audio), codec2 (David Rowe, open source for ham radio, 1.2–3.2 kbps), Lyra v2 (Google, neural, 3.2 kbps).

Communications

• OFDM (Orthogonal Frequency Division Multiplexing): parallel narrowband subcarriers via inverse FFT; cyclic prefix absorbs multipath; equaliser per-subcarrier is one complex multiply. Used in Wi-Fi, LTE, 5G NR, DVB-T, DOCSIS, DSL.
• Channel equalisation: zero-forcing, MMSE linear, decision-feedback (DFE), turbo (LDPC + soft equaliser); adaptive via LMS/RLS for time-varying channels (mobile).
• Matched filter: correlate against the known transmit waveform to maximise SNR before symbol decision. Foundation of radar pulse compression, GPS code acquisition, DSSS.
• Forward error correction: convolutional + Viterbi (Viterbi 1967), turbo codes (Berrou 1993), LDPC (Gallager 1962, rediscovered 1996), polar codes (Arıkan 2009, 5G control channel).

Video and image

• JPEG (ISO/IEC 10918, 1992): 8×8 DCT-II per block, quantisation matrix, zig-zag scan, Huffman or arithmetic coding.
• H.264/AVC (ISO/IEC 14496-10, 2003): 4×4 / 8×8 integer DCT, CAVLC/CABAC entropy, intra-prediction, motion estimation, deblocking filter.
• H.265/HEVC (ISO/IEC 23008-2, 2013): variable transform sizes 4–32 with DCT-II and DST-VII, larger CTUs.
• AV1 (AOMedia, 2018): open-source successor; 64×64 superblocks, multiple transform types, neural-style intra prediction.
• JPEG2000 (ISO/IEC 15444, 2000): 5/3 lossless or 9/7 lossy DWT, EBCOT block coder; cinema (DCI digital cinema package), medical archival.

Industrial / sensor

• Vibration condition monitoring: FFT of accelerometer at 25.6 kHz for rotating machinery; envelope detection (Hilbert + magnitude) reveals bearing fault frequencies (BPFO, BPFI, BSF, FTF) that are hidden in raw vibration. See [[Engineering/vibration-dynamics]].
• Motor current signature analysis (MCSA): FFT of motor current finds broken rotor bars (sidebands around mains frequency at ±2·s·f_mains where s is slip), eccentricity, bearing wear — all without a dedicated sensor ([[Engineering/electric-motors]]).
• Ultrasonic flow meter: cross-correlate upstream and downstream transit times; sub-ns resolution via interpolation.
• Lock-in amplifier: modulate signal at f_ref, multiply received by sin(2π·f_ref·t) and cos(…), low-pass aggressively; extracts narrowband signal buried below the noise floor by 60+ dB.

Machine learning preprocessing

Front-ends for ML on audio, vibration, and biosignal data are almost always classical DSP — the network learns from features, not raw samples, because raw-sample networks need far more data and parameters to discover what a 50-line DSP block computes deterministically.

• Mel-spectrogram (audio): STFT (e.g., 25 ms window, 10 ms hop at 16 kHz → 400-sample frame, 160-sample hop, 512-point FFT) → magnitude squared → mel filterbank (40–128 triangular filters spaced on mel scale, 0–8 kHz) → log. Feeds keyword spotting (Hello Google / Alexa / Siri wake-word), speech recognition (Conformer, Whisper), speaker ID, audio event detection. Computed natively on Tensilica HiFi DSP, Hexagon, Cortex-M55 + CMSIS-NN; ~50 µs per 25-ms frame on Cortex-M55 at 480 MHz.
• MFCC (Mel-Frequency Cepstral Coefficients): mel-spectrogram → DCT-II → keep first 12–20 coefficients. Older speech-recognition standard (HMM-GMM era); still used in compact wake-word models.
• Log-mel features for vibration / acoustic monitoring: identical pipeline at different parameters (4 kHz sample rate, 32 mel bins, longer windows) — bearing fault classification, valve leak detection, predictive maintenance.
• PCG (constant-Q spectrogram): logarithmic frequency spacing aligned with musical octaves, used in music information retrieval and tonal analysis.
• Whisper / Conformer / Wav2Vec2: modern speech models still consume STFT or mel-spectrogram inputs even though end-to-end raw-audio variants exist — the inductive bias from DSP preprocessing remains beneficial at practical scales.

The growing field of differentiable DSP (DDSP, Engel et al. 2020) flips the relationship: filter coefficients, oscillator parameters, and modulators are learned via backprop through analytic DSP primitives. The result combines the data-efficiency of structured signal models with the flexibility of neural training; applied to instrument synthesis, voice conversion, and audio effect modelling.

Medical

• ECG R-peak detection (Pan & Tompkins 1985): band-pass 5–15 Hz → differentiate → square → moving-window integrate (150 ms) → threshold with adaptive bias. Detects 99.3 % of beats on MIT-BIH benchmark; baseline of every clinical ECG monitor (IEC 60601-2-25). Modern replacements use 1-D CNNs.
• EEG quantification: band-power in δ (0.5–4 Hz), θ (4–8), α (8–13), β (13–30), γ (30–100); computed by Welch periodogram with 50 % overlap and Hann windows.
• MRI reconstruction: raw k-space samples are 2-D / 3-D FFTed into image space; modern accelerated MRI uses non-uniform FFT + compressed sensing (Lustig 2007).
• Ultrasound beamforming: delay-and-sum across N transducer elements; modern systems use plane-wave imaging with software delay-and-sum on GPU at 10 000 frames/s.

7p. Edge cases and gotchas

Forgotten anti-alias filter. A 1 kHz tone at the input of a 200 Hz-sampled accelerometer logger appears as a 200 Hz alias if no AAF is fitted. Symptom: ghost peaks in the spectrum at frequencies that don’t correspond to any real mechanical mode. Fix: passive RC ahead of the ADC pin sized for f_−3dB ≈ f_s/4, or active filter for sharper rolloff ([[Engineering/op-amps]]).

Spectral leakage masking weak signals. A 0 dBFS tone at a non-integer-bin frequency, sampled with a rectangular window, has side-lobes at −13 dB falling at only −6 dB/octave. A 1100 Hz weak signal 60 dB below a strong 1000 Hz tone is utterly hidden under rectangular-window leakage. Fix: window the data (Hann minimum, Blackman-Harris for spectrum analysers), or arrange synchronous sampling so all tones sit on bin centres.

IIR limit cycles in fixed-point. A 2nd-order section with poles very near the unit circle (high-Q resonator, narrow notch, slow integrator) can sustain a low-amplitude oscillation at the natural frequency forever after a transient, because the quantisation grid doesn’t include zero on the resonant orbit. Fix: dither, increase wordlength, use floating-point, or use a structure that explicitly preserves passivity (lattice / wave-digital).

Coefficient quantisation moving a pole outside the unit circle. A direct-form 8th-order IIR designed to have a pole at z = 0.998 + 0.063j (just inside the unit disc) can land at z = 1.001 + … when coefficients are quantised to Q15 — and the filter blows up on first input. Fix: factor into biquads (scipy.signal.sos); the poles of each biquad section move independently and modestly, and stability of each section is independently verifiable.

Accumulator overflow in FIR. A 128-tap FIR with Q15 coefficients summing 128 products of Q15 values needs ⌈log₂(128)⌉ = 7 guard bits of headroom above Q15 in the accumulator — a 16-bit accumulator overflows on the worst case. Use a Q1.30 or Q1.62 accumulator and saturate to Q15 on output. Cortex-M4 SMLAD instruction provides a 32-bit accumulator natively; CMSIS-DSP exploits.

DC offset and very-low-frequency drift. Sensor signals often carry a DC bias and slow drift (thermocouple cold-junction, op-amp Vos, electrode polarisation in ECG). A naïve FFT puts most energy in bin 0 and bin 1, masking signal. Fix: DC-blocking IIR (single real pole, H(z) = (1 − z^{−1})/(1 − α·z^{−1}), α ≈ 0.995 — corner ~10 Hz at 48 kHz) or running-mean subtraction.

Group-delay mismatch in cascade. Cascading two IIR filters with different group delays distorts transient signals (square waves, click impulses) — the high-frequency content arrives before the low-frequency. Fix: design as linear-phase FIR if phase matters, or apply a forward-backward (filtfilt) for zero-phase response in offline / non-real-time work.

FFT zero-padding misinterpreted as resolution increase. Zero-padding an N-sample record to 2N doesn’t add information — it interpolates the existing N/2 DTFT samples to 2N samples. Two unresolved tones in the N-sample record remain unresolved. Real frequency resolution comes only from longer record length (more time-domain samples of the actual signal).

Block-vs-sample real-time latency. A 1024-sample FFT at 48 kHz introduces 21.3 ms of buffering delay before the first output appears — fatal for live ANC, telephony, monitoring headphones (where 10 ms is the perceptual limit). Sample-by-sample processing or short overlapping blocks (32–128 samples) are required for low-latency work; long blocks suit batch / mastering / non-real-time analysis.

Real-time multitasking and sample jitter. A higher-priority ISR delaying the audio ISR introduces sample-time jitter; the FFT then has a frequency-modulation effect (sidebands at the jitter rate). Symptom: clean tones grow ±50 Hz sidebands around their fundamental when the system is loaded. Fix: priority the audio/DSP ISR above all else, use DMA double-buffering so the CPU never blocks waiting for the converter, instrument with GPIO scope toggles.

Endianness in cross-platform sample storage. WAV files are little-endian, AIFF big-endian; embedded codecs may stream MSB-first or LSB-first; FFT results are platform-native unless converted. A byte-order mismatch produces sample noise that is white-noise-like and easy to mistake for ADC noise. Fix: validate file headers, use portable serialisation (htons / htonl family, or explicit struct.pack('<i', ...) in Python).

Wavelet edge effects. DWT of a finite-length signal is implicitly periodic — the filter wraps around the ends, producing artifacts at signal boundaries. Mitigations: symmetric extension (pywt.Modes.symmetric), zero-padding, smooth reflection. Lossy compression standards (JPEG2000) use symmetric extension explicitly.

Decimation without anti-alias filter. Downsampling a digital signal by M without first low-pass filtering folds energy above f_s/(2M) back as alias — same phenomenon as un-AAF’d ADC sampling. Always pair decimation with a low-pass filter, even when “the input is already bandlimited” — most signals carry quantisation or processing noise that aliases.

Cycle budget for worst-case path. A DSP loop nominally running at 48 kHz with 80 % CPU load fails when an occasional control-plane task (USB enumeration, network ISR, OS housekeeping) takes the remaining 20 %. Real-time scheduling needs worst-case execution time (WCET) analysis, not average. Tools: aiT (AbsInt), Bound-T, manual loop profiling on the target with the DWT cycle counter.

8p. Tools and software

The dominant practical flow has three layers — algorithm prototyping, target-specific implementation, and verification — that engineering teams stitch together differently depending on industry.

Algorithm prototyping in 2026 is almost exclusively MATLAB or Python. MATLAB Signal Processing Toolbox + DSP System Toolbox is the industry default for established players (defence primes, automotive Tier 1s, telecom equipment vendors) — fir1, firpm, firpmord, butter, cheby1, ellip, besself, bilinear, sosfilt, zp2sos, freqz, fvtool, dspfpga codegen, fixed-point conversion tools, and HDL Coder for FPGA generation. Python (scipy.signal, numpy.fft, pyfftw for FFTW bindings, pywt for wavelets, librosa and torchaudio for audio, soundfile, python-sounddevice for I/O) has overtaken MATLAB in startups, academia, and ML-adjacent work. PyTorch / JAX enable differentiable DSP — filter coefficients and modulator parameters can be optimised by gradient descent (Engel et al. 2020, “DDSP: Differentiable Digital Signal Processing”). GNU Octave is a free MATLAB-compatible interpreter for budget projects. GNU Radio is a flowgraph-based DSP environment with hardware integration for software-defined radio — Ettus USRP, HackRF, LimeSDR, BladeRF — used in research, ham radio, and commercial SDR.

Target implementation depends on the chip family. ARM CMSIS-DSP (arm_fir_*, arm_biquad_*, arm_rfft_fast_*, arm_lms_*) covers Cortex-M0/M3/M4/M7/M55/M85 with hand-optimised assembly per architecture; the Helium SIMD path on M55/M85 gives 4× throughput on bulk FIR. TI DSPLIB / IMGLIB / VLIB target C2000 / C6000 with intrinsics and assembly kernels; TI Mathlib handles transcendentals. Analog Devices VisualDSP++ libraries for SHARC and Blackfin. Intel IPP (Integrated Performance Primitives) for x86 (auto-vectorises across SSE/AVX/AVX-512); Intel MKL for FFT and linear algebra. FFTW (Frigo & Johnson, GPL + commercial license from MIT) is the canonical portable FFT library; chooses optimal radix and code path at runtime via its planning step; benchmark-leader for portable code on x86 / ARM / POWER. CuFFT (NVIDIA) for GPU; rocFFT (AMD); MKL DFTI (Intel hardware-accelerated). For FPGA, Xilinx Vitis HLS and Vitis DSP Library for high-level synthesis of DSP kernels (FFT, FIR, polyphase) into RTL targeting DSP48E2 / DSP58 slices.

Audio real-time I/O layers: PortAudio (cross-platform C library), JACK (low-latency Linux/macOS pro audio), ALSA (Linux kernel-level), CoreAudio (macOS / iOS), WASAPI / ASIO (Windows), AAudio / OpenSL ES (Android). VST3 and Audio Unit (AU) are the plugin formats hosting DSP code inside DAWs (Pro Tools, Cubase, Logic, Ableton Live, Reaper, Bitwig, FL Studio). JUCE is the dominant cross-platform plugin framework — every commercial VST/AU shipped in the last decade likely uses it.

Test and measurement instrumentation: Audio Precision APx555 / APx500-series (audio analysis, distortion, frequency response, jitter — the de facto golden audio instrument), Brüel & Kjær PULSE LabShop (vibration and acoustic analysis with FFT, octave, modal), National Instruments LabVIEW + Sound and Vibration Toolkit + PXIe FPGA modules, Keysight 89600 VSA (vector signal analyser for RF DSP characterisation), Rohde & Schwarz FSW / SMW (modulation analysis), Tektronix RSA series (real-time spectrum analyser).

DSP development boards for embedded prototyping: TI LaunchPad family (F2837x, F28004x for motor/power; C66x EVMs for high-performance DSP), NXP i.MX RT1170-EVK (Cortex-M7 @ 1 GHz), STM32 Nucleo-H723/H743, STM32H7B3I-DK with Cortex-M7 + Chrom-ART, Analog Devices EZ-KIT for SHARC and Tigershark, Xilinx Zynq UltraScale+ ZCU102 / ZCU111 RFSoC, Intel Cyclone V SoC Devkit, Lattice CrossLink-NX evaluation kit, Hailo-8 M.2 for NN inference.

Static analysis and verification for safety-critical DSP code: LDRA Testbed, MathWorks Polyspace, Synopsys Coverity, AbsInt aiT (WCET analysis). MISRA C:2012 and CERT C compliance is standard for IEC 60601 medical, ISO 26262 automotive, IEC 61508 industrial, and DO-178C / DO-254 aerospace.

11. Cross-references

• [[Engineering/digital-control]] — same sampling theory, z-transform, and discretisation but oriented to closed-loop stability; pre-warping and ZOH discussed in both.
• [[Engineering/digital-logic]] — building blocks for FPGA / ASIC DSP datapaths; flip-flop timing, CDC.
• [[Engineering/op-amps]] — analog front-end, anti-alias filters, instrumentation amplifiers ahead of every ADC.
• [[Engineering/fpga-design]] — DSP slice (DSP48E2 / DSP58 / Variable-Precision DSP) implementation of FFT, polyphase filter banks, beamformers.
• [[Engineering/microcontrollers]] — substrate for embedded DSP; Cortex-M4F/M7F/M55, Helium SIMD, DMA double-buffering, ADC peripherals.
• [[Engineering/realtime-embedded]] — RTOS, ISR priorities, jitter and worst-case execution time, double-buffer audio.
• [[Engineering/electromagnetics-engineering]] — RF + antenna signal characteristics, source of communications-DSP inputs.
• [[Engineering/vibration-dynamics]] — FFT and envelope analysis of structural and rotating-machine signals; modal analysis.
• [[Engineering/electric-motors]] — motor current signature analysis, FOC sensor pipelines.
• [[Engineering/state-space-methods]] — Kalman filter, MIMO observers; companion linear-estimation methodology.
• [[Engineering/semiconductor-devices]] — ADC and DAC physics underlying every DSP input/output stage.
• [[Engineering/power-electronics]] — switching-converter ripple analysis, conducted-EMI spectrum, synchronous PWM sampling.
• planned [[Engineering/rf-design]] — RF front-end, IQ modulation, image rejection.
• planned [[Engineering/system-identification]] — fitting LTI and nonlinear models from data; companion of adaptive filtering.
• planned [[Engineering/rf-design]] — OFDM, channel codes, MIMO, end-to-end link design.
• planned [[Robotics/bayesian-estimation]] — Kalman / EKF / UKF / particle filters running at sample rate.
• [[Languages/Tier3/scientific]] — toolchain for DSP algorithm design and codegen (when promoted from Tier3).
• planned [[Languages/Tier3/scientific]] — scipy.signal, numpy.fft, pywt.

12. Citations

• Oppenheim, A. V. & Schafer, R. W. (2010). Discrete-Time Signal Processing (3rd ed.). Pearson. The canonical graduate textbook; DFT, FFT, sampling, multirate, parametric signal modelling, spectral estimation.
• Proakis, J. G. & Manolakis, D. K. (2007). Digital Signal Processing: Principles, Algorithms, and Applications (4th ed.). Pearson. Undergraduate-to-graduate; strong FIR / IIR design coverage and adaptive filters.
• Lyons, R. G. (2010). Understanding Digital Signal Processing (3rd ed.). Prentice Hall. Practical engineering treatment; tricks-of-the-trade, intuitive geometry of DFT.
• Mitra, S. K. (2010). Digital Signal Processing: A Computer-Based Approach (4th ed.). McGraw-Hill. MATLAB-centred treatment with worked exercises.
• Smith, S. W. (1997). The Scientist and Engineer’s Guide to Digital Signal Processing. California Technical Publishing. Free at dspguide.com; intuitive entry point, broad coverage.
• Mallat, S. (2008). A Wavelet Tour of Signal Processing: The Sparse Way (3rd ed.). Academic Press. Modern wavelet bible; multi-resolution analysis, sparse representations.
• Daubechies, I. (1992). Ten Lectures on Wavelets. SIAM CBMS-NSF 61. Foundational; defines the compactly supported orthonormal wavelet families.
• Vaidyanathan, P. P. (1993). Multirate Systems and Filter Banks. Prentice Hall. Definitive treatment of polyphase decomposition, perfect-reconstruction filter banks, wavelet–filter-bank duality.
• Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. LMS, RLS, Kalman, neural adaptive filtering.
• Widrow, B. & Stearns, S. D. (1985). Adaptive Signal Processing. Prentice Hall. Original LMS work consolidated.
• Cooley, J. W. & Tukey, J. W. (1965). “An Algorithm for the Machine Calculation of Complex Fourier Series.” Mathematics of Computation, 19(90), 297–301. DOI: 10.1090/S0025-5718-1965-0178586-1. The radix-2 FFT.
• Parks, T. W. & McClellan, J. H. (1972). “Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase.” IEEE Transactions on Circuit Theory, 19(2), 189–194. Optimal equiripple FIR design via Remez exchange.
• Kaiser, J. F. (1974). “Nonrecursive Digital Filter Design Using the I_0–sinh Window Function.” Proc. IEEE Int. Symp. Circuits and Systems. Kaiser window family and FIR order estimate.
• Widrow, B. & Hoff, M. E. (1960). “Adaptive Switching Circuits.” IRE WESCON Convention Record, Part 4, 96–104. LMS algorithm.
• Daubechies, I. (1988). “Orthonormal Bases of Compactly Supported Wavelets.” Communications on Pure and Applied Mathematics, 41(7), 909–996. DOI: 10.1002/cpa.3160410705.
• Ahmed, N., Natarajan, T. & Rao, K. R. (1974). “Discrete Cosine Transform.” IEEE Transactions on Computers, C-23(1), 90–93. DCT-II.
• Pan, J. & Tompkins, W. J. (1985). “A Real-Time QRS Detection Algorithm.” IEEE Transactions on Biomedical Engineering, BME-32(3), 230–236. R-peak detector still embedded in clinical monitors.
• Welch, P. D. (1967). “The Use of Fast Fourier Transform for the Estimation of Power Spectra.” IEEE Transactions on Audio and Electroacoustics, 15(2), 70–73. Welch periodogram.
• Schroeder, M. R. (1962). “Natural Sounding Artificial Reverberation.” Journal of the Audio Engineering Society, 10(3), 219–223. Algorithmic reverberation.
• Bristow-Johnson, R. (1994). “Audio EQ Cookbook.” Widely circulated technical note; biquad coefficient formulas for shelving and peaking EQ.
• Frigo, M. & Johnson, S. G. (2005). “The Design and Implementation of FFTW3.” Proceedings of the IEEE, 93(2), 216–231. Adaptive FFT planner.
• Engel, J., Hantrakul, L., Gu, C. & Roberts, A. (2020). “DDSP: Differentiable Digital Signal Processing.” ICLR 2020. arXiv:2001.04643. Combines DSP primitives with neural learning.
• Lustig, M., Donoho, D. & Pauly, J. M. (2007). “Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging.” Magnetic Resonance in Medicine, 58(6), 1182–1195.
• IEEE Std 754-2019. IEEE Standard for Floating-Point Arithmetic. Defines binary32, binary64, bfloat16, rounding modes, special values.
• IEC 60601-2-25:2011 + AMD1:2020. Medical electrical equipment — Part 2-25: Particular requirements for the basic safety and essential performance of electrocardiographs. Bandwidth, noise, R-peak detection accuracy.
• ITU-T Recommendation G.722 (2012). 7 kHz audio-coding within 64 kbit/s. Sub-band ADPCM speech codec.
• IETF RFC 6716 (2012). Definition of the Opus Audio Codec. Hybrid SILK + CELT codec.
• ISO/IEC 14496-10:2020 (H.264/AVC). Information technology — Coding of audio-visual objects — Part 10: Advanced Video Coding.
• ISO/IEC 23008-2:2020 (H.265/HEVC). Information technology — High efficiency coding and media delivery in heterogeneous environments — Part 2: High efficiency video coding.
• ISO/IEC 15444-1:2019 (JPEG2000 Core). Information technology — JPEG 2000 image coding system: Core coding system.
• ARM (2023). CMSIS-DSP Software Library Documentation. ARM Limited. Reference for arm_fir_*, arm_biquad_cascade_df2T_*, arm_rfft_fast_*, arm_lms_*, arm_cfft_* and the Helium SIMD path.
• MATLAB Signal Processing Toolbox and DSP System Toolbox documentation (R2025b).
• SciPy v1.14 documentation, scipy.signal and scipy.fft modules.
• PyWavelets (pywt) v1.6 documentation.
• FFTW 3.3.10 documentation (Frigo & Johnson, MIT).