Digital filter design: FIR, IIR and ultrasound filters
An ultrasound chain is a cascade of filters: the RF band-pass after the ADC, the decimation and IQ low-pass, the Hilbert FIR, the Doppler wall filter on 8–16 slow-time samples, the smoothing on parametric maps. This module turns filter design into geometry (poles and zeros), into order formulas you can size from memory, and into the failure modes that matter for imaging: dispersion, edge transients, short-ensemble transients and fragile coefficients.
You will be able to
- Read magnitude, phase and group delay off a pole-zero plot, and place poles and zeros to shape a response
- Size an FIR from a spec in your head (Kaiser window, equiripple estimate, the A/22 rule) and check it
- Design Butterworth, Chebyshev and elliptic filters by the bilinear transform with prewarping, and ship them as biquads
- Predict the delay and the pulse distortion of a design, and say when zero-phase filtering is allowed
- Build decimation chains (half-band, polyphase, CIC) and a Hilbert FIR
- Choose and defend a Doppler wall filter for an 8–16-sample ensemble
Assumes
- Module 03 · Convolution, LTI systems and filters: difference equations, \(H(z)\), poles and zeros, stability, phase and group delay, the windowing theorem
- Module 04 · Sampling, Fourier & LTI systems: sampling, windows and their figures of merit, which filter sits where on the ultrasound signal chain
- Module 06 · Eigendecomposition & SVD: SVD clutter filtering, the data-adaptive wall filter
- Module 08 · Analytic signal & IQ demodulation: IQ demodulation, decimation, the analytic signal
M is the FIR order: M + 1 taps, linear-phase delay M/2 samples (Module 03 counts taps, so its (M − 1)/2 is the same number). \(\omega = 2\pi f/f_s\) in rad/sample; \(\Delta f\) is a transition width normalized to \(f_s\) and \(\Delta\omega = 2\pi\Delta f\). \(A = -20\log_{10}\delta\) in dB (amplitude). For IIR specs, \(A_p\) is the peak-to-peak passband ripple and \(A_s\) the minimum stopband attenuation, both in dB. Default sampling rate for RF examples: \(f_s\) = 40 MHz.
Fast-time filters run on every channel before beamforming and must add the same known delay everywhere (Sections 09.2–09.6); slow-time wall filters run per pixel on 8–16 samples and are judged as matrices (Section 09.7). The order-72 band-pass is the minimum equiripple design for the 3–7 MHz spec of Demo 09.4.
09.1The z-plane: the response is geometry
Recap from Module 03. Section 03.7 derives the transfer function \(H(z) = B(z)/A(z)\) from the difference equation (03.21), the vector picture of the response (03.23), and the stability condition (03.24). Factored as \(H(z) = g\prod_k (1 - z_k z^{-1})/\prod_k (1 - p_k z^{-1})\), the magnitude on the unit circle is the product of distances from \(e^{j\omega}\) to the zeros over the product of distances to the poles, and the phase is the sum of the corresponding angles. A pole close to the circle makes a peak where the circle passes it; a zero on the circle makes an exact null. This module uses that geometry as a design tool: where to put roots to meet a spec, and what each root costs in delay.
Group delay root by root. \(\tau_g(\omega) = -d\varphi/d\omega\) is additive over factors. For a root at \(r e^{j\theta}\):
- \(\omega\)
- digital frequency, rad/sample
- \(r,\ \theta\)
- radius and angle of the root (zero or pole), \(r\) dimensionless, \(\theta\) rad
- \(\tau_z,\ \tau_p\)
- group delay of one zero and of one pole, samples
- \(\cos\)
- cosine of an angle in rad
| Property | Pole-zero condition | Consequence |
|---|---|---|
| Causal and stable | All poles inside |z| = 1 | Impulse response decays; an FIR (all poles at z = 0) is always stable |
| Real coefficients | Roots real or in conjugate pairs | |H| symmetric about ω = 0: the plots below show 0 to π only |
| Minimum phase | Zeros also inside | Least delay for a given |H|; energy front-loaded; stable inverse |
| All-pass | Each pole p paired with a zero at \(1/p^*\) | |H| = 1, only delay; used for phase equalization |
| Linear phase (FIR) | Zeros in reciprocal pairs \((z, 1/z^*)\) or on the circle | Constant delay M/2; cannot be minimum phase unless trivial |
Drag poles and zeros, watch the response
Exact formulaMagnitude, phase and group delay evaluated exactly from the roots by (03.23) and (09.1); impulse response by the difference equation; JS.Slow-time filters run at the PRF, so pole radii translate into pulses. At PRF = 4 kHz a pole at r = 0.95 rings for about 20 pulses (5 ms), longer than an 8–16-pulse color-flow ensemble. That single number explains why IIR wall filters never reach steady state in color flow (Section 09.7), while the same filter on a 4000-sample RF line at 40 MHz settles in 0.5 µs.
09.2FIR design by windowing: the Kaiser formulas
Start from the ideal response, delay it by M/2 so it can be causal, truncate to M + 1 taps, and taper the truncation with a window \(w[n]\):
- \(h[n]\)
- FIR impulse response, \(n = 0,\dots,M\)
- \(h_d[n]\)
- ideal low-pass impulse response, delayed by \(M/2\)
- \(w[n],\ W\)
- window and its DTFT
- \(n\)
- sample index
- \(M\)
- filter order (\(M + 1\) taps), even
- \(\omega_c\)
- cutoff, rad/sample
- \(H(e^{j\omega}),\ H_d\)
- designed and ideal frequency responses
- \(\circledast\)
- circular convolution in \(\omega\)
- \(\sin\)
- sine of an angle in rad
- \(w[n]\)
- Kaiser window, \(n = 0,\dots,M\)
- \(I_0\)
- zeroth-order modified Bessel function of the first kind
- \(\beta\)
- Kaiser shape parameter
- \(A\)
- attenuation \(-20\log_{10}\delta\), dB, with \(\delta = \min(\delta_p, \delta_s)\)
- \(\delta_p,\ \delta_s\)
- passband and stopband ripple (linear amplitude)
- \(M\)
- filter order (\(M + 1\) taps)
- \(n\)
- sample index
- \(\Delta\omega\)
- transition width \(\omega_s - \omega_p\), rad/sample
- \(\omega_p,\ \omega_s,\ \omega_c\)
- passband edge, stopband edge, cutoff, rad/sample
| Window | Window peak sidelobe | Filter min. stopband | Passband ripple | Transition width | Kaiser equivalent (beta, width) |
|---|---|---|---|---|---|
| Rectangular | −13.3 dB | −21.0 dB | 0.090 | 1.85π/M | 0 |
| Hann | −31.5 dB | −43.9 dB | 0.0064 | 6.25π/M | 3.85 (5.0π/M) |
| Hamming | −42.7 dB | −54.0 dB | 0.0019 | 6.66π/M | 5.0 (6.4π/M) |
| Blackman | −58.1 dB | −75.3 dB | 0.0002 | 11.2π/M | 7.3 (9.4π/M) |
| Kaiser, A = 60 | −41.4 dB | −59.8 dB | 0.0010 | 7.23π/M | β = 5.65 |
Measured on M = 200 low-pass filters with \(\omega_c = 0.4\pi\); transition width measured between the points where the error equals the ripple, written as \(\Delta\omega\). The Kaiser formula predicts \(\Delta\omega\cdot M = (A - 8)/2.285\) within 1 % over 30–90 dB. Kaiser is never worse than the fixed windows for the same attenuation, because β tunes the ripple and the length independently.
Worked example (RF anti-alias low-pass). \(f_s\) = 40 MHz, pass ≤ 8 MHz, stop ≥ 9 MHz, A = 60 dB. \(\Delta\omega = 2\pi\cdot 1/40 = 0.1571\) rad/sample, β = 5.653, \(M = 52/(2.285 \times 0.1571) = 144.9 \to 145\) (09.3). Measured: 59.6 dB, short by 0.4 dB (M = 145 is odd, so the delay is a half sample; the even order 146 gives an integer delay and measures 59.9 dB). The formulas are good to about ±0.5 dB, and raising M alone does not help because β sets the floor (−59.8 dB at any length). Ask for A = 61: β = 5.76, M = 148, measured 60.9 dB (passband deviation 0.0009). Delay M/2 = 74 samples = 1.85 µs, which is 1.42 mm of depth at c/2 = 770 m/s.
09.3Optimal FIRs: frequency sampling, least squares, equiripple
Write a type I linear-phase filter (M even, symmetric taps) as a delay times a real amplitude: \(H(e^{j\omega}) = e^{-j\omega M/2} A(\omega)\), with \(A(\omega) = a_0 + \sum_{k=1}^{M/2} a_k \cos k\omega\) and \(h[M/2 \pm k] = a_k/2\). Every design method chooses the M/2 + 1 numbers \(a_k\) by a different criterion on the error \(E(\omega) = W(\omega)\left(A(\omega) - D(\omega)\right)\) over the bands, ignoring the transition ("don't care") band.
| Method | Criterion | Solve by | Character |
|---|---|---|---|
| Window (Kaiser) | None: truncate and taper the ideal \(h_d\) | Closed form | \(\delta_p = \delta_s\); ripple largest near the edge; simplest, predictable |
| Frequency sampling | \(A(\omega_k) = D(\omega_k)\) at \(\omega_k = 2\pi k/(M + 1)\) | Inverse DFT | Exact at the samples, ripples between; 1–2 optimized transition samples lift the stopband from −16 dB to about −42 dB (one sample at 0.39) |
| Least squares | \(\min \int W(\omega)\lvert A - D\rvert^2\, d\omega\) | Normal equations (M/2 + 1 unknowns) | Lowest error energy; peak error stays at the band edges. With no transition band and W = 1 it reduces to the rectangular window |
| Equiripple (Parks–McClellan) | \(\min \max \lvert E(\omega)\rvert\) | Remez exchange on the alternation theorem | Optimal for the peak error: ≥ M/2 + 2 equal-height extrema; weight \(W_s/W_p = \delta_p/\delta_s\) sets the ripple ratio |
- \(M\)
- filter order (\(M + 1\) taps)
- \(\delta_p,\ \delta_s\)
- passband and stopband ripple (linear amplitude)
- \(\Delta f\)
- transition width normalized to \(f_s\), cycles/sample
- \(A_s\)
- stopband attenuation, dB
- \(\log_{10}\)
- base-10 logarithm
Worked comparison for the anti-alias spec of 09.2: pass ≤ \(0.2\,f_s\), stop ≥ \(0.225\,f_s\) (8/9 MHz at 40 MHz).
| Design | Estimate | True minimum order | Multiplies/output (folded) |
|---|---|---|---|
| Equiripple, \(\delta_p = 0.01\) (±0.09 dB), \(\delta_s = 0.001\) (60 dB) | 101.4 (Kaiser), 101.4 (Herrmann) | 104 | 53 |
| Equiripple, \(\delta_p = \delta_s = 0.001\) | 128.8 | 132 | 67 |
| Kaiser window, \(\delta = 0.001\) | 144.9 | 148 (with A = 61) | 75 |
| A/22 rule | 109 taps | — | — |
True minima found by scanning the order with a Remez design (SciPy, grid density 64) and checking both ripples on a 64k-point grid; order 104 meets the first spec with 1 % to spare. The demo and m09_fir_window_kaiser.m approach the equiripple solution with Lawson's iteratively reweighted least squares, which needs only a linear solver.
FIR design bench: window vs least squares vs equiripple
SimulationDesigns computed in the page: Kaiser window (09.3), weighted least squares by the normal equations, equiripple approached by Lawson's iteratively reweighted least squares; responses by direct DTFT of the taps; JS.09.4IIR design: analog prototypes, the bilinear transform and biquads
IIR design borrows a century of analog filter theory. Design a normalized analog low-pass prototype \(H_a(s)\), transform it (low-pass to low-pass, high-pass or band-pass by substituting s), then map s to z. The bilinear transform is the standard map:
- \(s,\ \Omega\)
- Laplace variable and analog frequency, rad/s
- \(z\)
- z-transform variable
- \(T\)
- sampling interval \(1/f_s\), s
- \(\omega\)
- digital frequency, rad/sample
- \(\Omega_p,\ \omega_p\)
- analog (prewarped) and digital passband edge
- \(j\)
- imaginary unit
- \(A_p,\ A_s\)
- passband ripple (peak to peak) and stopband attenuation, dB
- \(\varepsilon_p,\ \varepsilon_s\)
- ripple parameters (dimensionless)
- \(\Omega_p,\ \Omega_s\)
- prewarped passband and stopband edges, rad/s
- \(N\)
- analog prototype order
- \(K,\ K'\)
- complete elliptic integral of the first kind and its complement, \(K'(k) = K(\sqrt{1 - k^2})\)
- \(k,\ k_1\)
- selectivity and discrimination moduli
- \(\ln,\ \operatorname{acosh}\)
- natural logarithm, inverse hyperbolic cosine
The same spec as 09.3 in IIR terms: 8/9 MHz at 40 MHz, \(A_p\) = 0.17 dB (the ±0.01 passband), \(A_s\) = 60 dB.
| Design | Order | Biquads | Multiplies/output | Phase |
|---|---|---|---|---|
| Butterworth | 53 | 27 | ≈ 135 | Nonlinear, smoothest |
| Chebyshev I or II | 16 | 8 | ≈ 40 | Nonlinear |
| Elliptic | 8 | 4 | ≈ 20 | Most nonlinear near the edge |
| Equiripple FIR (\(\delta_p = 0.01\)) | 104 | — | 53 | Exactly linear, delay 52 samples |
Orders from the formulas above with prewarped edges (they agree with SciPy's buttord, cheb1ord, cheb2ord, ellipord). Five multiplies per biquad. The elliptic filter is cheaper than the FIR by less than a factor of 3 here, and the FIR's cost falls by the decimation factor in a polyphase decimator (Section 09.6), which is why front ends are mostly FIR.
Ship it as biquads. Factor H into second-order sections, pair each complex pole pair with the nearest zero pair, and order the sections by pole radius. The reason is root sensitivity: for a polynomial with roots \(p_i\), \(\partial p_i/\partial a_k = -p_i^{N-k} / \prod_{j \ne i}(p_i - p_j)\). Narrowband filters cluster their poles near z = 1, the denominators collapse, and tiny coefficient errors move poles across the circle. Each biquad has \(a_1 \in (-2, 2)\) and \(a_2 \in (-1, 1)\) and only two roots to protect (Figure 09.B).
Section transfer function \(H_k(z) = (b_0 + b_1 z^{-1} + b_2 z^{-2})/(1 + a_1 z^{-1} + a_2 z^{-2})\). The two states hold partial sums, so only two delays are needed per section; quantizing \(a_1, a_2\) moves only this section's pole pair.
| Butterworth low-pass | Coefficients | Max pole radius, direct form | Max pole radius, biquads |
|---|---|---|---|
| N = 8, \(f_c = 0.02\,f_s\) (true 0.9758) | 16-bit words | 1.299 (unstable) | 0.9758 |
| N = 8, \(f_c = 0.02\,f_s\) | 24-bit words | 1.166 (unstable) | 0.9758 |
| N = 12, \(f_c = 0.01\,f_s\) | double precision | output diverges (NaN) | stable, matches design |
Direct form: the \(a_k\) of the 8th-order filter reach 50 in magnitude, so a 16-bit word keeps only 9 fractional bits; biquads use 2 integer bits. Computed with m09_iir_bilinear_butter_cheby.m.
IIR families at the same order
Exact formulaClosed-form analog prototypes (elliptic by Landen iteration), prewarped bilinear transform (09.5), exact magnitude and group delay of the resulting zpk designs; JS.09.5Phase, group delay, pulse distortion and zero-phase filtering
Recap from Module 03. Group delay \(\tau_g = -d\varphi/d\omega\) and phase delay \(\tau_\varphi = -\varphi/\omega\) are defined in Section 03.5, (03.16): a narrowband pulse \(a(t)\cos(\omega_0 t)\) through H comes out as approximately \(\lvert H(\omega_0)\rvert\, a(t - \tau_g)\cos\left(\omega_0 (t - \tau_\varphi)\right)\), and linear phase \(\varphi = -\omega M/2\) gives \(\tau_g = \tau_\varphi = M/2\) at every \(\omega\), so the shape is preserved. The design question here is what happens when a filter is not linear phase.
The approximation holds when \(\tau_g\) is nearly constant across the pulse band. When it varies (an IIR near its band edge), different frequencies arrive at different times: the pulse spreads and its envelope becomes asymmetric, with a tail on the side of the steeper edge. A mismatch \(\tau_\varphi \ne \tau_g\) rotates the carrier under the envelope, which a phase-based delay estimator (Module 11) reads as displacement.
FIR vs IIR meeting the same RF band-pass spec
SimulationEquiripple FIR (IRLS) and bilinear IIR band-passes designed in the page; a Gaussian-envelope RF pulse filtered by convolution and by biquad recursion, correlation after best integer alignment; JS.Zero-phase filtering. Filter forward, reverse the result, filter again, reverse back. The overall response is \(H(e^{j\omega})H^*(e^{j\omega}) = \lvert H\rvert^2\): no phase, no delay, magnitude squared.
- \(x,\ y\)
- input and output records
- \(h\)
- filter impulse response
- \(*\)
- linear (causal) convolution
- \(\mathcal{R}\{\cdot\}\)
- time reversal of a record
- \(H(e^{j\omega}),\ H_{\mathrm{eff}}\)
- single-pass and forward-backward frequency responses
- \(\omega\)
- digital frequency, rad/sample
- \(N\)
- Butterworth prototype order
- \(\Omega_c\)
- original −3 dB frequency, rad/s
filtfilt do.
Causal vs forward–backward, and the edges
SimulationSynthetic RF record (three 5 MHz echoes, offset and drift) passed through a bilinear Butterworth band-pass as biquads, one pass and forward–backward with and without padding; the reference filters a record extended far beyond both ends; JS.Real-time RF and IQ chains use linear-phase FIRs and subtract the known delay M/2, identically on every channel; a channel-dependent filter delay is an aberration you introduced yourself. Zero-phase IIR filtering is fine offline: stored RF for QUS spectra, slow-time signals in shear-wave elastography post-processing, displacement maps. Never mix a zero-phase reference with a causal measurement in a delay or phase comparison.
09.6Multirate and special-purpose filters: decimation, half-band, polyphase, CIC, Hilbert
Decimation by D keeps every D-th output: \(y[m] = \sum_k h[k]\, x[mD - k]\). Everything between \(f_s/(2D)\) and \(f_s/2\) folds into the new band, so the low-pass must remove it first. If only \(0 \dots f_{\text{pass}}\) must stay clean, aliases may land in the transition band: \(f_{\text{stop}} = f_s/D - f_{\text{pass}}\). Polyphase splits h into D sub-filters \(h_i[m] = h[mD + i]\) that run at the output rate, so only kept samples are computed: \((M + 1)/D\) multiplies per input sample instead of M + 1.
- \(\omega_c\)
- cutoff, rad/sample
- \(\delta_p,\ \delta_s\)
- passband and stopband ripple
- \(h[n]\)
- impulse response
- \(M\)
- filter order, even
- \(k\)
- positive integer index of the zeroed even-offset taps
- \(H(z),\ H(f)\)
- CIC transfer function and its magnitude response
- \(z\)
- z-transform variable
- \(R\)
- decimation factor
- \(M\)
- differential delay of each comb (here \(M\) is not the FIR order)
- \(N\)
- number of integrator-comb stages
- \(f\)
- frequency normalized to the input rate, cycles/sample
- \(f_{\text{out}}\)
- output rate after decimation by \(R\)
- \(h[n]\)
- Hilbert FIR impulse response (type III, odd length \(M + 1\))
- \(M\)
- filter order, even
- \(n\)
- offset from the center tap; odd values carry the \(2/\pi n\) taps
- \(H\)
- frequency response of the 31-tap filter
Coefficient quantization moves poles and zeros (Section 09.4); signal rounding at each multiplier adds noise of variance \(q^2/12\) shaped by the transfer function from that node to the output, so a high-Q IIR section amplifies its own rounding noise by the gain of its pole pair. Scale each section so its internal nodes cannot overflow, put the highest-Q section last, and beware zero-input limit cycles in recursive sections. FIRs have none of the last two problems.
09.7Doppler wall filters on short ensembles
Color flow estimates velocity from N = 8–16 slow-time samples per pixel at the PRF. Stationary and slowly moving tissue can be 40 dB or more above blood and sits near f = 0; the wall (clutter) filter must remove it before the lag-one autocorrelation (Module 08). With so few samples, "filter" means a matrix: every linear wall filter is \(\mathbf{y} = \mathbf{A}\mathbf{x}\) with \(\mathbf{A}\) of size \(N_{\text{out}} \times N\), and the frequency response of a time-varying \(\mathbf{A}\) is defined through the average output power for a unit complex exponential input:
- \(\lvert H(f)\rvert^2\)
- average power gain of the filter at Doppler frequency \(f\)
- \(\mathbf{A}\)
- wall-filter matrix, \(N_{\text{out}} \times N\)
- \(N\)
- ensemble length (slow-time samples per pixel)
- \(N_{\text{out}}\)
- number of valid outputs
- \(\mathbf{e}_f,\ e_f[n]\)
- unit-amplitude complex exponential at Doppler frequency \(f\) and its \(n\)th sample
- \(\lVert\cdot\rVert\)
- Euclidean norm
- \(n\)
- slow-time sample index
- \(f,\ \mathrm{PRF}\)
- Doppler frequency and pulse repetition frequency, Hz
- \(j\)
- imaginary unit
| Filter | Matrix A | Outputs | Behaviour on short ensembles |
|---|---|---|---|
| FIR, order K | Toeplitz band of h, valid outputs only | N − K | Exact steady state, but every order costs one output; sharp cutoffs need K ≈ N |
| IIR, zero initial state | Lower-triangular Toeplitz of the impulse response | N | The first samples are transient: on N = 10, a 2nd-order Butterworth HP at 0.1 PRF suppresses DC by only 12.6 dB |
| IIR, step initialization | State set to steady state for constant x[0] | N | Exact null at DC; at 0.01 PRF still only −32 dB |
| IIR, projection initialization | \(\mathbf{A} = (\mathbf{I} - \mathbf{G}\mathbf{G}^{+})\mathbf{F}\); \(\mathbf{G}\) = zero-input responses | N | Initial state chosen to minimize output energy; −43 dB at 0.01 PRF |
| Polynomial regression, degree p | \(\mathbf{A} = \mathbf{I} - \mathbf{Q}\mathbf{Q}^H\), \(\mathbf{Q}\) = orthonormal polynomials of degree ≤ p | N | Removes p + 1 dimensions exactly; −63 dB at 0.01 PRF for p = 2; −3 dB at \(\approx (0.44 + 0.41p)\,\mathrm{PRF}/N\) |
| SVD / eigen-filter (Module 06) | \(\mathbf{A} = \mathbf{I} - \mathbf{U}_c\mathbf{U}_c^H\) from the data | N | Regression with data-derived basis; separates by spatial coherence when tissue moves |
Projection initialization written out: the IIR output from initial state \(\mathbf{s}_0\) is \(\mathbf{y} = \mathbf{F}\mathbf{x} + \mathbf{G}\mathbf{s}_0\), where \(\mathbf{F}\) is the zero-state convolution matrix and the columns of \(\mathbf{G}\) are the zero-input responses to unit initial states. Minimizing \(\lVert\mathbf{y}\rVert^2\) over \(\mathbf{s}_0\) gives \(\mathbf{s}_0 = -\mathbf{G}^{+}\mathbf{F}\mathbf{x}\), hence \(\mathbf{A} = (\mathbf{I} - \mathbf{G}\mathbf{G}^{+})\mathbf{F}\): the IIR filter followed by a projection that removes everything its own transient could explain. Regression filters are the same idea with a polynomial basis, and SVD filters with an empirical one. Numbers in the table: N = 10, cutoff 0.1 PRF, computed with m09_wall_filters_ensemble.m.
Wall filters on an N-sample ensemble: response and transient
SimulationWall filters built as matrices (FIR, IIR with zero, step or projection initialization, Legendre regression); power responses from (09.11) are exact for each matrix; the lower panel filters a synthetic ensemble (clutter tone + blood tone); JS.Convert the stopband to velocity with \(v = c f/(2 f_0\cos\theta)\). At \(f_0\) = 5 MHz, PRF = 4 kHz and N = 10, a degree-2 regression filter is −3 dB at 1.27 PRF/N = 507 Hz, i.e. 7.8 cm/s at θ = 0: slower flow is gone. The levers are a longer ensemble (frame rate), a lower PRF (aliasing), or a filter that separates tissue by something other than Doppler frequency, the spatial coherence exploited by the SVD filters of Module 06 on the hundreds of frames of ultrafast Doppler and ULM.
09.8Visual summary
Thirteen pictures that carry the module. Each card is a question to answer in five seconds, the answer with its number, and a figure computed from the formula or from the actual design.
09.9How experts think about it
A filter expert sees three things at once: the pole-zero map (where the response comes from), the error budget (\(\delta_p\), \(\delta_s\) and the transition width, which set the order before any design runs), and the implementation (delay, transient, word length). FIRs are chosen for linear phase, polyphase cost and robustness; IIRs for low order when phase does not matter or can be cancelled offline. On short records, such as Doppler ensembles, the steady-state frequency response is a fiction and the filter must be judged as a matrix.
Each factor \((1 - z_k e^{-j\omega})\) has the magnitude of the vector from \(z_k\) to \(e^{j\omega}\) ((03.23), Module 03). Walk the unit circle and watch the distances. For the resonator with poles at \(0.95\,e^{\pm j\pi/4}\) (no zeros), at \(\omega = \pi/4\) the near pole is 0.05 away and the far pole 1.379 away, so \(\lvert H\rvert = 1/(0.05 \times 1.379) = 14.5\); at DC both are \(\sqrt{0.559}\) away, \(\lvert H\rvert = 1.79\). The peak stands 18.2 dB above DC, and you read it off the picture without computing a transform.
Near the pole, \(\lvert e^{j\omega} - p\rvert^2 \approx (1 - r)^2 + (\omega - \theta)^2\), so the power falls by half at \(\lvert\omega - \theta\rvert = 1 - r\): full −3 dB width \(2(1 - r)\) rad/sample, i.e. \((1 - r)f_s/\pi\) in Hz. The impulse response decays as \(r^n\), time constant \(-1/\ln r \approx 1/(1 - r)\). The product of width and duration is ≈ 2, the time–bandwidth trade-off in pole form. Figure: normalized one-pole peaks for r = 0.8, 0.9, 0.95 against the −3 dB line.
A pole contributes \(\tau = \frac{r\cos(\omega - \theta) - r^2}{1 - 2r\cos(\omega - \theta) + r^2}\) (09.1), which at \(\omega = \theta\) is \(r(1 - r)/(1 - r)^2 = r/(1 - r)\). Sharp transitions need poles close to the circle, so the delay piles up exactly at the band edges, while the passband center can have little delay. The order-10 elliptic band-pass of Demo 09.4 delays 5 MHz by only 8.8 samples (0.22 µs) but its 3 MHz edge by 86 samples; the equiripple FIR delays everything by 36. Figure: delay of a pole pair at θ = π/4.
Symmetry gives \(H(z) = z^{-M}H(z^{-1})\), so if \(z_0\) is a zero, so is \(1/z_0\) (and the conjugates, for real taps): zeros sit on the circle or in quadruples \(z_0,\ 1/z_0,\ z_0^*,\ 1/z_0^*\). Evaluating at z = −1 for odd M gives \(H(-1) = -H(-1) = 0\): an even-length symmetric FIR cannot be high-pass. Antisymmetric taps force \(H(1) = 0\), plus \(H(-1) = 0\) for odd length: Hilbert transformers and differentiators. Figure: zeros of a 21-tap Kaiser low-pass: stopband zeros on the circle, passband zeros in reciprocal pairs.
Kaiser fitted these formulas (09.3) to designs with his \(I_0\)-sinh window, which approximates the prolate spheroidal sequence (the most energy in a given main lobe). The example: \(\Delta\omega = 2\pi/40 = 0.157\) rad/sample, \(M = 52/0.359 = 144.9\). The fit is good to about ±0.5 dB; the floor is set by β (−59.8 dB here at any length), so the fix is A = 61: β = 5.76, M = 148, measured 60.9 dB. Delay 74 samples, 1.85 µs, 1.42 mm. Figure: that design against the mask.
The designed response is the ideal edge convolved with \(W(\omega)\) (09.2). Lengthening the window compresses W along ω but keeps the ratio of sidelobe area to main-lobe area, so the overshoot and the stopband floor stay put while the transition narrows as 1/M. That is the Gibbs phenomenon: a truncated sinc overshoots by 8.9 % (−21 dB stopband) whether it has 21 or 81 taps. By symmetry of the convolution, every window design passes through 0.5 at \(\omega_c\). Figure: rectangular-window low-pass at M = 20 and M = 80, linear amplitude near the edge.
A window design spends its error unevenly (largest at the edges, decaying away) and ties \(\delta_p\) to \(\delta_s\). Parks–McClellan spreads the error into M/2 + 2 equal extrema (the alternation theorem) and weights the bands, so a loose passband buys stopband. Order estimates: Kaiser's \((-20\log_{10}\sqrt{\delta_p\delta_s} - 13)/(14.6\,\Delta f)\) (09.4) gives 101.4 for the 8/9 MHz spec; the true minimum is 104. Figure: same spec with a 2 MHz transition: equiripple order 54 (blue; true minimum 53) vs Kaiser order 74 (pink).
A transition of width \(\Delta f\) needs an impulse response lasting about \(1/\Delta f\) samples (time–bandwidth), and every extra 20 dB of stopband costs a fixed multiple of that. All three estimates (Kaiser window, Kaiser/Herrmann equiripple, A/22) are of the form \(A/\Delta f\). Consequence for system design: the same analog transition costs more taps at a higher sampling rate, 109 taps at 40 MHz but 170 at 62.5 MHz, unless you decimate first. Figure: A/22 rule for 40, 60, 80 dB versus Δf.
Decimation keeps \(y[mD]\). Writing h as D phases \(h_i[m] = h[mD + i]\) and the input as D interleaved streams, each kept output is a sum of D short convolutions running at the output rate, so the work per input sample is \((M + 1)/D\). A half-band filter adds two more factors of 2: every other tap is zero and symmetry folds the rest. Figure: multiplies per input sample for a 64-tap filter in polyphase form (bars) against direct form (dashed line at 64), versus D.
Without prewarping, an analog frequency \(f_a\) lands at \(f_d = (f_s/\pi)\operatorname{atan}(\pi f_a/f_s)\). At \(f_a\) = 4 MHz (\(f_s/10\)) the error is −3.1 %; at 15 MHz the edge lands at 11.04 MHz. Prewarping designs the analog prototype at \(\Omega_p = (2/T)\tan(\omega_p/2)\) (09.5), so the chosen edges land exactly; the shape between edges is warped, which is harmless for piecewise-constant specs. Figure: \(f_d\) versus \(f_a\) at \(f_s\) = 40 MHz, and the identity line.
With \(\varepsilon_p, \varepsilon_s\) from \(A_p, A_s\) (09.6), the stopband edge relative to the passband edge is \((\varepsilon_s/\varepsilon_p)^{1/N}\) for Butterworth, \(\cosh\left(\operatorname{acosh}(\varepsilon_s/\varepsilon_p)/N\right)\) for both Chebyshev types and \(1/k\) for elliptic. For N = 6, \(A_p\) = 1 dB, \(A_s\) = 60 dB: 3.54, 2.11 and 1.34. Allowing ripple in a band is what buys selectivity, and the poles that make it move closer to the jΩ axis, which is the group-delay peak. Figure: the three analog prototypes, N = 6.
Root sensitivity \(\partial p_i/\partial a_k = -p_i^{N-k}/\prod_{j \ne i}(p_i - p_j)\) blows up when poles cluster, and narrowband low-pass or wall filters cluster them near z = 1. The coefficients of the 8th-order filter at \(0.02\,f_s\) reach 50, so a 16-bit word keeps 9 fractional bits and the poles scatter outside the circle. A biquad has only two roots, coefficients below 2 in magnitude, and gets 14 fractional bits. Figure: true poles (blue) and 16-bit direct-form poles (pink) near z = 1.
Filtering a time-reversed signal applies \(H^*(e^{j\omega})\), so the cascade is \(\lvert H\rvert^2\) (09.7): zero phase, every dB figure doubled. For Butterworth, \(\lvert H\rvert^2 = 1/(1 + (\Omega/\Omega_c)^{2N})\) reaches \(1/\sqrt{2}\) at \(\Omega = (\sqrt{2} - 1)^{1/(2N)}\Omega_c\), \(0.896\,\Omega_c\) for N = 4. Design for the squared spec (half the dB in each pass). The two passes start transients at both record edges; reflection padding plus steady-state initial states removes most of it. Figure: \(\lvert H\rvert\) (pink) and \(\lvert H\rvert^2\) (blue) for a 4th-order Butterworth at \(0.1\,f_s\).
The steady-state null at DC is reached only after the poles' ringing time, about \(1/(1 - r)\) pulses, which is comparable to the whole ensemble. Treat the filter as the matrix it is (09.11): with zero state the clutter step at n = 0 rings through every output. Step initialization nulls exact DC; projection initialization, \(\mathbf{A} = (\mathbf{I} - \mathbf{G}\mathbf{G}^{+})\mathbf{F}\), also removes what the transient could explain and keeps all N outputs. Figure: power responses for N = 10, \(f_c\) = 0.1 PRF: zero state (gray), step (gold), projection (pink).
Over N samples, a slow complex exponential looks like a low-order polynomial, so projecting out polynomials of degree ≤ p removes it; the response is \(1 - \lVert\mathbf{Q}^H\mathbf{e}_f\rVert^2/N\). A tone escapes once it completes enough cycles over the ensemble, so the cutoff scales as 1/N, and each extra degree moves it up by about 0.4 PRF/N (computed: 0.44, 0.86, 1.26, 1.66 PRF/N for p = 0…3 at N = 12). Figure: N = 12, p = 0…3.
N integrators at the input rate and N combs \((1 - z^{-M})\) at the output rate implement \(\left[(1 - z^{-RM})/(1 - z^{-1})\right]^N\) (09.9), a boxcar of length RM raised to the N-th power: nulls at multiples of \(f_{\text{in}}/(RM)\), exactly where the aliases of DC land. At \(f_{\text{out}}/4\), \(\mathrm{sinc}(1/4) = 0.900\) per stage, −0.91 dB, so 3 stages droop 2.7 dB and need a compensating FIR. Registers need \(N\log_2(RM)\) extra bits; two's-complement wraparound in the integrators is harmless if the width is right. Figure: one to four stages, R = 8, M = 1.
The ideal response \(-j\,\mathrm{sgn}(\omega)\) has impulse response \(2/(\pi n)\) for odd n and 0 for even n (09.10). A type III FIR (odd length, antisymmetric) is forced to zero at \(z = \pm 1\), so the passband must stay away from both ends; RF at 2–8 MHz sampled at 40 MHz sits comfortably inside. 31 taps with 16 nonzero antisymmetric values cost 8 multiplies per output, and the real path must be delayed by the same 15 samples to form \(x + j\mathcal{H}\{x\}\).
09.10Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Response from geometry | \(\lvert H\rvert = \lvert g\rvert \prod\lvert e^{j\omega} - z_k\rvert / \prod\lvert e^{j\omega} - p_k\rvert\) | Placing notches, wall-filter zeros at DC |
| Pole near the circle | width \(2(1 - r)\) rad/sample, ringing \(1/(1 - r)\), delay \(r/(1 - r)\) | Ringing in pulses vs ensemble length |
| Linear phase | \(h[n] = \pm h[M - n] \Rightarrow \tau_g = M/2\) | Channel-consistent RF/IQ filtering |
| Kaiser β | \(0.1102(A - 8.7)\), A > 50; \(0.5842(A - 21)^{0.4} + 0.07886(A - 21)\), 21–50 | Quick FIR design without a toolbox |
| Kaiser order | \(M = (A - 8)/(2.285\,\Delta\omega)\) | Sizing RF band-pass and anti-alias filters |
| Equiripple order | \(M \approx (-20\log_{10}\sqrt{\delta_p\delta_s} - 13)/(14.6\,\Delta f)\); taps \(\approx A/(22\,\Delta f)\) | Front-end FIR budgets (FPGA multipliers) |
| Bilinear transform | \(s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}\), \(\Omega = \frac{2}{T}\tan(\omega/2)\) | Digitizing analog prototypes, prewarped edges |
| IIR orders | Butterworth \(\ln(\varepsilon_s/\varepsilon_p)/\ln(\Omega_s/\Omega_p)\); Chebyshev \(\operatorname{acosh}(\varepsilon_s/\varepsilon_p)/\operatorname{acosh}(\Omega_s/\Omega_p)\) | Choosing the family for a wall or smoothing filter |
| Biquads | \(H = \prod \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}\) | Any IIR above order 2, fixed point or not |
| Group vs phase delay | \(\tau_g = -d\varphi/d\omega\), \(\tau_\varphi = -\varphi/\omega\) | Envelope vs carrier shift; phase-based displacement bias |
| Forward–backward | \(H_{\mathrm{eff}} = \lvert H\rvert^2\), zero phase | Offline RF, QUS, SWE displacement smoothing |
| Decimation | \(f_{\text{stop}} = f_s/D - f_{\text{pass}}\); polyphase \((M + 1)/D\) per input | RF → IQ chains, data-rate reduction |
| Half-band | \(h[M/2] = \tfrac{1}{2}\), \(h[M/2 \pm 2k] = 0\) | Cascaded ↓2 stages |
| CIC | \(\left[(1 - z^{-RM})/(1 - z^{-1})\right]^N\), bit growth \(N\log_2(RM)\) | First decimation stage after oversampled ADCs |
| Hilbert FIR | \(h[n] = 2/(\pi n)\), n odd (centered), type III | Analytic signal and envelope on RF |
| Wall-filter response | \(\lvert H(f)\rvert^2 = \lVert\mathbf{A}\mathbf{e}_f\rVert^2/N_{\text{out}}\) | Comparing FIR, IIR, regression, SVD filters |
| Projection init | \(\mathbf{A} = (\mathbf{I} - \mathbf{G}\mathbf{G}^{+})\mathbf{F}\) | IIR wall filters on 8–16 pulses |
| Regression cutoff | \(f_{-3\,\mathrm{dB}} \approx (0.44 + 0.41p)\,\mathrm{PRF}/N\) | Lowest detectable velocity \(v = c f/(2 f_0\cos\theta)\) |
09.11Interview and qualifying-exam questions
I1Why are RF and IQ filters in a beamformer almost always FIR? · screen
Linear phase: a constant, known delay \(M/2\) that is identical on every channel, so filtering never looks like an aberration and the delay is subtracted once. Polyphase decimation divides the cost by D, which erases most of an IIR's per-output advantage (for the 8/9 MHz spec, 53 folded multiplies vs about 20 for an 8th-order elliptic, before decimation). No stability, limit-cycle or coefficient-sensitivity problems in fixed point. Follow-up: how many taps for 3–7 MHz at 40 MHz with 1 MHz transitions and 40 dB? The low-pass estimate says about 61; two transitions push the true equiripple minimum to 72.
I2A 129-tap linear-phase FIR runs at 50 MHz on RF. What delay do you remove, in samples, µs and mm? · screen
Order 128, delay 64 samples = 1.28 µs. In pulse-echo depth that is \(c\,t/2\) = 770 m/s × 1.28 µs = 0.99 mm. Follow-up: what if the filter has an even number of taps? Then the delay is a half sample; either pad to odd length or account for the half-sample shift in the delay law, otherwise IQ phase is biased by \(\pi f_0/f_s\).
I3Derive the bandwidth and the peak group delay of a pole at \(r e^{j\theta}\). · onsite
\(\lvert 1 - r e^{j\theta}e^{-j\omega}\rvert^2\) \(= 1 - 2r\cos(\omega - \theta) + r^2\) \(\approx (1 - r)^2 + r(\omega - \theta)^2\) near θ, so half power at \(\lvert\omega - \theta\rvert \approx 1 - r\): full width \(2(1 - r)\) rad/sample. The phase of \(1/(1 - re^{-j\Delta})\) is \(-\operatorname{atan}\left(r\sin\Delta/(1 - r\cos\Delta)\right)\); differentiating at Δ = 0 gives the group delay \(r/(1 - r)\). r = 0.95: 0.1 rad/sample and 19 samples. Follow-up: convert to a color-flow PRF of 4 kHz: 64 Hz wide, 4.75 ms delay, longer than a 12-pulse ensemble (3 ms).
I4Size a Kaiser-window low-pass: 50 dB, 0.5 MHz transition at \(f_s\) = 20 MHz. Then estimate equiripple with \(\delta_p = 0.01\). · onsite
\(\Delta f = 0.025\), \(\Delta\omega = 0.157\) rad/sample. \(\beta = 0.5842 \times 29^{0.4} + 0.07886 \times 29 = 4.53\). \(M = 42/(2.285 \times 0.157) = 117\), round to 118 for an integer delay. Equiripple: \(\delta_s = 0.00316\), \(-20\log_{10}\sqrt{0.01 \times 0.00316} = 45.0\) dB, \(M \approx (45.0 - 13)/(14.6 \times 0.025) = 88\), so budget about 90 after the few-percent underestimate. Follow-up: why can the window design not take advantage of the 0.01 passband? Because its \(\delta_p\) equals \(\delta_s\) by construction.
I5Why does the bilinear transform need prewarping, and why not use impulse invariance instead? · onsite
The bilinear map sends Ω to \(\omega = 2\operatorname{atan}(\Omega T/2)\), so analog edges land at lower digital frequencies (15 MHz lands at 11.0 MHz for \(f_s\) = 40 MHz). Prewarping designs the prototype at \((2/T)\tan(\omega_p/2)\) so the edges land where specified. Impulse invariance samples \(h_a(t)\), which aliases its spectrum: fine for steep low-pass prototypes, impossible for high-pass and band-stop, whose response does not fall off. Follow-up: what does the bilinear transform do to an analog zero at infinity? It maps to z = −1, a null at Nyquist.
I6Your 8th-order elliptic low-pass works in MATLAB but oscillates on the FPGA. Diagnose. · onsite
Almost certainly direct-form coefficients: the \(a_k\) of a narrowband 8th-order filter span tens, so after quantization the clustered poles move outside the unit circle (16-bit direct form moves an 8th-order Butterworth pole at \(0.02\,f_s\) to radius 1.30). Re-implement as four biquads with poles paired to nearest zeros, scale each section against overflow, put the highest-Q section last. If it still hums with zero input, it is a limit cycle: add more state bits or use magnitude truncation / error feedback. Follow-up: why are FIRs immune? No feedback: quantization moves zeros, which changes the stopband floor but cannot destabilize.
I7Color flow shows strong clutter bleeding into the velocity map near vessel walls. The wall filter is a 2nd-order Butterworth high-pass on 10 pulses. What do you check first? · research
The initialization. With zero initial state, the clutter step at the first pulse rings through the ensemble: on N = 10 at \(f_c\) = 0.1 PRF the DC suppression is only 12.6 dB instead of infinite. Step initialization nulls DC exactly but gives only −32 dB at 0.01 PRF; projection initialization gives −43 dB with all 10 outputs kept; a degree-2 regression filter −63 dB. Follow-up: what does projection initialization compute? \(\mathbf{s}_0 = -\mathbf{G}^{+}\mathbf{F}\mathbf{x}\), minimizing output energy over the initial state, so \(\mathbf{A} = (\mathbf{I} - \mathbf{G}\mathbf{G}^{+})\mathbf{F}\).
I8When is a polynomial regression wall filter enough, and when do you need an SVD filter? · research
Regression assumes clutter is a low-order polynomial in slow time at every pixel independently: fine for slowly moving tissue and short ensembles where its cutoff (\(\approx (0.44 + 0.41p)\,\mathrm{PRF}/N\)) sits below the slowest flow of interest. When tissue moves at velocities overlapping blood (cardiac, breathing, handheld probe) or when slow microvascular flow must be kept, frequency separation fails; an SVD filter on a spatiotemporal Casorati matrix separates by spatial coherence (Module 06, Demené et al. 2015) and needs long ensembles (ultrafast Doppler, ULM). Follow-up: what is the regression filter's matrix? \(\mathbf{I} - \mathbf{Q}\mathbf{Q}^H\) with \(\mathbf{Q}\) orthonormal polynomials; an SVD filter is the same with \(\mathbf{Q}\) from the data.
I9You smooth shear-wave displacement vs time with filtfilt and a 4th-order Butterworth at 1 kHz. What did you actually apply? · onsite
\(\lvert H\rvert^2\): zero phase (so arrival times are not biased, which matters for time-of-flight SWE), but −6 dB at 1 kHz, the −3 dB point at 0.896 kHz, and an 8th-order magnitude roll-off. The first and last few periods of the cutoff are affected by edge transients unless the record is reflection-padded and the filter started from steady state. Follow-up: could you do this in real time? Only block-wise with overlap and a latency of at least the block length.
I10You need envelopes from RF on an FPGA at 40 MHz. FFT Hilbert or FIR Hilbert, and how big? · research
FIR: streaming, no block edges, no circular wrap. A 31-tap Hamming-windowed type III Hilbert is within −0.19/+0.05 dB from 2 to 18 MHz; only 16 taps are nonzero and they are antisymmetric, so 8 multiplies per output. Delay the in-phase path by 15 samples to match. Better still, demodulate to IQ and decimate (Module 08), which also cuts the data rate. Follow-up: why not an even-length Hilbert? Type IV has a half-sample delay, awkward to match in the real path.
09.12MATLAB lab
Base MATLAB only (no toolboxes): the Kaiser window, Butterworth/Chebyshev prototypes, the bilinear transform, biquads, IRLS equiripple and the wall-filter matrices are implemented in each script with filter, conv, fft, roots, poly and backslash. Each runs as-is in MATLAB R2019b+ and GNU Octave 8. Runtimes are Octave timings.
m09_pole_zero_response.mResponse from pole-zero geometry, pole width and delay rules, linear-phase zeros · change r, θ · runtime ≈ 1 s
%% m09_pole_zero_response.m
% Module 09 - The response of a filter read off its pole-zero map.
% 1) |H| = |g| prod|e^jw - z_k| / prod|e^jw - p_k|: geometry vs polynomial evaluation
% 2) A pole at radius r: -3 dB width ~ 2(1-r) rad/sample, peak group delay r/(1-r),
% ringing time -1/ln(r) ~ 1/(1-r) samples
% 3) Linear-phase FIR: zeros on the unit circle or in reciprocal pairs (z, 1/z*)
% What to change: r_list, theta, the FIR length and cutoff.
% Expected output: geometric and polynomial |H| agree to ~1e-15; for r = 0.95 the
% width is ~0.102 rad/sample (2(1-r) = 0.1) and the peak delay ~18.5 samples for the
% pole pair (19 for a single pole); FIR zeros print as reciprocal pairs.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; clc;
rng(1);
%% Helpers (anonymous only)
% frequency response and group delay of b/a on a grid w (rad/sample), via
% tau = Re{ sum(n c_n e^-jwn) / sum(c_n e^-jwn) } for numerator minus denominator
freqresp = @(c, w) exp(-1j * w(:) * (0:numel(c)-1)) * c(:);
dfreq = @(c, w) exp(-1j * w(:) * (0:numel(c)-1)) * ((0:numel(c)-1)' .* c(:));
grpdel = @(b, a, w) real(dfreq(b, w) ./ freqresp(b, w)) - real(dfreq(a, w) ./ freqresp(a, w));
%% 1) Resonator: geometry vs polynomial
r = 0.95; theta = pi/4; % pole radius and angle
p = [r*exp(1j*theta); r*exp(-1j*theta)]; % conjugate pole pair, no zeros
a = real(poly(p)); b = 1;
w = linspace(0, pi, 4001)';
Hpoly = abs(freqresp(b, w) ./ freqresp(a, w));
Hgeo = 1 ./ (abs(exp(1j*w) - p(1)) .* abs(exp(1j*w) - p(2)));
fprintf('Resonator r = %.2f, theta = pi/4\n', r);
fprintf(' max |Hgeo - Hpoly| / max|H| = %.2e\n', max(abs(Hgeo - Hpoly)) / max(Hpoly));
fprintf(' |H| at DC = %.3f, at theta = %.2f, peak/DC = %.1f dB\n', Hpoly(1), 1/(abs(exp(1j*theta)-p(1))*abs(exp(1j*theta)-p(2))), 20*log10(max(Hpoly)/Hpoly(1)));
%% 2) Width, delay and ringing versus r
r_list = [0.8 0.9 0.95 0.99];
fprintf('\n r -3dB width 2(1-r) peak tau(pair) r/(1-r) -1/ln r\n');
for r = r_list
p = [r*exp(1j*theta); r*exp(-1j*theta)]; a = real(poly(p));
H = abs(freqresp(1, w) ./ freqresp(a, w)); [Hm, im] = max(H);
above = find(H >= Hm/sqrt(2)); width = w(max(above)) - w(min(above));
tau = grpdel(1, a, w);
fprintf(' %.2f %8.4f %7.4f %10.2f %8.2f %7.2f\n', r, width, 2*(1-r), tau(im), r/(1-r), -1/log(r));
end
%% 3) Linear-phase FIR zeros (21 taps, Kaiser beta = 5, cutoff 0.42 pi)
M = 20; wc = 0.42*pi; beta = 5;
k = 0:40; I0 = @(x) sum(((x(:)/2).^(2*k)) ./ (factorial(k).^2), 2); % modified Bessel I0 by series
n = (0:M)'; m = n - M/2;
hd = wc/pi * ones(M+1, 1); nz = m ~= 0; hd(nz) = sin(wc*m(nz)) ./ (pi*m(nz));
win = I0(beta*sqrt(1 - (2*n/M - 1).^2)) / I0(beta);
h = hd .* win;
z = roots(h);
fprintf('\nLinear-phase FIR: %d taps, %d zeros\n', M+1, numel(z));
fprintf(' zeros on the unit circle (| |z| - 1 | < 1e-6): %d\n', sum(abs(abs(z) - 1) < 1e-6));
zin = z(abs(z) < 1 - 1e-6); zout = z(abs(z) > 1 + 1e-6);
for i = 1:numel(zin)
[dmin, j] = min(abs(zout - 1/conj(zin(i))));
fprintf(' |z| = %.4f paired with |z| = %.4f product = %.6f (match err %.1e)\n', abs(zin(i)), abs(zout(j)), abs(zin(i))*abs(zout(j)), dmin);
end
fprintf(' symmetry check max|h - flip(h)| = %.1e, delay M/2 = %d samples\n', max(abs(h - flipud(h))), M/2);
%% Figures
figure('Name', 'Pole-zero geometry');
r = 0.95; p = [r*exp(1j*theta); r*exp(-1j*theta)]; a = real(poly(p));
subplot(2,2,1); tt = linspace(0, 2*pi, 400); plot(cos(tt), sin(tt), 'k-'); hold on;
plot(real(p), imag(p), 'rx', 'MarkerSize', 10, 'LineWidth', 2); plot(cos(theta), sin(theta), 'o', 'MarkerFaceColor', [0.8 0.6 0]);
plot([real(p(1)) cos(theta)], [imag(p(1)) sin(theta)], '-', 'Color', [0.8 0.6 0]); plot([real(p(2)) cos(theta)], [imag(p(2)) sin(theta)], '-', 'Color', [0.8 0.6 0]);
axis image; axis([-1.3 1.3 -1.3 1.3]); grid on; title('Resonator poles, r = 0.95'); xlabel('Re z'); ylabel('Im z');
subplot(2,2,2); plot(w/pi, 20*log10(abs(freqresp(1, w) ./ freqresp(a, w))), 'LineWidth', 1.5); grid on;
xlabel('\omega / \pi'); ylabel('|H| (dB)'); title('Magnitude');
subplot(2,2,3); hold on; cols = lines(4);
for i = 1:numel(r_list)
rr = r_list(i); pp = [rr*exp(1j*theta); rr*exp(-1j*theta)]; plot(w/pi, grpdel(1, real(poly(pp)), w), 'Color', cols(i,:), 'LineWidth', 1.3);
end
grid on; xlabel('\omega / \pi'); ylabel('\tau_g (samples)'); title('Group delay, r = 0.8 ... 0.99'); ylim([-5 100]);
subplot(2,2,4); plot(cos(tt), sin(tt), 'k-'); hold on; plot(real(z), imag(z), 'bo', 'LineWidth', 1.5);
axis image; axis([-2.5 2.5 -2.5 2.5]); grid on; title('Zeros of a 21-tap linear-phase FIR'); xlabel('Re z'); ylabel('Im z');
What to look for: the geometric \(\lvert H\rvert\) equals the polynomial evaluation to machine precision; for r = 0.95 the printed −3 dB width is close to 0.1 rad/sample and the peak delay close to 19 samples; the linear-phase FIR zeros print in reciprocal pairs (\(\lvert z\rvert\cdot\lvert 1/z\rvert = 1\)).
m09_fir_window_kaiser.mKaiser window design from the formulas vs Lawson IRLS equiripple vs order estimates · change A, transition · runtime ≈ 3 s
%% m09_fir_window_kaiser.m
% Module 09 - FIR design without a toolbox: Kaiser window formulas vs equiripple.
% 1) Kaiser: beta(A) and M = (A - 8)/(2.285 dw) for an RF anti-alias low-pass
% (fs = 40 MHz, pass <= 8 MHz, stop >= 9 MHz, A = 60 dB); measure what it achieves
% 2) Equiripple order estimates (Kaiser 1974, Herrmann-Rabiner-Chan 1973) for
% dp = 0.01, ds = 0.001, and an equiripple design by Lawson's iteratively
% reweighted least squares (IRLS), which needs only backslash
% What to change: fp, fst, A, dp, ds, the IRLS orders.
% Expected output: beta = 5.653, M = 145 -> about 59.6 dB; A = 61 -> M = 148, 60.9 dB;
% estimates 101.4 (both); IRLS fails at order 102 and meets dp <= 0.01, ds <= 0.001 from 104.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; clc;
rng(1);
%% Spec
fs = 40; % sampling rate (MHz)
fp = 8; fst = 9; % passband and stopband edges (MHz)
wp = 2*pi*fp/fs; ws = 2*pi*fst/fs; dw = ws - wp; % rad/sample
wgrid = linspace(0, pi, 16384)'; % dense evaluation grid
%% Helpers (anonymous only)
k = 0:40; I0 = @(x) sum(((x(:)/2).^(2*k)) ./ (factorial(k).^2), 2); % Bessel I0, series
kbeta = @(A) (A > 50) .* 0.1102 .* (A - 8.7) + (A >= 21 & A <= 50) .* (0.5842*max(A-21,0).^0.4 + 0.07886*(A-21));
kord = @(A, dw) (A - 8) ./ (2.285 * dw);
kwin = @(M, b) I0(b * sqrt(max(0, 1 - (2*(0:M)'/M - 1).^2))) / I0(b);
ideal = @(M, wc) (wc/pi) * ((0:M)' == M/2) + ((0:M)' ~= M/2) .* sin(wc*((0:M)' - M/2)) ./ (pi*((0:M)' - M/2) + ((0:M)' == M/2));
Hmag = @(h) abs(exp(-1j * wgrid * (0:numel(h)-1)) * h(:));
ripple_p = @(H) max(abs(H(wgrid <= wp) - 1));
ripple_s = @(H) max(H(wgrid >= ws));
%% 1) Kaiser window design
fprintf('Kaiser window, pass %g MHz, stop %g MHz at fs = %g MHz\n', fp, fst, fs);
designs = {};
for A = [60 61]
b = kbeta(A); Mf = kord(A, dw); M = ceil(Mf);
h = ideal(M, (wp + ws)/2) .* kwin(M, b);
H = Hmag(h); dp_ = ripple_p(H); ds_ = ripple_s(H);
fprintf(' A = %d dB: beta = %.3f, M = %.1f -> %d (%d taps); measured stopband %.2f dB, passband dev %.4f\n', ...
A, b, Mf, M, M+1, 20*log10(ds_), dp_);
designs{end+1} = h; %#ok<SAGROW>
end
hK = designs{2}; MK = numel(hK) - 1;
fprintf(' delay M/2 = %d samples = %.3f us = %.2f mm (c/2 = 770 m/s)\n', MK/2, MK/2/fs, MK/2/fs*0.770);
%% 2) Equiripple order estimates
dp = 0.01; ds = 0.001; df = (fst - fp)/fs;
Nk = (-20*log10(sqrt(dp*ds)) - 13) / (14.6*df);
a6 = [0.005309 0.07114 -0.4761 0.00266 0.5941 0.4278]; lp = log10(dp); ls = log10(ds);
Dinf = ls*(a6(1)*lp^2 + a6(2)*lp + a6(3)) - (a6(4)*lp^2 + a6(5)*lp + a6(6));
fK = 11.01217 + 0.51244*(lp - ls); Nh = Dinf/df - fK*df;
fprintf('\nEquiripple estimates (dp = %.3f, ds = %.4f, df = %.4f): Kaiser %.1f, Herrmann %.1f\n', dp, ds, df, Nk, Nh);
fprintf(' A/22 rule: %.0f taps for A = %d dB\n', -20*log10(ds)/(22*df), round(-20*log10(ds)));
%% 3) Equiripple by Lawson IRLS (type I linear phase, M even)
% A(w) = sum_k a_k cos(k w); weighted LS on a grid of the two bands, weights updated
% as u <- u .* |E| so the error envelope flattens toward the minimax solution.
Wst = dp/ds; % stopband weight
fprintf('\nIRLS equiripple, weight %.0f in the stopband:\n', Wst);
hE = [];
for M = [102 104 106 108]
L = M/2; ng = 16*M;
gp = linspace(0, wp, round(ng*wp/(pi - dw)))'; gs = linspace(ws, pi, round(ng*(pi - ws)/(pi - dw)))';
g = [gp; gs]; d = [ones(size(gp)); zeros(size(gs))]; W = [ones(size(gp)); Wst*ones(size(gs))];
C = cos(g * (0:L)); u = W;
for it = 1:80
a = (C .* u) \ (d .* u); % weighted least squares
e = abs(W .* (C*a - d));
u = u .* sqrt(e + 1e-14); u = u / sqrt(mean(u.^2));
end
h = [flipud(a(2:end))/2; a(1); a(2:end)/2];
H = Hmag(h); dp_ = ripple_p(H); ds_ = ripple_s(H);
ok = dp_ <= dp && ds_ <= ds;
fprintf(' M = %d: passband dev %.5f, stopband %.6f (%.2f dB) -> %s\n', M, dp_, ds_, 20*log10(ds_), char('fails' .* ~ok + 'meets' .* ok));
if M == 106, hE = h; end
end
fprintf(' (true Remez minimum for this spec: order 104; equiripple with dp = ds = 0.001 needs 132)\n');
%% Figures
f = wgrid/(2*pi)*fs;
figure('Name', 'Kaiser vs equiripple');
subplot(2,1,1); plot(f, 20*log10(Hmag(hK) + 1e-12), 'b', 'LineWidth', 1); hold on;
plot([fst fs/2], [-60 -60], 'r--', [fp fp], [-120 5], 'r:', [fst fst], [-120 5], 'r:');
ylim([-110 5]); grid on; xlabel('frequency (MHz)'); ylabel('|H| (dB)'); title(sprintf('Kaiser window, M = %d, beta = %.2f', MK, kbeta(61)));
subplot(2,1,2); plot(f, 20*log10(Hmag(hE) + 1e-12), 'Color', [0.8 0.2 0.4], 'LineWidth', 1); hold on;
plot([fst fs/2], [-60 -60], 'r--', [fp fp], [-120 5], 'r:', [fst fst], [-120 5], 'r:');
ylim([-110 5]); grid on; xlabel('frequency (MHz)'); ylabel('|H| (dB)'); title('Equiripple (IRLS), M = 106, dp = 0.01, ds = 0.001');
What to look for: β = 5.653 and M = 145 for the 60 dB spec with a measured stopband of about 59.6 dB, then 60.9 dB with A = 61 (M = 148); Kaiser and Herrmann estimates of 101.4 for the \(\delta_p = 0.01\) equiripple and an IRLS design that fails at order 102 and meets the spec from order 104 (\(\delta_s = 0.000992\)), matching the Remez minimum.
m09_iir_bilinear_butter_cheby.mButterworth and Chebyshev I by analog prototype + prewarped bilinear transform, order formulas, biquads vs direct-form quantization · runtime ≈ 3 s
%% m09_iir_bilinear_butter_cheby.m
% Module 09 - IIR design without a toolbox: analog prototype + prewarped bilinear
% transform, order formulas, biquads, and why direct form is fragile.
% 1) Orders for the RF anti-alias spec (fs = 40 MHz, pass 8 MHz with Ap = 0.17 dB,
% stop 9 MHz with As = 60 dB): Butterworth, Chebyshev, elliptic (via ellipke)
% 2) Butterworth and Chebyshev I designs from their analog poles, bilinear with
% prewarping, implemented as second-order sections; check the edges land exactly
% 3) Coefficient quantization: 8th-order Butterworth at fc = 0.02 fs in direct form
% (16- and 24-bit words) vs biquads; 12th-order at 0.01 fs in double precision
% What to change: fp, fst, Ap, As; the word lengths; N and fc of the fragile filter.
% Expected output: orders 53 / 16 / 8; |H(8 MHz)| = -0.174 dB and |H(9 MHz)| <= -60 dB;
% direct-form max pole radius 1.299 (16 bit) and 1.166 (24 bit), biquads 0.9758;
% the N = 12 direct-form impulse response diverges, the biquad version decays.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; clc;
rng(1);
%% Spec and order formulas
fs = 40; fp = 8; fst = 9; % MHz
Ap = 20*log10(1.01/0.99); As = 60; % dB: the +-0.01 passband of the FIR spec, 60 dB stopband
ep = sqrt(10^(Ap/10) - 1); es = sqrt(10^(As/10) - 1);
Wp = 2*tan(pi*fp/fs); Ws = 2*tan(pi*fst/fs); % prewarped edges (T = 1, 2/T = 2)
Nb = ceil(log(es/ep) / log(Ws/Wp));
Nc = ceil(acosh(es/ep) / acosh(Ws/Wp));
k = Wp/Ws; k1 = ep/es; % ellipke takes m = k^2
[K, ~] = ellipke(k^2); [Kp, ~] = ellipke(1 - k^2);
[K1, ~] = ellipke(k1^2); [K1p, ~] = ellipke(1 - k1^2);
Ne = ceil(K*K1p / (Kp*K1));
fprintf('Spec: pass %g MHz (Ap = %.3f dB), stop %g MHz (As = %g dB), fs = %g MHz\n', fp, Ap, fst, As, fs);
fprintf(' prewarped edges: Wp = %.4f, Ws = %.4f (selectivity %.4f)\n', Wp, Ws, Ws/Wp);
fprintf(' orders: Butterworth %d, Chebyshev I/II %d, elliptic %d\n', Nb, Nc, Ne);
%% Helpers (anonymous only)
bil = @(s) (2 + s) ./ (2 - s); % bilinear map s -> z (2/T = 2)
evalz = @(zz, pp, g, w) g * prod(exp(1j*w(:).') - zz(:), 1) ./ prod(exp(1j*w(:).') - pp(:), 1);
wgrid = linspace(0, pi, 8001); f = wgrid/(2*pi)*fs;
%% Butterworth (passband edge exact) and Chebyshev I from analog poles
for fam = 1:2
if fam == 1
N = Nb; Wc = Wp / ep^(1/N); kk = 1:N;
pa = Wc * exp(1j*pi*(2*kk + N - 1)/(2*N)); name = 'Butterworth';
else
N = Nc; mu = asinh(1/ep)/N; kk = 1:N; th = pi*(2*kk - 1)/(2*N);
pa = Wp * (-sinh(mu)*sin(th) + 1j*cosh(mu)*cos(th)); name = 'Chebyshev I';
end
pd = bil(pa); zd = -ones(1, N); % N analog zeros at infinity -> z = -1
dc = 1; if fam == 2 && mod(N, 2) == 0, dc = 1/sqrt(1 + ep^2); end
g = real(dc * prod(1 - pd) / prod(1 - zd)); % set the DC gain
H = abs(evalz(zd, pd, g, wgrid));
% biquads: conjugate pole pairs with a double zero at z = -1
pu = pd(imag(pd) > 1e-12); pr = pd(abs(imag(pd)) <= 1e-12);
y = [1 zeros(1, 2999)]; gs = g;
for i = 1:numel(pu)
a2 = [1, -2*real(pu(i)), abs(pu(i))^2]; y = filter([1 2 1], a2, y);
end
for i = 1:numel(pr), y = filter([1 1], [1 -pr(i)], y); end
y = gs * y;
Hs = abs(exp(-1j*wgrid(:)*(0:numel(y)-1)) * y(:)); % response of the SOS impulse response
fprintf('\n%s, N = %d (%d biquads):\n', name, N, ceil(N/2));
fprintf(' |H| at %g MHz = %.3f dB (target -%.3f), at %g MHz = %.2f dB (target <= -%g)\n', fp, 20*log10(abs(evalz(zd, pd, g, 2*pi*fp/fs))), Ap, fst, 20*log10(abs(evalz(zd, pd, g, 2*pi*fst/fs))), As);
fprintf(' max pole radius %.5f; SOS impulse response vs zpk: max |dH| = %.1e\n', max(abs(pd)), max(abs(Hs(:) - H(:))));
if fam == 1, Hb = H; else, Hc = H; end
end
%% Coefficient quantization: direct form vs biquads
fprintf('\nCoefficient quantization (Butterworth low-pass, -3 dB at fc):\n');
N = 8; fc = 0.02; Wc = 2*tan(pi*fc); kk = 1:N;
pd = bil(Wc * exp(1j*pi*(2*kk + N - 1)/(2*N)));
a = real(poly(pd));
ib = ceil(log2(max(abs(a)))) + 1; % integer bits incl. sign
fprintf(' N = %d, fc = %.2f fs: true max pole radius %.4f, max |a_k| = %.1f (%d integer bits)\n', N, fc, max(abs(pd)), max(abs(a)), ib);
for W = [16 24]
fr = W - ib; aq = round(a * 2^fr) / 2^fr;
rq = max(abs(roots(aq)));
pu = pd(imag(pd) > 0); r2 = 0;
for i = 1:numel(pu)
sec = [1, -2*real(pu(i)), abs(pu(i))^2]; secq = round(sec * 2^(W-2)) / 2^(W-2); r2 = max(r2, max(abs(roots(secq))));
end
fprintf(' %d-bit words: direct form max pole radius %.4f %s; biquads (%d fractional bits) %.4f\n', W, rq, char('(unstable)' .* (rq > 1) + ' ' .* (rq <= 1)), W-2, r2);
if W == 16, pq16 = roots(aq); end
end
N = 12; fc = 0.01; Wc = 2*tan(pi*fc); kk = 1:N;
pd = bil(Wc * exp(1j*pi*(2*kk + N - 1)/(2*N))); a = real(poly(pd)); g = real(prod(1 - pd) / 2^N); b = g * real(poly(-ones(1, N)));
x = [1 zeros(1, 39999)]; yd = filter(b, a, x);
ys = x; pu = pd(imag(pd) > 0); for i = 1:numel(pu), ys = filter([1 2 1], [1, -2*real(pu(i)), abs(pu(i))^2], ys); end; ys = g * ys;
fprintf(' N = 12, fc = 0.01 fs, double precision: direct-form tail max |h| = %.3g (NaN/Inf = diverged); biquad tail max |h| = %.3g\n', max(abs(yd(end-1999:end))), max(abs(ys(end-1999:end))));
%% Figures
figure('Name', 'IIR by bilinear transform');
subplot(2,1,1); plot(f, 20*log10(Hb + 1e-15), 'k', f, 20*log10(Hc + 1e-15), 'b', 'LineWidth', 1.2); hold on;
plot([fst fs/2], [-As -As], 'r--', [fp fp], [-120 5], 'r:', [fst fst], [-120 5], 'r:');
ylim([-100 5]); xlim([0 20]); grid on; xlabel('frequency (MHz)'); ylabel('|H| (dB)'); title(sprintf('Butterworth N = %d (black), Chebyshev I N = %d (blue)', Nb, Nc));
subplot(2,1,2); N = 8; Wc = 2*tan(pi*0.02); kk = 1:N; pt = bil(Wc * exp(1j*pi*(2*kk + N - 1)/(2*N)));
tt = linspace(-0.6, 0.6, 200); plot(cos(tt), sin(tt), 'k-'); hold on;
plot(real(pt), imag(pt), 'bx', 'MarkerSize', 9, 'LineWidth', 1.5); plot(real(pq16), imag(pq16), 'rx', 'MarkerSize', 9, 'LineWidth', 1.5);
axis image; axis([0.4 1.45 -0.5 0.5]); grid on; xlabel('Re z'); ylabel('Im z'); title('N = 8 at 0.02 fs: true poles (blue), 16-bit direct form (red)');
What to look for: orders 53 (Butterworth) and 16 (Chebyshev I) for the 8/9 MHz, 0.17/60 dB spec; passband edges landing exactly at 8 MHz after prewarping; direct-form pole radius 1.299 at 16 bits and 1.166 at 24 bits for the 8th-order filter at \(0.02\,f_s\), unchanged poles for biquads; a diverging direct-form output for N = 12 at \(0.01\,f_s\).
m09_fir_vs_iir_pulse.m5 MHz RF pulse through an equiripple FIR and a Chebyshev I band-pass; group delay, envelope shift, forward–backward · runtime ≈ 1 s
%% m09_fir_vs_iir_pulse.m
% Module 09 - The same RF band-pass spec met by an equiripple FIR and by a
% Chebyshev I IIR: order, delay, group-delay dispersion and pulse distortion.
% Spec: pass 3-7 MHz (0.5 dB), stop <= 2 and >= 8 MHz (40 dB), fs = 40 MHz.
% FIR: order 72 equiripple by Lawson IRLS (linear phase, delay 36 samples)
% IIR: Chebyshev I prototype N = 8 -> band-pass order 16, prewarped bilinear, biquads
% Pulse: Gaussian-enveloped 5 MHz, 40 % (-6 dB) fractional bandwidth.
% What to change: f0, bw, the IIR prototype order, forward-backward on/off.
% Expected output: FIR group delay 36.0 samples flat; IIR group delay ~21 samples at
% 5 MHz rising to >100 near the 3 MHz edge; rho(FIR) ~ 1.000, rho(IIR causal) clearly
% lower than rho(IIR magnitude only); forward-backward IIR aligned with zero shift.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; clc;
rng(1);
fs = 40; % MHz
W = @(f) 2*pi*f/fs; % MHz -> rad/sample
Ap = 0.5; As = 40;
dp = (10^(Ap/20) - 1)/(10^(Ap/20) + 1); ds = 10^(-As/20);
%% FIR: type I equiripple via IRLS, order 72
M = 72; L = M/2; ng = 2000;
tot = W(2) + (W(7) - W(3)) + (pi - W(8)); % grid per band, edges included exactly
g = [linspace(0, W(2), round(ng*W(2)/tot))'; linspace(W(3), W(7), round(ng*(W(7)-W(3))/tot))'; linspace(W(8), pi, round(ng*(pi-W(8))/tot))'];
d = double(g >= W(3) - 1e-12 & g <= W(7) + 1e-12); Wt = ones(size(g)); Wt(d == 0) = dp/ds;
C = cos(g * (0:L)); u = Wt;
for it = 1:80
a = (C .* u) \ (d .* u); e = abs(Wt .* (C*a - d)); u = u .* sqrt(e + 1e-14); u = u / sqrt(mean(u.^2));
end
h = [flipud(a(2:end))/2; a(1); a(2:end)/2];
%% IIR: Chebyshev I prototype -> band-pass -> bilinear, as biquads
N = 8; ep = sqrt(10^(Ap/10) - 1); mu = asinh(1/ep)/N; th = pi*(2*(1:N) - 1)/(2*N);
pl = -sinh(mu)*sin(th) + 1j*cosh(mu)*cos(th); % prototype poles, passband edge 1 rad/s
W1 = 2*tan(W(3)/2); W2 = 2*tan(W(7)/2); W0 = sqrt(W1*W2); B = W2 - W1;
pb = [pl*B/2 + sqrt((pl*B/2).^2 - W0^2), pl*B/2 - sqrt((pl*B/2).^2 - W0^2)]; % s -> (s^2 + W0^2)/(B s)
pz = (2 + pb) ./ (2 - pb); % bilinear (T = 1)
pu = pz(imag(pz) > 1e-12); % one of each conjugate pair (16 poles -> 8)
% each biquad: zeros at z = +1 and z = -1 (analog zeros at s = 0 and s = inf)
evalH = @(w) prod((exp(-2j*w(:)) - 1) * ones(1, numel(pu)) ./ (1 - 2*real(pu(:).') .* exp(-1j*w(:)) + abs(pu(:).').^2 .* exp(-2j*w(:))), 2);
wg = linspace(0, pi, 8001)'; fg = wg/(2*pi)*fs;
Hraw = evalH(wg); gIIR = 1/max(abs(Hraw(fg >= 3 & fg <= 7))); % Chebyshev I: passband peak = 1
%% Responses and group delays
Hfir = exp(-1j*wg*(0:M)) * h;
Hiir = gIIR * Hraw;
gdf = @(H) -gradient(unwrap(angle(H)), wg(2) - wg(1));
gdFIR = gdf(Hfir); gdIIR = gdf(Hiir);
sel = fg >= 3 & fg <= 7;
fprintf('FIR order %d: passband dev %.4f (spec %.4f), stopband %.2f dB (spec -%g)\n', M, max(abs(abs(Hfir(sel)) - 1)), dp, 20*log10(max(abs(Hfir(fg <= 2 | fg >= 8)))), As);
fprintf('IIR Chebyshev I band-pass order %d: passband loss %.3f dB, stopband %.2f dB\n', 2*N, -20*log10(min(abs(Hiir(sel)))), 20*log10(max(abs(Hiir(fg <= 2 | fg >= 8)))));
[~, i5] = min(abs(fg - 5));
fprintf('Group delay: FIR %.1f samples everywhere in 3-7 MHz (%.2f to %.2f)\n', gdFIR(i5), min(gdFIR(sel)), max(gdFIR(sel)));
fprintf(' IIR %.1f samples at 5 MHz, %.1f to %.1f over 3-7 MHz\n', gdIIR(i5), min(gdIIR(sel)), max(gdIIR(sel)));
fprintf('Multiplies per output: FIR %d (folded), IIR %d (5 per biquad)\n', M/2 + 1, 5*numel(pu));
%% Pulse through both
Ns = 1024; n0 = 200; f0 = 5; bw = 0.40;
st = sqrt(2*log(2)) / (pi*bw*f0) * fs; % envelope sigma (samples)
n = (0:Ns-1)'; x = exp(-(n - n0).^2/(2*st^2)) .* cos(2*pi*f0*(n - n0)/fs);
yF = filter(h, 1, x);
yI = x; for i = 1:numel(pu), yI = filter([1 0 -1], [1, -2*real(pu(i)), abs(pu(i))^2], yI); end; yI = gIIR * yI;
yB = flipud(yI); % forward-backward: reverse the causal output,
for i = 1:numel(pu), yB = filter([1 0 -1], [1, -2*real(pu(i)), abs(pu(i))^2], yB); end % filter again,
yB = flipud(gIIR * yB); % and reverse back: |H|^2, zero phase
% magnitude-only version (|H| applied with zero phase) via FFT
X = fft(x); kf = (0:Ns-1)'; wk = 2*pi*min(kf, Ns - kf)/Ns; Hm = abs(gIIR * evalH(wk)); yM = real(ifft(X .* Hm));
% envelope by FFT analytic signal
anal = @(v) ifft(fft(v) .* [1; 2*ones(Ns/2 - 1, 1); 1; zeros(Ns/2 - 1, 1)]);
env = @(v) abs(anal(v));
% correlation coefficient after best integer alignment (circular is safe: signals are zero-padded)
rho = @(y) max(real(ifft(fft(y) .* conj(fft(x))))) / (norm(y) * norm(x));
[~, kF] = max(env(yF)); [~, kI] = max(env(yI)); [~, kB] = max(env(yB)); [~, kx] = max(env(x));
fprintf('\nPulse f0 = %g MHz, %d %% bandwidth:\n', f0, round(100*bw));
fprintf(' envelope peak delay: FIR %d, IIR %d, IIR forward-backward %d samples\n', kF - kx, kI - kx, kB - kx);
fprintf(' rho: FIR %.3f | IIR causal %.3f | IIR magnitude only %.3f | IIR forward-backward %.3f\n', rho(yF), rho(yI), rho(yM), rho(yB));
%% Figures
figure('Name', 'FIR vs IIR band-pass');
subplot(3,1,1); plot(fg, 20*log10(abs(Hfir) + 1e-12), 'b', fg, 20*log10(abs(Hiir) + 1e-12), 'r', 'LineWidth', 1.1); hold on;
plot([0 2], [-As -As], 'k--', [8 20], [-As -As], 'k--'); ylim([-80 5]); grid on; xlabel('frequency (MHz)'); ylabel('|H| (dB)'); title('Equiripple FIR (blue) vs Chebyshev I (red)');
gdFp = gdFIR; gdFp(fg < 2.5 | fg > 7.5) = NaN; % FIR phase jumps at its stopband zeros
subplot(3,1,2); plot(fg, gdFp, 'b', fg, gdIIR, 'r', 'LineWidth', 1.1); xlim([1 9]); ylim([0 150]); grid on; xlabel('frequency (MHz)'); ylabel('\tau_g (samples)'); title('Group delay');
subplot(3,1,3); t = (n - n0)/fs; plot(t, circshift(x, 36), 'k', t, yF, 'b', t, yI, 'r', t, yB, 'm--', 'LineWidth', 1); xlim([-1 3]); grid on;
xlabel('time (\mus)'); ylabel('amplitude'); title('Input (shifted 36 samples), FIR, IIR causal, IIR forward-backward');
What to look for: a constant FIR group delay of M/2 samples against an IIR delay that rises steeply toward the 3 and 7 MHz edges; correlation ρ after alignment lower for the causal IIR than for the zero-phase version of the same filter, which isolates the phase damage; forward–backward output aligned with the input.
m09_wall_filters_ensemble.mWall filters as matrices on N = 10 pulses: FIR, IIR with zero/step/projection initialization, polynomial regression; power responses and clutter transients · runtime ≈ 1 s
%% m09_wall_filters_ensemble.m
% Module 09 - Doppler wall (clutter) filters on a short slow-time ensemble, judged as
% matrices y = A x. Power response |H(f)|^2 = ||A e_f||^2 / N_out (Bjaerum et al. 2002).
% FIR : order K windowed-sinc high-pass (Hamming), N - K valid outputs
% IIR : 2nd-order Butterworth high-pass (bilinear), zero / step / projection init
% Regr.: polynomial regression of degree p, A = I - Q Q'
% What to change: N, fc, K, p, the tissue Doppler ft and the clutter level.
% Expected output (N = 10, fc = 0.1 PRF): at DC the zero-state IIR gives only -12.6 dB;
% at 0.01 PRF step init -32 dB, projection -43 dB, regression p = 2 -63 dB;
% regression -3 dB cutoffs ~ (0.44 + 0.41 p) PRF/N.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; clc;
rng(1);
N = 10; % ensemble length (pulses)
fc = 0.1; % IIR / FIR cutoff, normalized to the PRF
K = 4; % FIR order
p = 2; % regression degree
%% Filters
% IIR: Butterworth low-pass prototype -> high-pass (s -> Wc/s) -> bilinear (2/T = 2)
no = 2; Wc = 2*tan(pi*fc); kk = 1:no;
pl = exp(1j*pi*(2*kk + no - 1)/(2*no)); % prototype poles (Wc = 1)
ph = Wc ./ pl; % high-pass poles; zeros at s = 0
pz = (2 + ph) ./ (2 - ph); a = real(poly(pz)); b = real(poly(ones(1, no)));
b = b * abs(sum(a .* (-1).^(0:no))) / abs(sum(b .* (-1).^(0:no))); % unit gain at Nyquist
% steady-state initial state for a unit step (direct form II transposed, as filter() uses)
nz = numel(a) - 1; Mz = eye(nz); for i = 1:nz, Mz(i,1) = Mz(i,1) + a(i+1); if i < nz, Mz(i,i+1) = Mz(i,i+1) - 1; end; end
zi = Mz \ (b(2:end) - a(2:end)*b(1))';
% FIR: windowed-sinc high-pass, low-pass normalized to unit DC gain (exact zero at DC)
nf = (0:K)'; mf = nf - K/2; lp = 2*fc*ones(K+1, 1); q = mf ~= 0; lp(q) = sin(2*pi*fc*mf(q)) ./ (pi*mf(q));
lp = lp .* (0.54 - 0.46*cos(2*pi*nf/K)); hfir = -lp / sum(lp); hfir(K/2+1) = hfir(K/2+1) + 1;
%% Matrices
I = eye(N);
F = zeros(N); S = zeros(N); for k = 1:N, F(:,k) = filter(b, a, I(:,k)); S(:,k) = filter(b, a, I(:,k), zi*I(1,k)); end
G = zeros(N, nz); for i = 1:nz, e0 = zeros(nz, 1); e0(i) = 1; G(:,i) = filter(b, a, zeros(N,1), e0); end
P = (I - G*pinv(G)) * F; % projection initialization
t = linspace(-1, 1, N)'; V = zeros(N, p+1); for d = 0:p, V(:,d+1) = t.^d; end; [Q, ~] = qr(V, 0);
R = I - Q*Q'; % polynomial regression
Af = zeros(N-K, N); for i = 1:N-K, Af(i, i:i+K) = fliplr(hfir'); end % FIR, valid outputs only
presp = @(A, f) sum(abs(A * exp(2j*pi*(0:N-1)'*f(:)')).^2, 1) / size(A, 1);
fr = [0 0.01 0.02 0.05 0.25];
names = {'FIR K=4', 'IIR zero init', 'IIR step init', 'IIR projection', 'regression p=2'};
mats = {Af, F, S, P, R};
fprintf('Power response (dB), N = %d, fc = %.2f PRF\n', N, fc);
fprintf('%-16s %6s %7s %7s %7s %7s Nout\n', 'filter', 'f=0', '0.01', '0.02', '0.05', '0.25');
for i = 1:numel(mats)
r = 10*log10(presp(mats{i}, fr) + 1e-30); r(r < -99) = -99;
fprintf('%-16s %6.1f %7.1f %7.1f %7.1f %7.1f %d\n', names{i}, r, size(mats{i}, 1));
end
%% Regression cutoff versus degree and ensemble length
ff = linspace(0, 0.5, 20001);
fprintf('\nRegression -3 dB cutoff x N / PRF (approx. 0.44 + 0.41 p):\n N p=0 p=1 p=2 p=3\n');
for NN = [8 10 12 16 32]
tt = linspace(-1, 1, NN)'; row = zeros(1, 4);
for pp = 0:3
Vv = zeros(NN, pp+1); for d = 0:pp, Vv(:,d+1) = tt.^d; end; [Qq, ~] = qr(Vv, 0);
resp = 1 - sum(abs(Qq' * exp(2j*pi*(0:NN-1)'*ff)).^2, 1) / NN;
row(pp+1) = ff(find(resp >= 0.5, 1)) * NN;
end
fprintf(' %4d %5.3f %5.3f %5.3f %5.3f\n', NN, row);
end
fprintf(' N = 10, p = 2 at PRF 4 kHz, f0 = 5 MHz: cutoff %.0f Hz = %.1f cm/s\n', 1.267*4000/10, 100*1540*1.267*4000/10/(2*5e6));
%% Clutter + blood ensemble
ft = 0.01; n = (0:N-1)';
c = 100 * exp(2j*pi*ft*n); % clutter, 40 dB above blood
bl = exp(2j*pi*0.25*n + 0.7j); % blood at 0.25 PRF
fprintf('\nClutter-to-blood ratio after filtering (input +40.0 dB), tissue at %.2f PRF:\n', ft);
for i = 1:numel(mats)
A = mats{i}; fprintf(' %-16s %6.1f dB\n', names{i}, 10*log10(sum(abs(A*c).^2) / sum(abs(A*bl).^2)));
end
%% Figures
figure('Name', 'Wall filters');
subplot(2,1,1); fg = linspace(0, 0.5, 400); cols = [0.7 0.45 0; 0.5 0.5 0.5; 0.85 0.65 0.1; 0.8 0.2 0.4; 0.1 0.35 0.8];
hold on; for i = 1:numel(mats), plot(fg, 10*log10(presp(mats{i}, fg) + 1e-12), 'Color', cols(i,:), 'LineWidth', 1.3); end
ylim([-80 5]); grid on; xlabel('f / PRF'); ylabel('power response (dB)'); title(sprintf('N = %d, f_c = %.2f PRF', N, fc)); legend(names, 'Location', 'southeast');
subplot(2,1,2); hold on; x = c + bl;
for i = 1:numel(mats)
y = mats{i} * x; off = N - size(mats{i}, 1); plot((0:numel(y)-1) + off, 20*log10(abs(y)), 'o-', 'Color', cols(i,:), 'LineWidth', 1.2);
end
plot([0 N-1], [0 0], 'k:'); ylim([-40 45]); grid on; xlabel('slow-time sample n'); ylabel('|y[n]| (dB)'); title('Clutter (40 dB, 0.01 PRF) + blood (0 dB, 0.25 PRF)');
What to look for: the table of power responses at 0, 0.01, 0.02, 0.05 and 0.25 PRF: −12.6 dB at DC for the zero-state IIR, −32 dB (step) and −43 dB (projection) at 0.01 PRF, −63 dB for degree-2 regression; regression −3 dB cutoffs \(\approx (0.44 + 0.41p)\cdot\mathrm{PRF}/N\).
09.13Self-check
Q1A pole pair sits at radius 0.9 and angle 0.3π, \(f_s\) = 40 MHz. Roughly how wide (−3 dB) is the resonance, and how long does it ring?
\(2(1 - r) = 0.2\) rad/sample → \(0.2/(2\pi) \times 40\) MHz = 1.27 MHz. Ringing \(\approx 1/(1 - r) = 10\) samples (exactly \(-1/\ln 0.9 = 9.5\)) = 0.25 µs. At θ = 0.3π the conjugate pole is far enough that the approximation holds.
Q2Why can you not build a linear-phase high-pass with an even number of symmetric taps?
Even length means odd order M; symmetry gives \(H(-1) = \sum_n h[n](-1)^n = -H(-1)\), so \(H(-1) = 0\): a forced zero at Nyquist. Use odd length (type I) or antisymmetric taps with even length (type IV, zero at DC is fine for a high-pass).
Q3Kaiser design for 70 dB with a transition of \(0.02\,f_s\): β and M?
\(\beta = 0.1102 \times (70 - 8.7) = 6.76\). \(\Delta\omega = 2\pi \times 0.02 = 0.1257\) rad/sample, \(M = 62/(2.285 \times 0.1257) = 216\). Expect about ±0.5 dB from the formula; add 1 dB of margin if the spec is hard.
Q4Analog spec: 1 dB ripple up to \(\Omega_p\), 40 dB from \(2\Omega_p\). Butterworth and Chebyshev orders?
\(\varepsilon_p = 0.509\), \(\varepsilon_s = 99.99\), ratio 196.5. Butterworth: \(\ln 196.5/\ln 2 = 7.62 \to 8\). Chebyshev: \(\operatorname{acosh} 196.5/\operatorname{acosh} 2 = 5.97/1.317 = 4.54 \to 5\). Elliptic needs 4.
Q5You run a 2nd-order Butterworth low-pass (\(f_c\) = 1 MHz) forward and backward. Attenuation at 1 MHz? Delay?
−6 dB (\(\lvert H\rvert^2 = 0.5\) at the original −3 dB point), zero delay at all frequencies. The new −3 dB point is \((\sqrt{2} - 1)^{1/4} = 0.80\) of 1 MHz.
Q6Decimate RF from 40 to 10 MHz while keeping 0–3 MHz (after mixing to baseband) alias-free at 60 dB. Stopband edge, rough taps, cost?
\(f_{\text{stop}} = f_s/D - f_{\text{pass}}\) = 10 − 3 = 7 MHz, so \(\Delta f = 4/40 = 0.1\). Taps \(\approx 60/(22 \times 0.1) \approx 27\). Polyphase: 27/4 ≈ 7 multiplies per input sample (per real channel). Aliases landing in 3–7 MHz are removed later or ignored.
Q7Color flow with N = 12, PRF = 2 kHz, \(f_0\) = 3 MHz and a degree-1 regression wall filter. Where is the −3 dB cutoff, in Hz and cm/s?
\(0.86\,\mathrm{PRF}/N = 0.86 \times 2000/12\) = 143 Hz. \(v = c f/(2 f_0) = 1540 \times 143/(6 \times 10^6)\) = 3.7 cm/s at zero beam-to-flow angle.
09.14Figures from the literature
Published figures to hold the designs of this page against: in vivo wall-filter results and the classic design-theory plots. Links only, no images reproduced; each figure number was checked against an author-posted or open copy.
Regression versus eigen-based wall filters on in vivo color-flow ensembles
A. C. H. Yu, L. Løvstakken, "Eigen-based clutter filter design for ultrasound color flow imaging: a review," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 57(5), 1096–1111, 2010. Figs. 5 and 6 (also Figs. 9 and 10). doi:10.1109/TUFFC.2010.1521 (open copy: HKU Scholars Hub)
Doppler power maps from an intra-operative coronary frame in systole (Fig. 5, N = 10) and a thyroid nodule (Fig. 6, N = 12), filtered by polynomial regression with and without down-mixing and by three eigen-filter variants. These are the ensemble lengths of Demo 09.6: look at where regression leaves tissue flash when the clutter moves, the failure mode that motivates the data-adaptive filters of Module 06. Figs. 9 and 10 track filtered power and velocity frame by frame.
IIR wall filter: what the initialization does to the response on short ensembles
S. Bjærum, H. Torp, "Statistical evaluation of clutter filters in color flow imaging," Ultrasonics 38, 376–380, 2000. Figs. 1, 2 and 3. Open copy: NTNU (author page)
Fig. 2 shows the frequency responses of a fourth-order Chebyshev high-pass with different initializations, Fig. 1 those of fifth-order FIR filters and Fig. 3 polynomial regression filters of order 0–3 for N = 9, all defined as in (09.11). Compare Fig. 2 with rule 13's figure (zero, step, projection) and Fig. 3 with the re-plot below; the paper then evaluates the filters on thyroid recordings.
The regression filter in closed form, and its frequency response
H. Torp, "Clutter rejection filters in color flow imaging: a theoretical approach," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 44(2), 417–424, 1997. Figs. 2 and 3. Author copy: NTNU (author page)
Fig. 2 shows the power spectra of the first Legendre basis functions and the response of a third-order regression filter; Fig. 3 the responses of regression filters of order 0 to 4 for N = 10. The figure below is re-computed from the same definition, \(1 - \lVert\mathbf{Q}^H\mathbf{e}_f\rVert^2/N\) with Legendre (orthonormal polynomial) columns, so the curves should coincide with the published ones.
Extraripple and equiripple low-pass filters: the alternation theorem on a plot
L. R. Rabiner, O. Herrmann, "On the design of optimum FIR lowpass filters with even impulse response duration," IEEE Trans. Audio Electroacoust. AU-21(4), 329–336, 1973. Fig. 3 and Fig. 7. Author reprint: UCSB (Rabiner reprints)
Fig. 3 puts an extraripple design (seven error maxima) beside an equiripple one (six maxima) at the same short length, which makes the alternation count of Section 09.3 visible; Fig. 7 shows how the optimum moves as the ripple ratio is scaled, ending at an equiripple filter with M/2 + 1 peaks. The same authors' companion paper (Herrmann, Rabiner and Chan, 1973) turned such designs into the \(D_\infty(\delta_p, \delta_s)\) length formula quoted after (09.4).
09.15Go deeper
- A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson, 2010. Chapters 5–7: pole-zero analysis, linear phase types, window and optimal FIR design, bilinear transform.
- J. F. Kaiser, "Nonrecursive digital filter design using the I0-sinh window function," Proc. IEEE Int. Symp. Circuits and Systems, 20–23, 1974. The β and order formulas.
- J. H. McClellan, T. W. Parks and L. R. Rabiner, "A computer program for designing optimum FIR linear phase digital filters," IEEE Trans. Audio Electroacoustics AU-21(6), 506–526, 1973. doi:10.1109/TAU.1973.1162525.
- O. Herrmann, L. R. Rabiner and D. S. K. Chan, "Practical design rules for optimum finite impulse response lowpass digital filters," Bell System Technical Journal 52(6), 769–799, 1973.
- E. B. Hogenauer, "An economical class of digital filters for decimation and interpolation," IEEE Trans. Acoustics, Speech, and Signal Processing 29(2), 155–162, 1981. doi:10.1109/TASSP.1981.1163535. The CIC filter.
- R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall, 1983. Polyphase structures and multistage decimation.
- S. J. Orfanidis, "Lecture notes on elliptic filter design," Rutgers University, 2006. The Landen-transformation algorithm used in Demo 09.3.
- S. Bjaerum, H. Torp and K. Kristoffersen, "Clutter filter design for ultrasound color flow imaging," IEEE Trans. Ultrasonics, Ferroelectrics, and Frequency Control 49(2), 204–216, 2002. FIR, IIR with initialization, regression and their power responses on short ensembles.
- C. Demené et al., "Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity," IEEE Trans. Medical Imaging 34(11), 2271–2285, 2015. doi:10.1109/TMI.2015.2428634. The SVD filter that replaces frequency-based wall filters for long ensembles (Module 06).
- R. G. Lyons, Understanding Digital Signal Processing, 3rd ed., Prentice Hall, 2010. Practical chapters on CIC filters, half-bands and fixed-point effects.