Part I · Module 03

Convolution, LTI systems and filters, from superposition

Every linear step in the ultrasound chain (transducer, analog front end, RF and IQ filters, the pulse-echo model of an RF line, wall filters, persistence) is a convolution with an impulse response. This module builds that from two assumptions, linearity and time invariance, then shows why the same operation is a plain multiplication in frequency, what the real and the imaginary parts of the frequency response do to a sinusoid, and why undoing a filter is harder than applying it.

You will be able to

Test a system for linearity and time invariance, and name the ultrasound blocks that fail each test
Derive the convolution sum from superposition and compute a convolution by hand, as a table and as polynomial multiplication
Prove that \(e^{j\omega n}\) is an eigenfunction of every LTI system and read \(H(e^{j\omega})\) as a gain and a rotation
Say what \(\mathrm{Re}\,H\) and \(\mathrm{Im}\,H\) do to a cosine, and read phase delay and group delay off \(\angle H\)
Derive the convolution theorem and check it numerically at one frequency
Run difference equations by hand, read \(\lvert H\rvert\) from pole and zero distances, and state the stability condition
Explain why band-pass filtering an RF line helps SNR and why inverse filtering amplifies noise

Assumes

Module 01 · Complex numbers and phasors: phasors, \(e^{j\theta}\), \(j\) as a 90° rotation
Module 02 · The DFT and FFT by hand: the DFT as correlation with rotating phasors, conjugate symmetry, leakage
Next: Module 04 (sampling, the ultrasound signal chain as LTI) and Module 09 (filter design) build on this one

Exact formulaclosed form or exact finite sum (convolution, DTFT, recursion) evaluated as writtenSimulationnumerical solver run here (not used on this page)Approx. modelsimplified physics, assumptions stated under the badge (Gaussian pulse, 1D pulse-echo convolution, white noise)Schematicdrawing for intuition (block diagrams, flip-and-slide, z-plane anatomy), not computed

Conventions in this module

Discrete time \(n\) (integer samples). \(\omega\) is digital frequency in rad/sample, \(-\pi < \omega \le \pi\); \(f = \omega f_s/(2\pi)\) in Hz. \(\delta[n]\) is the unit impulse (1 at \(n = 0\), else 0). \(x\) has length \(L\), \(h\) has length \(M\). \(H(e^{j\omega}) = \sum_n h[n]\, e^{-j\omega n}\) is the frequency response (the DTFT of \(h\)). Convolution is written \(\ast\). RF examples use \(f_s = 40\) MHz, \(f_0 = 5\) MHz, \(c = 1540\) m/s.

03.1Systems as boxes: linearity and time invariance

A system \(T\) takes a whole input sequence \(x\) and returns an output sequence \(y = T\{x\}\). Nothing else is assumed. Two properties, tested separately, turn that black box into something you can describe with one sequence.

\[ T\{a\,x_1 + b\,x_2\} = a\,T\{x_1\} + b\,T\{x_2\} \qquad \text{for every } a,\ b,\ x_1,\ x_2 \]

(03.1)

\(T\{\cdot\}\): the system: maps a whole input sequence to a whole output sequence
\(x_1, x_2\): any two input sequences
\(a, b\): any two (real or complex) constants, dimensionless
\(T\{\cdot\}\): as above: the system operator

Linearity = additivity and homogeneity together. A quick necessary test: zero in must give zero out (take \(a = b = 0\)).

\[ T\{x[n]\} = y[n] \;\Longrightarrow\; T\{x[n-k]\} = y[n-k] \qquad \text{for every shift } k \]

(03.2)

\(n\): sample index (integer)
\(k\): delay applied to the input, samples
\(x[n],\ y[n]\): input and output sequences
\(T\{\cdot\}\): the system operator, as in (03.1)
\(\Longrightarrow\): implies

Time invariance: delaying the input only delays the output. The box does not know what time it is.

System	Rule	Linear?	Time-invariant?	The counterexample, by hand
3-point average	\(y[n] = (x[n] + x[n-1] + x[n-2])/3\)	yes	yes	None: this is an LTI filter
Delay	\(y[n] = x[n-2]\)	yes	yes	None
Squarer	\(y[n] = x[n]^2\)	no	yes	\(x = \delta\) gives \(\delta\); \(x = 2\delta\) gives \(4\delta\), not \(2\delta\)
Offset	\(y[n] = x[n] + 0.5\)	no	yes	\(x = 0\) gives 0.5, not 0
TGC (time-varying gain)	\(y[n] = g[n]\,x[n]\), \(g[n] = 1 + n/8\)	yes	no	\(\delta[n-2] \to 1.25\,\delta[n-2]\), but \(\delta[n-10] \to 2.25\,\delta[n-10]\), not \(1.25\,\delta[n-10]\)
Time reversal	\(y[n] = x[-n]\)	yes	no	Delaying the input advances the output
Envelope, log compression	\(\lvert x_a[n]\rvert\), \(20\log_{10}\lvert\cdot\rvert\)	no	yes	Envelope of \(-x\) equals envelope of \(x\)

Schematic Figure 03.A · The receive chain as a row of boxes, each tagged by which of the two tests it passes (block sizes not meaningful)

Demo 03.1

System tester: apply the two definitions literally

Exact formulaEach system rule evaluated sample by sample on two short test sequences x₁, x₂; no approximation.

input to the box\(T\{\)combined input\(\}\)combination of separate outputs (open circles)

System

Test

a 2.0

b −1.0

Shift k 6

Try: TGC passes the linearity test for every \(a, b\) and fails time invariance for every \(k > 0\); the squarer does the opposite. Set the offset system to \(a = 0.5\), \(b = 0.5\): the linearity test passes by accident because \(a + b = 1\). One passing test proves nothing; one failing test is a proof.

Ultrasound link

From the element to the beamformer output the receive chain is close to linear until something clips (quantization adds noise but is not a gain). TGC makes it time-varying, but the gain is known and can be divided out, which is why QUS methods either work on pre-TGC data or remove the TGC curve first. Nonlinear propagation at high mechanical index breaks linearity in the medium (the second harmonic that harmonic imaging uses), and envelope detection and log compression are nonlinear on purpose. Rule: filter, beamform and deconvolve before the envelope, never after.

03.2The impulse response, and the convolution sum from superposition

Four small steps take you from the two definitions to the formula that describes every LTI system.

Any sequence is a stack of scaled, shifted impulses, (03.3). For \(x = [1, 2, 3]\): \(x = 1\cdot\delta[n] + 2\cdot\delta[n-1] + 3\cdot\delta[n-2]\).
Name the response to one impulse. \(h[n] = T\{\delta[n]\}\), the impulse response.
Time invariance moves it. \(T\{\delta[n-k]\} = h[n-k]\).
Linearity adds the pieces. \(T\{\sum_k x[k]\,\delta[n-k]\}\) \(= \sum_k x[k]\,T\{\delta[n-k]\}\) \(= \sum_k x[k]\,h[n-k]\).

\[ x[n] = \sum_{k} x[k]\,\delta[n-k], \qquad h[n] = T\{\delta[n]\} \]

(03.3)

\(\delta[n]\): unit impulse: 1 at \(n = 0\), 0 elsewhere
\(x[k]\): the value of the input at sample \(k\) (a number, the weight of the impulse at \(k\))
\(h[n]\): impulse response: the output when the input is \(\delta[n]\)
\(x[n]\): input sequence
\(n,\ k\): sample indices
\(\sum_k\): sum over all \(k\)
\(T\{\cdot\}\): the system operator, as in (03.1)

\[ y[n] = (x \ast h)[n] = \sum_{k} x[k]\, h[n-k] \]

(03.4)

\(y[n]\): output sample at index \(n\)
\(x[k]\): input sample at index \(k\), length \(L\)
\(h[n-k]\): impulse response delayed by \(k\) samples, length \(M\)
\(k\): summation index over all input samples
\(\ast\): convolution
\(x,\ h,\ y\): input, impulse response, output
\(n\): output sample index

Input-side reading: input sample \(x[k]\) launches a copy of \(h\), scaled by \(x[k]\), starting at \(n = k\). The output is the sum of all the copies. One sequence, \(h\), describes the system completely.

Demo 03.2

Superposition builder: every input sample launches a copy of h

Exact formulaSum (03.4) accumulated copy by copy; the 9-tap RF pulse is a sampled Gaussian-windowed cosine used as a test sequence.

\(x[k]\) and its copy \(x[k]\,h[n-k]\) share a colorrunning sumfinal \(y\) (open)

Input x

Impulse response h

Copies added 3

Try: step the by-hand pair one copy at a time and compare with the table below. With "Four reflectors" and the RF pulse, the two reflectors one sample apart overlap and their copies interfere: this is the 1D pulse-echo model, RF line \(= \sum_k \rho_k\, p[n - n_k]\).

By hand. \(x = [1, 2, 3]\), \(h = [2, 1, -1]\). Each row is one launched copy; the output is the column sum.

Copy (columns: n)	0	1	2	3	4
\(x[0]\,h[n] = 1\cdot[2, 1, -1]\)	2	1	−1
\(x[1]\,h[n-1] = 2\cdot[2, 1, -1]\)		4	2	−2
\(x[2]\,h[n-2] = 3\cdot[2, 1, -1]\)			6	3	−3
\(y[n]\) = column sum	2	5	7	1	−3

\[ \left(1 + 2z^{-1} + 3z^{-2}\right)\left(2 + z^{-1} - z^{-2}\right) = 2 + 5z^{-1} + 7z^{-2} + z^{-3} - 3z^{-4} \]

(03.5)

\(z^{-1}\): a one-sample delay marker: the coefficient of \(z^{-n}\) is the sample at index \(n\)
\(z\): complex variable of the z-transform; \(z^{-k}\) marks a delay of \(k\) samples
\(\text{coefficients}\): sample values of the sequences \((1,2,3)\) and \((2,1,-1)\); the product is their convolution

Convolution is polynomial multiplication: the powers of \(z^{-1}\) keep track of the shifts. It is also long multiplication without carries: conv([1 2 3], [2 1 1]) = [2 5 9 5 3] and 123 × 211 = 25 953.

\[ L_y = L + M - 1, \qquad \sum_n y[n] = \Big(\sum_n x[n]\Big)\Big(\sum_n h[n]\Big) \]

(03.6)

\(L, M\): lengths of \(x\) and \(h\), samples
\(L_y\): length of \(y = x \ast h\), samples
\(x[n],\ h[n],\ y[n]\): input, impulse response, output
\(n\): sample index
\(\sum_n\): sum over all samples

Two checks you can do in your head.

Length: \(L + M - 1 = 3 + 3 - 1 = 5\) samples, because the last copy starts at \(n = L - 1\) and lasts \(M\) samples. Sum: \(\sum y = 12 = \sum x \cdot \sum h = 6 \cdot 2\). Setting \(z = 1\) in (03.5) proves it in general; in 03.6 you will recognize it as \(Y = XH\) at \(\omega = 0\).

Scatter and gather, as in DAS

The table can be filled by rows (each input scatters a copy of \(h\)) or by columns (each output gathers the inputs, weighted by \(h\) read backwards). Both give the same numbers. It is the same choice as writing a beamformer element-driven (each channel sample is spread to every pixel it affects) or pixel-driven (each pixel gathers along its delay curve). Section 03.3 is the gather view.

03.3Flip, slide, multiply, sum: the output-side view, and the properties

Fix the output index \(n\) and read \(y[n] = \sum_k x[k]\,h[n-k]\) as a function of \(k\). The sequence \(h[n-k]\) is \(h\) reversed (\(k \to -k\)) and moved right by \(n\). So one output sample is an inner product between \(x\) and a flipped, shifted copy of \(h\).

\(y[2]\) step by step. Flip \(h = [2, 1, -1]\) to get \([-1, 1, 2]\) at \(k = -2, -1, 0\). Shift right by 2: \(h[2-k] = (-1, 1, 2)\) at \(k = 0, 1, 2\). Multiply with \(x = (1, 2, 3)\) and add, (03.7): the middle column of the table.

\[ \begin{aligned} y[2] &= \sum_k x[k]\,h[2-k] = x[0]\,h[2] + x[1]\,h[1] + x[2]\,h[0] \\ &= 1\cdot(-1) + 2\cdot 1 + 3\cdot 2 = 7 \end{aligned} \]

(03.7)

\(x[k],\ h[n-k]\): input \(x=(1,2,3)\), impulse response \(h=(2,1,-1)\), as in (03.5)
\(y[2]\): third output sample
\(k\): summation index
\(\sum_k\): sum over all \(k\)

Gather view: the indices inside each product always add up to the output index, \(k + (n-k) = n\).

Schematic Figure 03.B · Flip, slide, multiply, sum for \(y[2]\) of the by-hand pair (stem heights to scale)

Demo 03.3

The same convolution twice: flip-and-slide in time, multiply in frequency

Exact formulaDirect convolution sum, and the DTFT sums evaluated on 257 frequencies in [0, π]; the dots are the DFT of the computed y. Noise in the "noisy step" pair is seeded.

Time domain: \(y[n] = \sum_k x[k]\,h[n-k]\)

Frequency domain: \(Y(e^{j\omega}) = X(e^{j\omega})\,H(e^{j\omega})\)

\(x\)\(h\) (flipped, shifted on the left)product and \(y\)dots: DFT of the computed \(y\)

Pair

Frequency view

Output index n 2

Probe frequency ω/π 0.50

Try: in the by-hand pair at \(\omega/\pi = 0.5\), the readout multiplies \(X = -2 - 2j\) by \(H = 3 - j\) and gets \(-8 - 4j\), which equals \(Y\) computed straight from \(y = [2, 5, 7, 1, -3]\). Switch to "Real and imaginary": \(\mathrm{Re}\,Y\) is not \(\mathrm{Re}\,X \cdot \mathrm{Re}\,H\). In "Rect 6 ∗ rect 4" the flat top of the trapezoid is the range of \(n\) where the short rect is entirely inside the long one, and the zeros of \(\lvert Y\rvert\) are the zeros of \(\lvert X\rvert\) and \(\lvert H\rvert\) together.

Properties, each a one-line substitution in the sum:

Property	Statement	Proof idea	Ultrasound meaning
Commutative	\(x \ast h = h \ast x\)	substitute \(m = n - k\)	"Signal" and "system" are interchangeable roles
Associative	\((x \ast h_1) \ast h_2 = x \ast (h_1 \ast h_2)\)	swap the order of two sums	Pulse-echo response = excitation ∗ transmit response ∗ receive response; a cascade is one filter
Distributive	\(x \ast (h_1 + h_2) = x \ast h_1 + x \ast h_2\)	linearity of the sum	Parallel paths add: direct echo plus a reverberation path
Identity, delay	\(x \ast \delta = x\), \(\;x \ast \delta[n-d] = x[n-d]\)	only one term survives	A delay is a convolution: DAS delays are filters
Length	\(L + M - 1\)	first copy starts at 0, last ends at \(L + M - 2\)	Filtering lengthens a record by \(M - 1\) samples

\[ y_{\mathrm{circ}}[n] = \sum_{r=-\infty}^{\infty} y[n + rN], \qquad 0 \le n \le N-1 \]

(03.8)

\(N\): DFT (FFT) length, samples
\(y\): the linear convolution \(x \ast h\), length \(L + M - 1\)
\(r\): integer: which period the tail wraps from
\(y_{\mathrm{circ}}[n]\): circular (wrap-around) convolution output
\(y[n+rN]\): linear convolution result at shifted index
\(n\): sample index
\(\sum_r\): sum over all \(r\)

Inverse DFT of the product of two N-point DFTs: the linear result folded onto one period. No wrap if \(N \ge L + M - 1\).

Pitfall: FFT convolution is circular

Multiplying two N-point DFTs and inverting gives circular convolution (03.8): output samples past \(N - 1\) wrap to the start. For the by-hand pair with \(N = 4\), the −3 at \(n = 4\) lands on \(n = 0\) and the result is \([-1, 5, 7, 1]\). Pad both sequences to \(N \ge L + M - 1 = 5\) and the wrap disappears (Module 04, Demo 04.7). For long RF lines, overlap-add applies the same rule block by block.

03.4Why complex exponentials are eigenfunctions of every LTI system

Feed one rotating phasor \(x[n] = e^{j\omega n}\) into the convolution, using the commuted form \(y[n] = \sum_k h[k]\,x[n-k]\):

\[ \begin{aligned} y[n] &= \sum_k h[k]\, e^{j\omega (n-k)} = e^{j\omega n} \sum_k h[k]\, e^{-j\omega k} = H(e^{j\omega})\, e^{j\omega n} \\ H(e^{j\omega}) &= \sum_k h[k]\, e^{-j\omega k} \end{aligned} \]

(03.9)

\(e^{j\omega n}\): input phasor: unit length, turning \(\omega\) rad per sample
\(\omega\): digital frequency, rad/sample (\(\omega = 2\pi f/f_s\))
\(h[k]\): impulse response taps
\(H(e^{j\omega})\): frequency response: one complex number per \(\omega\), dimensionless gain
\(y[n]\): output for the complex exponential input
\(n,\ k\): sample indices
\(h[k]\): impulse response taps
\(j\): imaginary unit, \(j^2=-1\)

The factor \(e^{j\omega n}\) does not depend on \(k\), so it leaves the sum. What remains is one complex number per frequency: the frequency response, the correlation of \(h\) with a probe phasor (Module 02).

"Eigenfunction" means the same thing as eigenvector of a matrix: the output is the input times a number. For a convolution that number is \(H(e^{j\omega})\), and the only things it can do to the phasor are scale it by \(\lvert H\rvert\) and rotate it by \(\angle H\). Consequence: an LTI system never creates a frequency that was not in the input. New frequencies need a nonlinearity (the squarer turns \(\cos\omega n\) into \(\tfrac12 + \tfrac12\cos 2\omega n\)) or time variation (TGC multiplies, which convolves spectra).

A real cosine is two counter-rotating phasors, \(\cos\omega n = (e^{j\omega n} + e^{-j\omega n})/2\). Each passes through scaled by its own eigenvalue. For real \(h\), \(H(e^{-j\omega}) = H(e^{j\omega})^{*}\), the two outputs are conjugates and their sum is real:

\[ \cos\omega n \;\longrightarrow\; \tfrac12\left[H e^{j\omega n} + H^{*} e^{-j\omega n}\right] = \mathrm{Re}\{H e^{j\omega n}\} = \lvert H\rvert \cos\!\left(\omega n + \angle H\right) \]

(03.10)

\(H^{*}\): complex conjugate of \(H(e^{j\omega})\), equal to \(H(e^{-j\omega})\) for real \(h\)
\(\lvert H\rvert,\ \angle H\): gain (dimensionless) and phase (rad) at \(\omega\)
\(H\): shorthand for \(H(e^{j\omega})\)
\(n\): sample index
\(\omega\): digital frequency, rad/sample
\(\mathrm{Re}\{\cdot\}\): real part
\(j\): imaginary unit, \(j^2=-1\)
\(\longrightarrow\): passes through the system to give

Same frequency, amplitude multiplied by \(\lvert H(e^{j\omega})\rvert\), phase advanced by \(\angle H(e^{j\omega})\) (a negative angle is a delay).

By hand. Two-point average \(h = [\tfrac12, \tfrac12]\) at \(\omega = \pi/2\) (a quarter of the sampling rate), (03.11). The input \(\cos(\pi n/2) = 1, 0, -1, 0, 1\) for \(n = 0\ldots4\) gives \(y[n] = (x[n] + x[n-1])/2\) = 0.5, 0.5, −0.5, −0.5, 0.5, and indeed \(0.707\cos(\pi n/2 - \pi/4)\) = 0.5, 0.5, −0.5, −0.5, 0.5. The −45° at 90° per sample is a delay of half a sample, the center of the two taps.

\[ \begin{gathered} H(e^{j\pi/2}) = \tfrac12 + \tfrac12 e^{-j\pi/2} = \tfrac12 - \tfrac12 j \\ \lvert H\rvert = \frac{1}{\sqrt 2} = 0.707, \qquad \angle H = -45^\circ \end{gathered} \]

(03.11)

\(H(e^{j\omega})\): frequency response of the two-tap average \(h=(\tfrac12,\tfrac12)\) evaluated at \(\omega=\pi/2\)
\(\pi/2\): quarter of the sampling rate, rad/sample
\(\angle H\): phase of \(H\), degrees
\(j\): imaginary unit, \(j^2=-1\)

Demo 03.4

Eigenfunction: H is a sum of rotated taps, and it scales and rotates the input

Exact formulaTap arrows and resultant from (03.9); output from the convolution sum. The one-pole filter is truncated at 60 taps (0.8⁶⁰ ≈ 1.5 × 10⁻⁶).

Complex plane: \(H = \sum_k h[k]\, e^{-j\omega k}\), tap arrows tip to tail

Time domain: \(\cos\omega n\) in, \(\lvert H\rvert\cos(\omega n + \angle H)\) out

tap arrows \(h[k]\,e^{-j\omega k}\)\(H\) (resultant)input \(\cos\omega n\)output \(y[n]\)\(\mathrm{Re}\,H\cos\omega n\) (dashed)\(-\mathrm{Im}\,H\sin\omega n\) (dotted)

Filter h

ω/π 0.50

Show output = Re H·cos ωn − Im H·sin ωn

Try: the one-pole filter draws a spiral of shrinking arrows; its resultant is \(H\). The centered smoother always gives a horizontal resultant (real \(H\), zero phase); the centered difference always gives a vertical one (imaginary \(H\), ±90°). For "Delay by 3" the single arrow just turns: \(\lvert H\rvert = 1\) and \(\angle H = -3\omega\). Near \(\omega/\pi = 0.667\) the 3-point average has its three arrows closing into a triangle: \(\lvert H\rvert \to 0\), the null at exactly \(2\pi/3\).

Ultrasound link: a delay is a phase ramp

A delay \(\tau\) is the LTI system \(h(t) = \delta(t - \tau)\) with \(H(f) = e^{-j2\pi f\tau}\): unit gain, phase \(-2\pi f\tau\). Applying a beamforming delay in the frequency domain is exactly this multiplication, and the narrowband shortcut of rotating IQ data by \(-2\pi f_0\tau\) keeps only the value at \(f_0\). Aberration correction estimates the same per-channel phase ramps.

03.5What the imaginary part of H means

Write \(H = \mathrm{Re}\,H + j\,\mathrm{Im}\,H\) and expand \(e^{-j\omega k} = \cos\omega k - j\sin\omega k\):

\[ \mathrm{Re}\,H(e^{j\omega}) = \sum_k h[k]\cos\omega k, \qquad \mathrm{Im}\,H(e^{j\omega}) = -\sum_k h[k]\sin\omega k \]

(03.12)

\(H(e^{j\omega})\): frequency response, as in (03.9)
\(h[k]\): impulse response (real) taps
\(k\): tap index
\(\omega\): digital frequency, rad/sample
\(\mathrm{Re},\ \mathrm{Im}\): real and imaginary parts

\(\mathrm{Re}\,H\) is the correlation of \(h\) with a cosine probe, \(\mathrm{Im}\,H\) minus its correlation with a sine probe (Module 02). For real \(h\), \(\mathrm{Re}\,H\) is even in \(\omega\) and \(\mathrm{Im}\,H\) is odd.

Now push a cosine through and keep both parts. Multiply \((\mathrm{Re}\,H + j\,\mathrm{Im}\,H)(\cos\omega n + j\sin\omega n)\) and take the real part:

\[ \cos\omega n \;\longrightarrow\; y[n] = \mathrm{Re}\,H \cdot \cos\omega n \;-\; \mathrm{Im}\,H \cdot \sin\omega n \]

(03.13)

\(\mathrm{Re}\,H\): in-phase gain at \(\omega\), dimensionless
\(\mathrm{Im}\,H\): quadrature gain at \(\omega\), dimensionless (its negative multiplies the late copy)
\(\sin\omega n\): \(= \cos(\omega n - \pi/2)\): the input a quarter period late
\(H\): shorthand for \(H(e^{j\omega})\), as in (03.10)
\(n\): sample index
\(\omega\): digital frequency, rad/sample
\(y[n]\): output sequence
\(\longrightarrow\): passes through the system to give

The output is a real signal made of two real pieces: an in-phase copy of the input with gain \(\mathrm{Re}\,H\), and a quadrature copy, a quarter period late, with gain \(-\mathrm{Im}\,H\).

\(\mathrm{Re}\,H\) is the in-phase gain: how much of the output stays in step with the input.
\(\mathrm{Im}\,H\) is the quadrature gain: how much of the output comes out shifted by 90°. \(\mathrm{Im}\,H < 0\) adds a copy a quarter period late (the filter delays); \(\mathrm{Im}\,H > 0\) adds one a quarter period early (a differentiator leads). Nothing in the output is imaginary; \(\mathrm{Im}\,H\) scales a real sine.
In phasor language: \(H e^{j\omega n} = \mathrm{Re}\,H \cdot e^{j\omega n} + \mathrm{Im}\,H \cdot (j\,e^{j\omega n})\), and \(j\,e^{j\omega n}\) is the input arrow turned 90° ahead (Module 01: \(j\) is the quarter-turn operator).
Polar form packs the two gains into one, (03.14).

\[ \lvert H\rvert = \sqrt{(\mathrm{Re}\,H)^2 + (\mathrm{Im}\,H)^2}, \qquad \angle H = \operatorname{atan2}\!\left(\mathrm{Im}\,H,\ \mathrm{Re}\,H\right) \]

(03.14)

\(H\): frequency response \(H(e^{j\omega})\)
\(\mathrm{Re}\,H,\ \mathrm{Im}\,H\): its real and imaginary parts
\(\lvert H\rvert,\ \angle H\): gain (dimensionless) and phase (rad)
\(\operatorname{atan2}\): four-quadrant arctangent

By hand. Two-point average at \(\omega = \pi/2\): \(\mathrm{Re}\,H = 0.5\), \(\mathrm{Im}\,H = -0.5\). The output is \(0.5\cos(\pi n/2) + 0.5\sin(\pi n/2)\). At \(n = 1\): \(0 + 0.5 = 0.5\); at \(n = 3\): \(0 - 0.5 = -0.5\). Same numbers as in 03.4, now explained as half an in-phase copy plus half a late quadrature copy.

\[ \begin{aligned} h_e[n] &= \tfrac12\left(h[n] + h[-n]\right) \;\longleftrightarrow\; \mathrm{Re}\,H \\ h_o[n] &= \tfrac12\left(h[n] - h[-n]\right) \;\longleftrightarrow\; j\,\mathrm{Im}\,H \end{aligned} \]

(03.15)

\(h_e, h_o\): even and odd parts of a real \(h\) about \(n = 0\); \(h = h_e + h_o\)
\(h[n]\): real impulse response
\(h[-n]\): time-reversed impulse response
\(\mathrm{Re}\,H,\ \mathrm{Im}\,H\): real and imaginary parts of \(H(e^{j\omega})\)
\(\longleftrightarrow\): is the transform pair of
\(j\): imaginary unit, \(j^2=-1\)

The even part of \(h\) (about \(n = 0\)) makes the real part of \(H\), the odd part makes the imaginary part. A causal \(h\) is never even unless it is a single tap at \(n = 0\), so causal filters always have some phase.

Phase delay and group delay. Two delays can be read off the phase \(\varphi(\omega) = \angle H(e^{j\omega})\):

\[ \tau_p(\omega) = -\frac{\varphi(\omega)}{\omega}, \qquad \tau_g(\omega) = -\frac{d\varphi}{d\omega} \]

(03.16)

\(\varphi(\omega)\): unwrapped phase of \(H(e^{j\omega})\), rad
\(\tau_p\): phase delay: delay of the carrier at \(\omega\), samples (divide by \(f_s\) for seconds)
\(\tau_g\): group delay: delay of the envelope of a narrowband pulse centered at \(\omega\), samples
\(\omega\): digital frequency, rad/sample
\(d\varphi/d\omega\): slope of the phase

For a narrowband pulse at \(\omega_0\), the carrier is delayed by \(\tau_p(\omega_0)\) and the envelope by \(\tau_g(\omega_0)\). Linear phase \(\varphi = -\omega d\) gives \(\tau_p = \tau_g = d\) at every \(\omega\), so the pulse shape is preserved.

Worked numbers. Two-point average: \(\varphi = -\omega/2\), so \(\tau_p = \tau_g = 0.5\) sample at all \(\omega\). A symmetric FIR with \(M\) taps: \(\tau = (M-1)/2\); a 41-tap band-pass at 40 MHz delays by 20 samples = 0.5 µs, which is 0.385 mm of depth at \(c/2\). One-pole smoother \(y[n] = a\,y[n-1] + (1-a)\,x[n]\), (03.17): 4.00 samples at DC for \(a = 0.8\) and −0.44 at \(\omega = \pi\). Group delay of a non-linear-phase filter varies with frequency and distorts the pulse.

\[ \tau_g(\omega) = \frac{a\cos\omega - a^2}{1 - 2a\cos\omega + a^2}, \qquad \tau_g(0) = \frac{a}{1-a}, \qquad \tau_g(\pi) = -\frac{a}{1+a} \]

(03.17)

\(a\): feedback coefficient (pole position) of the one-pole smoother, \(0 < a < 1\)
\(\tau_g(\omega)\): group delay of the one-pole filter, samples
\(\omega\): digital frequency, rad/sample
\(\tau_g(0),\ \tau_g(\pi)\): group delay at DC and at Nyquist, samples

Intuition: the same split as impedance

An impedance \(Z = R + jX\) has the same anatomy. The resistive part \(R\) gives current in phase with voltage, the reactive part \(X\) gives current 90° out of phase. \(\mathrm{Re}\,H\) and \(\mathrm{Im}\,H\) are the "resistive" and "reactive" parts of a filter's response at one frequency.

Ultrasound link

The Hilbert transformer used to build IQ data is the filter with \(H = -j\,\mathrm{sgn}(\omega)\): \(\mathrm{Re}\,H = 0\) everywhere, so its output is all quadrature, cos becomes sin (Module 08). The other place \(\mathrm{Im}\,H\) bites is channel matching: if two receive channels have filters whose phases differ by 10° at 5 MHz, that is a phase-delay mismatch of \(10/360 \times 200\) ns = 5.6 ns, which the beamformer cannot tell from an aberration.

03.6The convolution theorem, step by step

The eigenfunction property already contains the theorem. Three steps:

Write \(x\) as a sum of phasors (inverse DTFT), (03.18). For an N-point DFT it is a sum of N phasors (Module 02).
Each phasor is multiplied by \(H\) at its own frequency (03.4): \(X(e^{j\omega})\,e^{j\omega n} \to H(e^{j\omega})\,X(e^{j\omega})\,e^{j\omega n}\).
Linearity adds the outputs: \(y[n] = \frac{1}{2\pi}\int H X\, e^{j\omega n}\, d\omega\), which says \(Y = H X\).

\[ x[n] = \frac{1}{2\pi}\int_{-\pi}^{\pi} X(e^{j\omega})\, e^{j\omega n}\, d\omega, \qquad X(e^{j\omega}) = \sum_n x[n]\, e^{-j\omega n} \]

(03.18)

\(X(e^{j\omega})\): DTFT of \(x\): complex amplitude of the phasor at \(\omega\)
\(x[n]\): sequence
\(n\): sample index
\(\omega\): digital frequency, rad/sample
\(\int_{-\pi}^{\pi}\): integral over one period of the spectrum
\(j\): imaginary unit, \(j^2=-1\)

The direct algebra is one substitution, \(m = n - k\):

\[ \begin{aligned} Y(e^{j\omega}) &= \sum_n \sum_k x[k]\,h[n-k]\, e^{-j\omega n} \\ &= \sum_k x[k]\, e^{-j\omega k} \cdot \sum_m h[m]\, e^{-j\omega m} = X(e^{j\omega})\, H(e^{j\omega}) \end{aligned} \]

(03.19)

\(X, H, Y\): DTFTs of input, impulse response and output, complex
\(m\): \(= n - k\), the index along \(h\)
\(x[k],\ h[n-k]\): input and shifted impulse response
\(k,\ n\): sample indices
\(\sum\): sum over all indices shown
\(j\): imaginary unit, \(j^2=-1\)

Split \(e^{-j\omega n} = e^{-j\omega k}\, e^{-j\omega(n-k)}\). The double sum factors because each exponential depends on only one index. Convolution in time is multiplication in frequency.

By hand at \(\omega = \pi/2\), where \(e^{-j\pi n/2} = 1, -j, -1, j, 1\) for \(n = 0\ldots4\):

Quantity	Sum	Rectangular	Polar
\(X = \sum x[n]\, e^{-j\pi n/2}\)	\(1 - 2j - 3\)	\(-2 - 2j\)	2.83 ∠ −135°
\(H = \sum h[n]\, e^{-j\pi n/2}\)	\(2 - j + 1\)	\(3 - j\)	3.16 ∠ −18.4°
\(X H\)	\((-2)(3) - (-2)(-1) + j[(-2)(-1) + (-2)(3)]\)	\(-8 - 4j\)	8.94 ∠ −153.4°
\(Y = \sum y[n]\, e^{-j\pi n/2}\)	\(2 - 5j - 7 + j - 3\)	\(-8 - 4j\)	8.94 ∠ −153.4°

In polar form, magnitudes multiply (2.83 × 3.16 = 8.94) and angles add (−135° − 18.4° = −153.4°). In rectangular form, the parts mix: \(\mathrm{Re}\,Y = \mathrm{Re}\,X\,\mathrm{Re}\,H - \mathrm{Im}\,X\,\mathrm{Im}\,H\). Filtering in frequency is complex multiplication, one complex number per bin.

Pitfall

Multiplying only the real parts, or only the magnitudes, of two spectra is not a convolution. Magnitude-only multiplication gives a zero-phase result that is centered on \(n = 0\) and is generally a different signal. Keep the complex values, and keep the padding.

\[ \begin{aligned} \text{DTFT:}\quad x[n]\, w[n] &\;\longleftrightarrow\; \frac{1}{2\pi}\int_{-\pi}^{\pi} X(e^{j\theta})\, W(e^{j(\omega-\theta)})\, d\theta \\ \text{DFT:}\quad x[n]\, w[n] &\;\longleftrightarrow\; \frac{1}{N}\sum_{m=0}^{N-1} X[m]\, W[(k-m) \bmod N] \end{aligned} \]

(03.20)

\(w[n]\): a window or gate multiplying the signal sample by sample
\(\theta\): integration variable, rad/sample
\(X[m], W[m]\): N-point DFTs
\(x[n]\): signal
\(n\): sample index
\(X(e^{j\omega})\): DTFT of \(x\), as in (03.18)
\(W(e^{j\omega})\): DTFT of \(w\)
\(\omega\): digital frequency, rad/sample
\(N\): DFT length, samples
\(k,\ m\): DFT bin indices
\(\bmod N\): remainder after division by \(N\) (circular index)
\(\longleftrightarrow\): is the transform pair of
\(j\): imaginary unit, \(j^2=-1\)

Duality, (03.20): multiplication in time is (circular) convolution in frequency, divided by \(2\pi\) for the DTFT or by \(N\) for the DFT. Gating a tone with a window smears its spectral line into the window's transform: that is leakage (Module 02) and apodization (Module 04).

03.7Difference equations, H(z), poles and zeros, causality and stability

Most filters are implemented as a recursion with a few coefficients rather than a long \(h\):

\[ \begin{aligned} y[n] &= \sum_{k=0}^{Q} b_k\, x[n-k] \;-\; \sum_{k=1}^{P} a_k\, y[n-k] \\ \Longleftrightarrow\quad H(z) &= \frac{b_0 + b_1 z^{-1} + \cdots + b_Q z^{-Q}}{1 + a_1 z^{-1} + \cdots + a_P z^{-P}} = \frac{B(z)}{A(z)} \end{aligned} \]

(03.21)

\(b_k\): feedforward coefficients, \(k = 0\ldots Q\)
\(a_k\): feedback coefficients, \(k = 1\ldots P\) (\(a_0 = 1\))
\(z^{-1}\): one-sample delay operator; \(z\) is a complex variable
\(B(z), A(z)\): numerator and denominator polynomials in \(z^{-1}\)
\(y[n],\ x[n]\): output and input sequences
\(n\): sample index
\(Q,\ P\): number of feedforward and feedback taps (orders)
\(H(z)\): transfer function
\(z\): complex variable of the z-transform
\(\Longleftrightarrow\): is equivalent to

Without feedback (all \(a_k = 0\)) the impulse response is the \(b_k\) themselves: FIR. With feedback it never ends: IIR. On the unit circle \(z = e^{j\omega}\), \(H(z)\) is the frequency response.

Schematic Figure 03.C · A first-order difference equation as a block diagram (direct form I): delays, gains and one adder. The one-pole smoother is b₀ = 1 − a, b₁ = 0, a₁ = −a.

Filter	Recursion	Impulse response	Frequency response	What it does
3-point moving average	\(y[n] = (x[n] + x[n-1] + x[n-2])/3\)	⅓, ⅓, ⅓	\(\tfrac13 e^{-j\omega}(1 + 2\cos\omega)\)	Low-pass, null at \(\omega = 2\pi/3\) (\(f_s/3\)), delay 1 sample
First difference	\(y[n] = x[n] - x[n-1]\)	1, −1	\(1 - e^{-j\omega} = 2j\sin(\omega/2)\, e^{-j\omega/2}\)	High-pass, \(\lvert H\rvert = 2\lvert\sin(\omega/2)\rvert \approx \omega\): a derivative (× \(f_s\)), phase \(+90^\circ - \omega/2\)
One-pole smoother	\(y[n] = a\,y[n-1] + (1-a)\,x[n]\)	\((1-a)\,a^n\): for \(a = \tfrac12\), 0.5, 0.25, 0.125, …	\(\dfrac{1-a}{1 - a e^{-j\omega}}\)	Low-pass, DC gain 1, \(\lvert H(\pi)\rvert = (1-a)/(1+a)\)

Step response of the one-pole smoother by hand (\(a = 0.8\), \(x[n] = 1\) for \(n \ge 0\), \(y[-1] = 0\)): \(y[0] = 0.2\), \(y[1] = 0.8\cdot 0.2 + 0.2 = 0.36\), \(y[2] = 0.488\), \(y[3] = 0.5904\), \(y[4] = 0.67232\). In closed form, (03.22), the output creeps toward 1 with time constant \(-1/\ln a = 4.48\) samples.

\[ y[n] = (1-a)\sum_{m=0}^{n} a^m = 1 - a^{n+1}, \qquad \tau = -\frac{1}{\ln a} \approx \frac{1}{1-a} \]

(03.22)

\(\tau\): time constant, samples: \(a^n\) falls to \(1/e\) after \(\tau\) samples; the approximation holds for \(a\) near 1
\(y[n]\): step response of the one-pole filter \(y[n]=a\,y[n-1]+(1-a)\,x[n]\)
\(a\): pole, \(0<a<1\)
\(n,\ m\): sample indices
\(\ln\): natural logarithm

Demo 03.5

Difference-equation stepper

Exact formulaRecursion (03.21) run sample by sample; |H| from (03.21) on the unit circle. The noisy step is seeded.

\(x[n]\)\(y[n]\) computed so far\(\lvert H(e^{j\omega})\rvert\)

Recursion

Input

a (one-pole) 0.80

Computed up to n 4

Try: step the one-pole filter from \(n = 0\) and check 0.2, 0.36, 0.488, 0.5904. Set \(a = 0.95\) and feed the noisy step: smoother output, slower rise; that is the noise-versus-lag trade of display persistence. Set \(a = 1.05\): the impulse response grows and the readout flags the system as unstable.

Poles and zeros as distances. Factor numerator and denominator: \(H(z) = g\,\prod_k(1 - z_k z^{-1}) / \prod_k(1 - p_k z^{-1})\). Multiply top and bottom by the same power of \(z\); on the unit circle each factor becomes a vector from a root to the point \(e^{j\omega}\):

\[ \begin{aligned} \lvert H(e^{j\omega})\rvert &= \lvert g\rvert\, \frac{\prod_k \lvert e^{j\omega} - z_k\rvert}{\prod_k \lvert e^{j\omega} - p_k\rvert} \\ \angle H(e^{j\omega}) &= \angle g + \sum_k \angle\!\left(e^{j\omega} - z_k\right) - \sum_k \angle\!\left(e^{j\omega} - p_k\right) + (P - Q)\,\omega \end{aligned} \]

(03.23)

\(z_k\): zeros: roots of \(B(z)\), \(Q\) of them
\(p_k\): poles: roots of \(A(z)\), \(P\) of them
\(g\): overall gain constant (\(= b_0\))
\(\lvert e^{j\omega} - z_k\rvert\): length of the vector from a root to the point at angle \(\omega\) on the unit circle
\(H(e^{j\omega})\): frequency response
\(\omega\): digital frequency, rad/sample
\(P,\ Q\): numbers of poles and zeros (orders of \(A\) and \(B\))
\(\prod_k\): product over all \(k\)
\(\angle\): phase (angle) of a complex number, rad
\(j\): imaginary unit, \(j^2=-1\)

Walk around the circle: near a zero the product on top is small (a dip or null); near a pole the bottom is small (a peak). For the one-pole smoother (pole at \(a\), zero at 0): \(\lvert H\rvert = (1-a)/\lvert e^{j\omega} - a\rvert\), which is 1 at \(\omega = 0\) (distance \(1 - a\)) and \((1-a)/(1+a)\) at \(\omega = \pi\) (distance \(1 + a\)).

Schematic Figure 03.D · Anatomy of the z-plane: where each frequency sits on the unit circle, and the distances that set |H| (geometry to scale, f_s = 40 MHz labels)

Causality: \(h[n] = 0\) for \(n < 0\), so the output never uses future inputs. Real-time filters must be causal; offline processing (a stored frame) need not be, which is how zero-phase forward-backward filtering is possible.

Stability (BIBO), (03.24): every bounded input gives a bounded output if and only if \(h\) is absolutely summable. For a causal rational \(H(z)\) the test is geometric: all poles strictly inside the unit circle. At \(a = 1.05\) the impulse response doubles every \(\ln 2/\ln 1.05 = 14.2\) samples.

\[ \begin{gathered} \sum_{n} \lvert h[n]\rvert < \infty \\ \text{one-pole: } \sum_{n=0}^{\infty} \left\lvert (1-a)\,a^n\right\rvert = \frac{\lvert 1-a\rvert}{1-\lvert a\rvert} \ \ (\lvert a\rvert < 1), \quad \text{divergent for } \lvert a\rvert \ge 1 \end{gathered} \]

(03.24)

\(h[n]\): impulse response
\(n\): sample index
\(a\): pole of the one-pole filter, real
\(\sum_n\): sum over all samples
\(\infty\): infinity

Demo 03.6

Pole-zero vectors give |H| and ∠H

Exact formula|H| and ∠H from the vector products (03.23) for a conjugate zero pair and pole pair; h[n] from the same poles and zeros.

z-plane: zeros ○, poles ×, point \(e^{j\omega}\) on the unit circle

\(\lvert H\rvert\) (dB), \(\angle H\) (rad) and \(h[n]\)

zero vectorspole vectors\(\lvert H\rvert\), \(h[n]\)\(\angle H\)

Zero radius 1.00

Zero angle θ_z/π 0.75

Pole radius 0.90

Pole angle θ_p/π 0.25

ω/π 0.25

Try: roots come in conjugate pairs so that \(h\) is real. With the pole pair at \(0.9\,e^{\pm j\pi/4}\), put \(\omega\) on the pole angle: the near pole vector is 0.10 long and \(\lvert H\rvert\) peaks. Push the pole radius to 1.00 and beyond: the peak becomes infinite and the impulse response stops decaying. Put a zero on the circle and watch \(\lvert H\rvert\) fall to −∞ dB as \(\omega\) crosses it.

Ultrasound link

Display persistence is a one-pole smoother across frames. For white frame-to-frame noise the output variance is \((1-a)/(1+a)\) of the input, so \(a = 0.8\) acts like averaging \((1+a)/(1-a) = 9\) frames, at the price of a time constant of 4.5 frames. A Doppler wall filter is a high-pass with zeros at or near \(z = 1\) (DC in slow time); its design on short ensembles is Module 09.

03.8Filtering an RF line step by step, and why deconvolution is hard

The 1D pulse-echo model of Module 04, in discrete time: \(r[n] = (p \ast \rho)[n] + w[n]\), with \(p\) the two-way pulse, \(\rho\) the reflectivity along the beam and \(w\) white noise. Apply a band-pass FIR \(h_{\mathrm{bp}}\) and use the properties of 03.3:

\[ y = h_{\mathrm{bp}} \ast r = \left(h_{\mathrm{bp}} \ast p\right) \ast \rho + h_{\mathrm{bp}} \ast w \quad\Longleftrightarrow\quad Y = H_{\mathrm{bp}} P\, S + H_{\mathrm{bp}} W \]

(03.25)

\(p, P\): two-way (pulse-echo) pulse and its spectrum
\(\rho, S\): reflectivity along the beam and its spectrum
\(w, W\): additive white noise and its spectrum
\(h_{\mathrm{bp}}, H_{\mathrm{bp}}\): band-pass FIR and its frequency response
\(y\): filtered RF line
\(r\): RF line, \(r=p\ast\rho+w\)
\(\ast\): convolution
\(Y,\ S,\ W,\ P,\ H_{\mathrm{bp}}\): DTFTs of \(y,\ \rho,\ w,\ p,\ h_{\mathrm{bp}}\)
\(\Longleftrightarrow\): is equivalent to

Associativity: the output is still a sum of shifted, scaled copies, now of the filtered pulse \(h_{\mathrm{bp}} \ast p\), one per reflector. Distributivity: the filter acts on signal and noise separately.

The pulse here is a 5 MHz Gaussian with 60 % −6 dB fractional bandwidth (\(\sigma_t = 5.0\) samples at 40 MHz); the filter is a 41-tap Hamming-windowed band-pass. White noise spreads evenly over 0–20 MHz, the pulse lives in roughly 3–7 MHz, so keeping only that band removes most of the noise power and little of the signal. Expected values for this pulse, from (03.26):

\[ \begin{gathered} G_N = \sum_n h_{\mathrm{bp}}[n]^2, \qquad G_S = \frac{\sum_n \left(h_{\mathrm{bp}} \ast s\right)[n]^2}{\sum_n s[n]^2} \\ \text{SNR gain} = 10\log_{10}\frac{G_S}{G_N}\ \mathrm{dB} \end{gathered} \]

(03.26)

\(G_N\): noise power gain of the filter for white noise (Parseval: mean of \(\lvert H_{\mathrm{bp}}\rvert^2\) over frequency)
\(s\): noise-free RF line \(p \ast \rho\)
\(G_S\): signal power gain, measured on the noise-free line
\(n\): sample index
\(\sum_n\): sum over all samples
\(\ast\): convolution
\(\mathrm{SNR\ gain}\): improvement of SNR by the filter, dB
\(\log_{10}\): base-10 logarithm
\(\mathrm{dB}\): decibel, power ratio \(10\log_{10}\)

Pass band	Signal power	Noise power	SNR gain	Cost
2–8 MHz	−0.05 dB	−5.8 dB	+5.8 dB	Least noise rejection
3–7 MHz	−0.6 dB	−7.9 dB	+7.3 dB	Pulse nearly unchanged
4–6 MHz	−3.6 dB	−12.1 dB	+8.5 dB	Longer pulse: axial resolution lost

Deconvolution tries to go back, (03.27). Wherever \(\lvert P\rvert\) is small the division multiplies the noise \(W/P\) by a large number. For this pulse \(\lvert P\rvert\) is −24.1 dB at 2 and 8 MHz and −67 dB at 10 MHz relative to its peak, so the naive inverse amplifies noise at 10 MHz by 67 dB more than at 5 MHz. The Wiener filter (03.28) caps the gain.

\[ \hat S = \frac{Y}{P} = S + \frac{W}{P} \]

(03.27)

\(\hat S\): estimate of the reflectivity spectrum
\(Y,\ P,\ S,\ W\): spectra of data, pulse, reflectivity, noise, as in (03.25)

Exact in the signal, catastrophic in the noise: the error term \(W/P\) has power \(\lvert W\rvert^2/\lvert P\rvert^2\) at each frequency.

\[ G(\omega) = \frac{P^{*}(\omega)}{\lvert P(\omega)\rvert^2 + \varepsilon}, \qquad \max_{\omega} \lvert G\rvert = \frac{1}{2\sqrt{\varepsilon}} \ \text{ at } \lvert P\rvert = \sqrt{\varepsilon} \]

(03.28)

\(P(\omega)\): pulse spectrum (complex), normalized so \(\lvert P\rvert_{\max} = 1\) in the tables below
\(P^{*}\): complex conjugate of \(P\): the matched filter
\(\varepsilon\): regularization, same units as \(\lvert P\rvert^2\); optimal value = noise-to-signal power ratio per frequency
\(G\): Wiener deconvolution filter; estimate \(\hat S = G\,Y\)
\(G(\omega)\): deconvolution filter response
\(\omega\): angular frequency, rad/sample
\(\lvert P\rvert^2\): pulse power spectrum
\(\max_\omega\): maximum over frequency

Where \(\lvert P\rvert^2 \gg \varepsilon\), \(G \approx 1/P\) (inverse filter). Where \(\lvert P\rvert^2 \ll \varepsilon\), \(G \approx P^{*}/\varepsilon\) (a matched filter, scaled down). With \(\varepsilon = 10^{-3}\lvert P\rvert_{\max}^2\), the gain never exceeds \(15.8/\lvert P\rvert_{\max}\) (+24 dB).

Lowering \(\varepsilon\) widens the band that is inverted, so the effective pulse \(P G\) gets shorter, and the noise gain grows. For this pulse (peak SNR change = effective pulse peak over output noise std, relative to the input, for white noise; gains measured on a 512-point frequency grid, which just misses the exact maximum \(1/(2\sqrt{\varepsilon})\) = +34.0 dB at \(10^{-4}\)):

\(\varepsilon / \lvert P\rvert^2_{\max}\)	Max gain re \(1/\lvert P\rvert_{\max}\)	Effective pulse −6 dB width	Peak SNR change
no deconvolution	—	0.294 µs (0.227 mm)	0 dB
\(10^{-1}\)	+4.0 dB	0.276 µs (0.212 mm)	+4.4 dB
\(10^{-2}\)	+14.0 dB	0.214 µs (0.165 mm)	−0.6 dB
\(10^{-3}\)	+24.0 dB	0.178 µs (0.137 mm)	−7.6 dB
\(10^{-4}\)	+33.9 dB	0.155 µs (0.119 mm)	−15.7 dB
\(10^{-5}\)	+44.0 dB	0.138 µs (0.106 mm)	−24.3 dB

Each decade of \(\varepsilon\) shortens the pulse by 6–22 % (11–17 % from \(10^{-3}\) down) and costs about 8 dB of SNR once the inversion reaches beyond the −10 dB band edges (5 to 7 dB in the first decades). At large \(\varepsilon\) the Wiener filter is nearly a matched filter and actually gains SNR. Two reflectors 0.15 µs apart need an effective width near 0.15 µs, \(\varepsilon \approx 10^{-4}\) to \(10^{-5}\), which is affordable only with an input peak SNR of roughly 40 dB.

Demo 03.7

An RF line through a band-pass filter, and back through an inverse filter

Approx. model1D pulse-echo convolution: Gaussian 5 MHz pulse (60 % bandwidth), five point reflectors, seeded white Gaussian noise; no attenuation, diffraction or TGC. Filtering and Wiener deconvolution are computed exactly on that model (512-point FFT).

Time domain (µs at \(f_s = 40\) MHz)

Frequency domain (dB re peak, 0–20 MHz)

reflectivity \(\rho\) / pulse spectrumnoisy RF line \(r\)filter \(H_{\mathrm{bp}}\) or \(G\)outputper-reflector copies

Operation

Pass band

Noise std (re pulse peak) 0.10

log₁₀(ε/|P|²_max) −2.0

Show per-reflector copies (band-pass)

Try: in band-pass mode the gold dashed copies (\(h_{\mathrm{bp}} \ast p\) shifted to each reflector and scaled) add up to the noise-free output; the two reflectors 0.15 µs apart (0.115 mm, half the axial resolution) merge into one wavelet. Switch to deconvolution and lower \(\varepsilon\): the readout shows the effective pulse narrowing and the SNR falling; at \(10^{-7}\) the estimate drowns in amplified out-of-band noise. With noise 0.10 (20 dB peak SNR) no \(\varepsilon\) separates the close pair cleanly; set the noise to 0.01 and \(\varepsilon\) to \(10^{-4.5}\) and the pair splits into two opposite-sign wavelets. The best \(\varepsilon\) is set by the noise, not by the pulse.

Ultrasound link

No filter recovers information at frequencies where \(P \approx 0\): deconvolution can at best flatten the band the transducer delivers, so axial resolution after deconvolution is still bounded by that band (Module 13 treats this as regularized inversion). ULM sidesteps the limit differently: it does not invert the PSF, it localizes isolated echoes, whose position can be estimated far more finely than the PSF width.

03.9Visual summary

Each card is a question an expert answers in five seconds, the answer, and the picture to keep. Every curve is computed from the formula in its caption; the badge on each card says which.

03.10How experts think about it

An expert treats any linear, time-invariant block as one sequence \(h\), or equivalently one complex number per frequency, \(H(e^{j\omega})\). In time, output is a sum of shifted copies of \(h\); in frequency, each sinusoid is scaled by \(\lvert H\rvert\) and rotated by \(\angle H\), and nothing else can happen. Every practical question (how long, how much delay, how much noise, can it be undone, is it stable) is answered by reading \(h\), \(\lvert H\rvert\), \(\angle H\) or the pole-zero picture.

Linearity and time invariance are two separate tests. TGC is linear but time-varying; a squarer is time-invariant but nonlinear.

\(T\{ax_1+bx_2\} = aT\{x_1\}+bT\{x_2\}\)\(T\{x[n-k]\} = y[n-k]\)

Only systems that pass both, (03.1) and (03.2), are described by one impulse response. With \(g[n] = 1 + n/8\), a unit echo at \(n = 2\) comes out 1.25 high and the same echo at \(n = 10\) comes out 2.25 high: the output to a shifted input is not the shifted output. Figure: TGC response to \(\delta[n-2]\) shifted by 8 (pink circle) versus its response to \(\delta[n-10]\) (blue). Zero-in-zero-out is the quickest necessary check for linearity; it catches offsets and DC bias.

Exact formula

Real hardware is linear only below saturation and clipping; "time-invariant" also assumes the system does not drift during the record.

An LTI system is its impulse response: knowing \(T\{\delta\}\) you know \(T\{x\}\) for every \(x\).

\(h = T\{\delta\}\)\(y = x \ast h\)

\(x = \sum_k x[k]\,\delta[n-k]\), time invariance moves \(h\), linearity adds the copies, (03.4). This is why transducers are characterized by a pulse-echo impulse response from a planar reflector, and why an imaging system is characterized by its PSF. Figure: a unit impulse in (blue), the 9-tap RF pulse out (pink).

Exact formula

A beamformed image is linear but not shift-invariant in space (the PSF changes with depth and angle), so one PSF describes only a small isoplanatic patch.

Convolution is superposition: every input sample launches a scaled, shifted copy of h, and the output is their sum.

\([1,2,3] \ast [2,1,-1] = [2,5,7,1,-3]\)

Rows of the by-hand table are the copies \(1\cdot[2, 1, -1]\), \(2\cdot[2, 1, -1]\) shifted by one, \(3\cdot[2, 1, -1]\) shifted by two; column sums are \(y\). The same picture is the RF line: one echo per reflector, overlapping where reflectors are closer than the pulse length. Figure: the three copies (gold, pink, green) and their sum (open circles).

Exact formula

Assumes the copies add linearly; multiple scattering and nonlinear propagation add terms that are not copies of \(h\).

Lengths add: L + M − 1 samples out. Widths add: a 6-sample rect convolved with a 4-sample rect is a 9-sample trapezoid with a 3-sample flat top.

\(L + M - 1\)flat top \(L - M + 1\)

The first copy starts at \(n = 0\) and the last starts at \(L - 1\) and lasts \(M\) samples, (03.6). The flat top is the range of shifts where the short rect lies entirely inside the long one: \(6 - 4 + 1 = 3\) samples at height 4. In continuous time the support of \(x \ast h\) is the sum of the supports, so a filter always lengthens a pulse by its own length. Figure: [1 1 1 1 1 1] ∗ [1 1 1 1] = [1 2 3 4 4 4 3 2 1].

Exact formula

Support adds exactly; effective (rms) widths add in quadrature for positive, bell-shaped pulses (variances add), not linearly.

Sanity check any convolution: the sum of the output equals the product of the sums.

\(\sum y = \sum x \cdot \sum h\)12 = 6 × 2

It is the convolution theorem (03.19) at \(\omega = 0\): \(Y(e^{j0}) = X(e^{j0})\,H(e^{j0})\), and the DTFT at zero frequency is just the sum. A smoothing filter normalized to \(\sum h = 1\) preserves the mean; a high-pass with \(\sum h = 0\) (first difference, wall filter) removes it exactly. Figure: \(\lvert X\rvert\), \(\lvert H\rvert\), \(\lvert Y\rvert\) for the by-hand pair, with the \(\omega = 0\) values 6, 2 and 12 marked.

Exact formula

Checks only the DC component; a wrong sign or a wrong shift with the right sums passes it. Check one more frequency (\(\omega = \pi\): alternating sums) to catch those.

Complex exponentials are eigenfunctions: an LTI system can only scale and rotate each frequency, never create a new one.

\(e^{j\omega n} \to H(e^{j\omega})\,e^{j\omega n}\)\(\cos \to \lvert H\rvert\cos(\omega n + \angle H)\)

Pull \(e^{j\omega n}\) out of \(\sum_k h[k]\,e^{j\omega(n-k)}\), (03.9). New spectral lines therefore prove nonlinearity or time variation: a squarer maps \(\cos\omega_0 n\) to \(\tfrac12 + \tfrac12\cos 2\omega_0 n\), the mechanism behind the second harmonic of nonlinear propagation. Figure: input line at \(\omega_0 = 0.3\pi\) (gray); the LTI filter \(h = [\tfrac12, \tfrac12]\) only changes its height to \(\cos(0.15\pi) = 0.89\) (blue); the squarer moves energy to 0 and \(2\omega_0\) (pink).

Exact formula

Holds in steady state; a sinusoid switched on at \(n = 0\) has a start-up transient lasting about the length of \(h\) (or the IIR time constant).

Re H is the in-phase gain; −Im H is the gain of a copy shifted a quarter period late.

\(y = \mathrm{Re}\,H\cos\omega n - \mathrm{Im}\,H\sin\omega n\)\(\mathrm{Im}\,H < 0 \to\) lag

Take the real part of \((\mathrm{Re}\,H + j\,\mathrm{Im}\,H)(\cos\omega n + j\sin\omega n)\), (03.13). The output is real; \(j\) only means "rotate by 90°". For \(h = [\tfrac12, \tfrac12]\) at \(\omega = \pi/2\), \(\mathrm{Re}\,H = 0.5\) and \(\mathrm{Im}\,H = -0.5\), so \(y = 0.5\cos(\pi n/2) + 0.5\sin(\pi n/2)\) \(= 0.707\cos(\pi n/2 - 45^\circ)\). Figure: output stems (green) equal the sum of the in-phase (blue) and quadrature (pink) parts.

Exact formula

Re and Im depend on where \(n = 0\) is placed: delaying \(h\) by \(d\) multiplies \(H\) by \(e^{-j\omega d}\) and reshuffles them, while \(\lvert H\rvert\) is unchanged.

Even impulse response gives real H (zero phase); odd gives purely imaginary H (±90°). Causal filters always have phase.

\(h_e \leftrightarrow \mathrm{Re}\,H\)\(h_o \leftrightarrow j\,\mathrm{Im}\,H\)Hilbert: \(-j\,\mathrm{sgn}\,\omega\)

\(\cos\omega k\) is even in \(k\) and \(\sin\omega k\) is odd, so odd taps cancel in \(\mathrm{Re}\,H\) and even taps cancel in \(\mathrm{Im}\,H\), (03.15). Centered smoother (¼, ½, ¼): \(H = \tfrac12 + \tfrac12\cos\omega\), real. Centered difference (½, 0, −½): \(H = j\sin\omega\), pure quadrature. A causal symmetric FIR of \(M\) taps is a centered one delayed by \((M-1)/2\), so \(H = e^{-j\omega(M-1)/2} A(\omega)\) with \(A\) real: linear phase. Figure: \(\mathrm{Re}\,H\) of the smoother (blue) and \(\mathrm{Im}\,H\) of the difference (pink).

Exact formula

Zero phase needs non-causal processing (stored data, forward-backward filtering), which squares \(\lvert H\rvert\) and doubles every dB figure.

Group delay is the slope of the phase: the envelope arrives τ_g = −dφ/dω late. A symmetric M-tap FIR delays everything by (M − 1)/2.

41 taps → 20 samples0.5 µs = 0.385 mm

A narrowband pulse at \(\omega_0\) sees \(\varphi(\omega) \approx \varphi(\omega_0) + \varphi'(\omega_0)(\omega - \omega_0)\): the constant term shifts the carrier (phase delay \(-\varphi/\omega\)), the slope shifts the envelope, (03.16). At 40 MHz a 41-tap band-pass delays by 0.5 µs, which is 0.385 mm of depth (\(c/2 \times 0.5\) µs) and must be removed before depths are assigned. Figure: the 5 MHz pulse in (gray) and out of the 3–7 MHz band-pass (green), shifted by exactly 20 samples.

Approx. modelGaussian 5 MHz pulse (60 % bandwidth) as the test input; the 41-tap Hamming FIR is applied exactly.

For non-linear phase (IIR, minimum phase) \(\tau_g\) varies across the band and the pulse is distorted; group delay can even be negative locally without violating causality.

Convolution in time is multiplication in frequency: magnitudes multiply, phases add.

\(Y = X H\)\((-2-2j)(3-j) = -8-4j\)

Each phasor in \(x\) is multiplied by \(H\) at its own frequency, (03.19). At \(\omega = \pi/2\) for the by-hand pair: \(\lvert X\rvert\lvert H\rvert = 2.83 \times 3.16 = 8.94 = \lvert Y\rvert\) and −135° − 18.4° = −153.4° = \(\angle Y\). In rectangular form \(\mathrm{Re}\,Y = \mathrm{Re}\,X\,\mathrm{Re}\,H - \mathrm{Im}\,X\,\mathrm{Im}\,H\); the real parts alone do not multiply. Figure: rect 6 (blue) times rect 4 (pink) spectra give the trapezoid spectrum (green); its nulls are the union of both sets of nulls.

Exact formula

With DFTs it is circular convolution; it equals linear convolution only with padding to \(N \ge L + M - 1\).

FFT convolution is circular: pad to N ≥ L + M − 1 or the tail wraps onto the start.

\(N = 4\): \([-1, 5, 7, 1]\)\(N \ge 5\): \([2, 5, 7, 1, -3]\)

Sampling the spectrum at \(N\) points makes the time signal periodic with period \(N\), so \(y_{\mathrm{circ}}[n] = \sum_r y[n + rN]\), (03.8). For the by-hand pair at \(N = 4\), \(y[4] = -3\) adds to \(y[0] = 2\). For long RF lines use overlap-add with blocks of \(B\) samples and FFT size \(\ge B + M - 1\). Figure: linear result (open circles) and the \(N = 4\) circular result (filled pink).

Exact formula

Circular convolution is exactly right when the data really are periodic (a full rotation of angles, a ring array) or when the edges are discarded anyway.

An L-point moving average is a low-pass with nulls at multiples of f_s/L, delay (L − 1)/2 and sidelobes near −13 dB.

\(\lvert H\rvert = \left\lvert \frac{\sin(L\omega/2)}{L\sin(\omega/2)} \right\rvert\)\(L = 5\): −12.0 dB

The sum of \(L\) unit arrows rotated by \(-\omega\) each closes into a polygon (zero resultant) when \(L\omega\) is a multiple of \(2\pi\). The first sidelobe is −9.5 dB for \(L = 3\), −12.0 dB for \(L = 5\), −12.8 dB for \(L = 8\), approaching the −13.3 dB of a sinc. Half-power point for \(L = 5\): \(\omega = 0.566\) rad/sample. Use it to null a known interference at \(f_s/L\), not as a good low-pass. Figure: \(\lvert H\rvert\) for \(L = 3\) and 5.

Exact formula

Poor stopband; a windowed or equiripple FIR of the same length does far better (Module 09).

Differencing is a high-pass and a derivative: |H| = 2 sin(ω/2) ≈ ω, with +90° of phase at low frequency.

\(1 - e^{-j\omega} = 2j\sin(\omega/2)\,e^{-j\omega/2}\)\(H(\pi/2) = 1 + j\)

The factor \(j\) is the 90° lead of a derivative (\(\frac{d}{dt} e^{j\omega t} = j\omega\, e^{j\omega t}\)), \(e^{-j\omega/2}\) is the half-sample offset of the two taps. At \(\omega = \pi/2\), \(\lvert H\rvert = \sqrt2\) (+3 dB) and \(\angle H = 45^\circ\). \(\sum h = 0\), so DC is removed exactly: the simplest wall filter. Figure: \(\lvert H\rvert = 2\sin(\omega/2)\) (blue) against the ideal differentiator \(\lvert j\omega\rvert = \omega\) (dashed).

Exact formula

The approximation to \(\omega\) is within 10 % only below \(\omega \approx \pi/2\) (\(2\sin(\pi/4) = 1.41\) vs 1.57); differencing boosts high-frequency noise by up to 6 dB at Nyquist.

A one-pole smoother y = a y[n−1] + (1 − a)x has time constant ≈ 1/(1 − a) samples, bandwidth ≈ 1 − a rad/sample, and cuts white-noise variance to (1 − a)/(1 + a).

\(a = 0.9\): \(\tau = 9.5\), \(\omega_{3\mathrm{dB}} = 0.105\)\(a = 0.8 \approx\) 9-frame average

\(h[n] = (1-a)\,a^n\) decays with time constant \(-1/\ln a\) (9.49 samples at \(a = 0.9\)), (03.22). \(\lvert H\rvert^2 = (1-a)^2/(1 - 2a\cos\omega + a^2)\) falls to ½ at \(\cos\omega = 1 - (1-a)^2/(2a)\): 0.105 rad/sample at \(a = 0.9\). \(\sum h^2 = (1-a)/(1+a)\) is the noise variance gain: 1/9 at \(a = 0.8\), 1/19 at \(a = 0.9\). Group delay at DC is \(a/(1-a)\): 9 samples at \(a = 0.9\). Figure: \(\lvert H\rvert\) for \(a = 0.8\) (blue) and 0.9 (pink).

Exact formula

The approximations need \(a\) close to 1; the lag (group delay) grows with the smoothing, which smears fast changes (motion, contrast wash-in).

Stable means Σ|h| < ∞; for a causal rational filter, every pole strictly inside the unit circle.

\(\lvert p\rvert < 1\)\(a = 1.05\): doubles every 14.2 samples

A pole \(p\) contributes a term \(p^n\) to \(h[n]\); it decays if \(\lvert p\rvert < 1\) and grows if \(\lvert p\rvert > 1\). With \(\lvert p\rvert = 1\) the response neither decays nor grows, and a bounded input at the pole frequency drives the output without bound. Bounded input, bounded output holds iff \(h\) is absolutely summable, (03.24). Figure: the pole term \(p^n\) for \(p = 0.9\) (blue, decays) and \(p = 1.05\) (pink, grows).

Exact formula

Poles very close to the circle are stable on paper and fragile in fixed point; coefficient quantization can push them out (use biquads, Module 09).

|H| is the product of distances to the zeros over the product of distances to the poles: a pole near the circle is a peak, a zero on it is a null.

pole \(0.9\,e^{\pm j\pi/4}\): \(\lvert H\rvert = 7.43\) at \(\pi/4\)+12.0 dB over DC

At \(\omega = \pi/4\) the near pole is 0.100 away and the far one 1.345, so \(\lvert H\rvert = 1/(0.100 \times 1.345) = 7.43\); at DC both are 0.733 away, \(\lvert H\rvert = 1.86\), (03.23). The ratio is 12.0 dB. Peak width (−3 dB, full) \(\approx -2\ln r \approx 2(1 - r)\) rad/sample. You read the response off the picture without computing a transform. Figure: \(\lvert H\rvert\) (dB re DC) for that pole pair.

Exact formula

Exact for the magnitude; the phase also needs the angle sum and the \((P - Q)\,\omega\) term from the factors of \(z\). The width rule needs \(r\) near 1 and the pole away from \(\omega = 0, \pi\).

Dividing by H amplifies noise by 1/|H|. Regularize: the Wiener filter P*/(|P|² + ε) never exceeds 1/(2√ε).

10 MHz: \(\lvert P\rvert\) = −67 dB\(\varepsilon = 10^{-3}\): max gain +24 dB

\(\hat S = Y/P = S + W/P\), and \(W/P\) is huge where the transducer passes nothing, (03.27). The Wiener filter (03.28) is the inverse inside the band and rolls off outside it; maximizing \(\lvert P\rvert/(\lvert P\rvert^2 + \varepsilon)\) over \(\lvert P\rvert\) gives \(1/(2\sqrt\varepsilon)\) at \(\lvert P\rvert = \sqrt\varepsilon\). Figure: gain of the naive inverse \(1/\lvert P\rvert\) (gray, clipped at +60 dB) and Wiener gains for \(\varepsilon = 10^{-2}\) (blue) and \(10^{-4}\) (pink), all in dB re \(1/\lvert P\rvert_{\max}\).

Approx. modelGaussian pulse spectrum, 5 MHz centre, 60 % −6 dB bandwidth; gains computed exactly from (03.28).

Wiener is optimal for Gaussian signal and noise with known spectra; sparse reflectivity priors (Module 13) can do better, but nothing recovers frequencies with no signal.

03.11Cheat sheet

Concept	Formula	Ultrasound use
Linearity	\(T\{ax_1 + bx_2\} = aT\{x_1\} + bT\{x_2\}\)	RF chain before envelope; fails under clipping, nonlinear propagation
Time invariance	\(T\{x[n-k]\} = y[n-k]\)	Fails with TGC and depth-dependent attenuation
Impulse response	\(h = T\{\delta\}\)	Pulse-echo response, PSF
Convolution sum	\(y[n] = \sum_k x[k]\,h[n-k]\)	RF line = pulse ∗ reflectivity; filtering
Output length	\(L + M - 1\)	Record growth; FFT padding
Circular convolution	\(y_{\mathrm{circ}}[n] = \sum_r y[n + rN]\)	Pad FFTs to \(N \ge L + M - 1\); overlap-add
Eigenfunction	\(e^{j\omega n} \to H(e^{j\omega})\,e^{j\omega n}\)	Delay = phase ramp; LTI makes no harmonics
Frequency response	\(H(e^{j\omega}) = \sum_k h[k]\,e^{-j\omega k}\)	Filter gain and phase per frequency
Real and imaginary parts	\(y = \mathrm{Re}\,H\cos\omega n - \mathrm{Im}\,H\sin\omega n\)	In-phase vs quadrature response; Hilbert = \(-j\,\mathrm{sgn}\,\omega\)
Delays from phase	\(\tau_p = -\varphi/\omega\), \(\tau_g = -d\varphi/d\omega\)	Carrier vs envelope delay; FIR delay \((M-1)/2\)
Convolution theorem	\(Y = X H\); multiplication in time ↔ convolution in frequency	FFT filtering, leakage, apodization
Difference equation	\(y = \sum b_k x[n-k] - \sum a_k y[n-k]\), \(H(z) = B(z)/A(z)\)	FIR/IIR filters, wall filters, persistence
One-pole smoother	\(h = (1-a)\,a^n\), \(\tau \approx 1/(1-a)\), var × \((1-a)/(1+a)\)	Frame persistence, slow-time smoothing
Pole-zero geometry	\(\lvert H\rvert = \lvert g\rvert \prod\lvert e^{j\omega} - z_k\rvert / \prod\lvert e^{j\omega} - p_k\rvert\)	Read peaks and nulls off the z-plane
BIBO stability	\(\sum\lvert h\rvert < \infty \Leftrightarrow\) poles inside \(\lvert z\rvert = 1\) (causal, rational)	IIR design; quantization safety
Wiener deconvolution	\(G = P^{*}/(\lvert P\rvert^2 + \varepsilon)\), \(\max\lvert G\rvert = 1/(2\sqrt\varepsilon)\)	Axial deconvolution, noise control

03.12Interview and qualifying-exam questions

I1Is time-gain compensation an LTI operation? Is envelope detection? · screen

TGC, \(y[n] = g[n]\,x[n]\), is linear (scaling and adding inputs scales and adds outputs) but time-varying: an echo at a later time gets a larger gain, so delaying the input does not just delay the output. Envelope detection is time-invariant but nonlinear (the envelope of \(-x\) equals the envelope of \(x\)). Neither has a transfer function. Follow-up: "So where may you filter?" Anywhere before the envelope; and TGC commutes with an LTI filter only approximately, when \(g\) varies slowly over the filter length.

I2Convolve [1, 2, 3] with [2, 1, −1] by hand and give two checks. · screen

\([2, 5, 7, 1, -3]\). Checks: length \(3 + 3 - 1 = 5\); \(\sum y = 12 = \sum x \cdot \sum h = 6 \cdot 2\); and the alternating sum (\(\omega = \pi\)): \(2 - 5 + 7 - 1 - 3 = 0\) \(= (1 - 2 + 3)(2 - 1 - 1) = 2 \cdot 0\). Follow-up: "Do it with a 4-point FFT." You get \([-1, 5, 7, 1]\) because the −3 wraps onto \(n = 0\); pad to 5 or more.

I3Why are complex exponentials the eigenfunctions of LTI systems? · screen

\(y[n] = \sum_k h[k]\,e^{j\omega(n-k)}\) \(= e^{j\omega n}\sum_k h[k]\,e^{-j\omega k}\) \(= H(e^{j\omega})\,e^{j\omega n}\). The shift-invariance turns a delay into multiplication by \(e^{-j\omega k}\), and exponentials are the functions whose shifts are scalar multiples of themselves. Matrix view: convolution matrices are Toeplitz (circulant with padding), and the DFT vectors diagonalize every circulant. Follow-up: "Are sines eigenfunctions?" Not individually: a sine comes out as a scaled sine plus a cosine unless \(H\) is real at that \(\omega\).

I4What does the imaginary part of a frequency response mean physically? · onsite

For a cosine input, \(y = \mathrm{Re}\,H\cos\omega n - \mathrm{Im}\,H\sin\omega n\). \(\mathrm{Re}\,H\) is the gain of the in-phase copy, \(-\mathrm{Im}\,H\) the gain of a copy shifted 90° (a quarter period late if \(\mathrm{Im}\,H < 0\)). Example: \(h = [\tfrac12, \tfrac12]\) at \(\omega = \pi/2\) has \(H = 0.5 - 0.5j\), so the output is half in-phase plus half lagging quadrature, net 0.707 at −45°, a half-sample delay. \(\mathrm{Im}\,H\) comes from the odd part of \(h\) about \(n = 0\). Follow-up: "Which filter has \(\mathrm{Re}\,H = 0\) everywhere?" Any odd \(h\), e.g. the Hilbert transformer (\(-j\,\mathrm{sgn}\,\omega\)) or the centered difference (\(j\sin\omega\)).

I5Phase delay versus group delay: which one tells you where the echo envelope lands? · onsite

Group delay \(\tau_g = -d\varphi/d\omega\) at the pulse center frequency gives the envelope delay; phase delay \(-\varphi/\omega\) gives the carrier delay. For linear phase they coincide; for a 41-tap symmetric FIR at 40 MHz both are 20 samples = 0.5 µs, 0.385 mm of depth. Follow-up: "Why does this matter for beamforming?" A phase-delay mismatch between channel filters looks like a time-delay error: 10° at 5 MHz is 5.6 ns.

I6Derive the convolution theorem and say what changes for the DFT. · onsite

\(Y = \sum_n\sum_k x[k]\,h[n-k]\,e^{-j\omega n}\); write \(e^{-j\omega n} = e^{-j\omega k}\,e^{-j\omega(n-k)}\) and set \(m = n - k\): \(Y = X H\). For the N-point DFT the indices are mod \(N\), so DFT products give circular convolution; linear convolution needs \(N \ge L + M - 1\). Follow-up: "And the dual?" Multiplication in time is circular convolution of the DFTs divided by \(N\): windowing smears lines (leakage).

I7A one-pole smoother has a = 0.9. Give its impulse response, DC gain, −3 dB frequency, delay and stability. · onsite

\(h[n] = 0.1 \cdot 0.9^n\), \(\sum h = 1\) (DC gain 1). \(\lvert H\rvert^2 = 0.01/(1 - 1.8\cos\omega + 0.81) = \tfrac12\) at \(\cos\omega = 1 - 0.01/1.8\), \(\omega = 0.105\) rad/sample (\(\approx 1 - a\)). Time constant \(-1/\ln 0.9 = 9.5\) samples; group delay at DC \(a/(1-a) = 9\) samples. Pole at 0.9, inside the circle: stable. Noise variance gain \(0.1/1.9 = 1/19\). Follow-up: "Frame rate 50 Hz?" Time constant 0.19 s, cutoff \(0.105 \times 50/(2\pi) = 0.84\) Hz.

I8Why can't you deconvolve an RF line by dividing spectra, and what do you do instead? · research

\(\hat S = Y/P = S + W/P\): outside the transducer band \(\lvert P\rvert\) is tiny (for a 5 MHz, 60 % pulse, −67 dB at 10 MHz), so noise is amplified enormously. Use Wiener \(G = P^{*}/(\lvert P\rvert^2 + \varepsilon)\) with \(\varepsilon\) the noise-to-signal ratio, whose gain is capped at \(1/(2\sqrt\varepsilon)\), or sparse or Bayesian inversion with priors on \(\rho\). Follow-up: "Limits?" No method restores frequencies with no signal, \(P\) must be known (it changes with depth through attenuation, so blind or depth-adaptive estimates are used), and phase errors in \(P\) bias the result.

I9Is the pulse-echo RF line really an LTI convolution of a pulse with reflectivity? · research

Only locally. Frequency-dependent attenuation low-passes the pulse progressively with depth (center frequency downshift), diffraction changes the pulse off focus, TGC adds a time-varying gain, and multiple scattering and nonlinear propagation add non-convolutional terms. Practice: assume shift invariance within short windows (a few pulse lengths), estimate a local pulse, and compensate attenuation. Follow-up: "How does this affect sound-speed estimation?" Methods that compare phase or delays between windows or channels assume the same local pulse on both sides; depth-dependent spectral change biases them unless both are measured over matched windows.

I10Can a causal filter have a purely real frequency response? · research

Real \(H\) for all \(\omega\) requires \(h\) even about \(n = 0\) (for real \(h\)); a causal \(h\) is zero for \(n < 0\), so evenness forces \(h = c\,\delta[n]\). Any useful causal filter has phase. The best you can do causally is linear phase (symmetric FIR, pure delay). Offline, filter forward and then backward: \(H H^{*} = \lvert H\rvert^2\), zero phase, magnitude squared. Follow-up: "Price of zero phase?" Doubled dB figures, non-causal transients at both ends, and it cannot run in real time.

I11An IIR filter has a pole at z = 1.02. What do you see, and how do you fix it? · onsite

The impulse response contains \(1.02^n\), doubling every \(\ln 2/\ln 1.02 = 35\) samples: the output eventually overflows whatever the input. Fix: redesign so all poles are inside the circle, or if the pole moved because of coefficient quantization, implement as cascaded biquads (second-order sections) where each pole pair depends on two coefficients. Follow-up: "Pole exactly at 1?" An integrator: a bounded input of constant sign (DC) makes the output grow linearly without bound, so it is not BIBO stable; zero-mean noise makes it wander like a random walk.

03.13MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.

Download the scripts: m03_convolution_by_hand.m · m03_lti_eigenfunctions.m · m03_difference_equations.m · m03_rf_filter_deconv.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m03_convolution_by_hand.mThe by-hand pair four ways (scatter loop, gather loop, conv, convolution matrix), circular wrap with a short FFT, the theorem at ω = π/2, and timing · change x, h, N · runtime ≈ 1 s

%% m03_convolution_by_hand.m
% Purpose : compute one convolution four ways (input-side scatter loop,
%           output-side flip-and-slide loop, conv, convolution matrix),
%           show the wrap-around of a too-short FFT, check the convolution
%           theorem at one frequency, and time loop vs conv vs FFT.
% Change  : x, h, Nshort; Lr and the filter length in the timing cell.
% Expect  : all four methods print [2 5 7 1 -3]; the 4-point FFT gives
%           [-1 5 7 1]; X*H = Y = -8 - 4j at w = pi/2; sum(y) = 12 = 6*2.
clear; close all; rng(1);

x = [1 2 3];          % input
h = [2 1 -1];         % impulse response
L = numel(x); M = numel(h); Ny = L + M - 1;

%% 1. Input side (scatter): every x(k) launches a scaled, shifted copy of h
y1 = zeros(1, Ny);
for k = 1:L
    y1(k:k+M-1) = y1(k:k+M-1) + x(k) * h;
end

%% 2. Output side (gather): y(n) = sum_k x(k) h(n-k), h flipped and slid
y2 = zeros(1, Ny);
for n = 1:Ny
    for k = 1:L
        m = n - k + 1;                 % index into h (1-based)
        if m >= 1 && m <= M
            y2(n) = y2(n) + x(k) * h(m);
        end
    end
end

%% 3. Built-in conv, 4. convolution (Toeplitz) matrix
y3 = conv(x, h);
T = zeros(Ny, L);
for k = 1:L
    T(k:k+M-1, k) = h(:);              % column k = h shifted down by k-1
end
y4 = (T * x(:)).';
fprintf('scatter loop : %s\n', mat2str(y1));
fprintf('gather loop  : %s\n', mat2str(y2));
fprintf('conv         : %s\n', mat2str(y3));
fprintf('matrix T*x   : %s\n', mat2str(y4));
fprintf('length L+M-1 = %d, sum(y) = %g = sum(x)*sum(h) = %g\n', Ny, sum(y3), sum(x)*sum(h));
fprintf('alternating sum (w = pi): %g = %g\n', sum(y3 .* (-1).^(0:Ny-1)), ...
        sum(x .* (-1).^(0:L-1)) * sum(h .* (-1).^(0:M-1)));

%% Circular convolution: FFT products wrap when N < L+M-1
Nshort = 4;
yc4 = real(ifft(fft(x, Nshort) .* fft(h, Nshort)));
yc8 = real(ifft(fft(x, 8) .* fft(h, 8)));
fprintf('FFT, N = %d   : %s   (y(5) = -3 wrapped onto y(1))\n', Nshort, mat2str(round(yc4 * 1e10) / 1e10));
fprintf('FFT, N = 8   : %s\n', mat2str(round(yc8(1:Ny) * 1e10) / 1e10));

%% Convolution theorem at w = pi/2
dtft = @(s, w) sum(s .* exp(-1j * w * (0:numel(s)-1)));
w0 = pi/2;
X = dtft(x, w0); H = dtft(h, w0); Y = dtft(y3, w0);
fprintf('w = pi/2: X = %s, H = %s\n', num2str(X), num2str(H));
fprintf('          X*H = %s, Y from y = %s\n', num2str(X*H), num2str(Y));
fprintf('          |X||H| = %.3f = |Y| = %.3f; angle X + angle H = %.1f deg, angle Y = %.1f deg\n', ...
        abs(X)*abs(H), abs(Y), (angle(X)+angle(H))*180/pi, angle(Y)*180/pi);

%% Timing on an RF-sized problem (101-tap filter, 20000 samples)
Lr = 20000; Mf = 101;
r = randn(1, Lr); hb = randn(1, Mf);
tic; ya = zeros(1, Lr+Mf-1);
for m = 1:Mf                            % scatter over taps: Mf vector adds
    ya(m:m+Lr-1) = ya(m:m+Lr-1) + hb(m) * r;
end
t1 = toc;
tic; yb = conv(r, hb); t2 = toc;
Nf = 2^nextpow2(Lr+Mf-1);
tic; yc = real(ifft(fft(r, Nf) .* fft(hb, Nf))); yc = yc(1:Lr+Mf-1); t3 = toc;
fprintf('timing (s): tap loop %.4f, conv %.4f, FFT (N = %d) %.4f\n', t1, t2, Nf, t3);
fprintf('max |difference|: loop-conv %.1e, FFT-conv %.1e\n', max(abs(ya-yb)), max(abs(yc-yb)));

%% Figures
figure('Name', 'Convolution by hand');
cols = [0.80 0.55 0.10; 0.80 0.20 0.45; 0.10 0.55 0.40];
subplot(2,2,1); stem(0:L-1, x, 'filled'); hold on; stem(0:M-1, h, 'r'); hold off;
xlim([-0.5 4.5]); title('x (blue) and h (red)'); xlabel('n (samples)');
subplot(2,2,2); hold on;
for k = 1:L
    stem((k-1:k+M-2) + 0.12*(k-2), x(k)*h, 'Color', cols(k,:), 'filled');
end
hold off; xlim([-0.5 4.5]); title('copies x(k) h(n-k)'); xlabel('n (samples)');
subplot(2,2,3); stem(0:Ny-1, y3, 'k', 'filled'); hold on; stem(0:Nshort-1, yc4, 'r'); hold off;
xlim([-0.5 4.5]); title('linear (k), 4-pt circular (r)'); xlabel('n (samples)');
subplot(2,2,4); w = linspace(0, pi, 400);
Xw = arrayfun(@(v) dtft(x, v), w); Hw = arrayfun(@(v) dtft(h, v), w); Yw = arrayfun(@(v) dtft(y3, v), w);
plot(w/pi, abs(Xw), 'b', w/pi, abs(Hw), 'r', w/pi, abs(Yw), 'k', w/pi, abs(Xw.*Hw), 'g--');
xlabel('\omega/\pi'); ylabel('magnitude'); title('|X|,|H|,|Y|; |X||H| dashed');

What to look for: all four methods print [2 5 7 1 −3]; the 4-point FFT prints [−1 5 7 1]; \(X H\) and \(Y\) at \(\omega = \pi/2\) both print \(-8 - 4j\); \(\sum y = 12 = 6 \times 2\).

m03_lti_eigenfunctions.mFeeds cosines and complex phasors through filters, measures in-phase and quadrature gains by least squares, compares with Re H and −Im H; shows a squarer creating 2ω · change h or the frequencies · runtime ≈ 1 s

%% m03_lti_eigenfunctions.m
% Purpose : show that sinusoids pass through an LTI filter at the same
%           frequency, scaled by |H| and shifted by angle(H); measure the
%           in-phase and quadrature gains by least squares and compare them
%           with Re H and -Im H; show a complex phasor coming out as H times
%           itself; show that a squarer (nonlinear) creates new frequencies.
% Change  : the filters in hs, the test frequencies ws.
% Expect  : measured in-phase gain = Re H and quadrature gain = -Im H to
%           ~1e-15; for [1/2 1/2] at w = pi/2: 0.5 and 0.5 (H = 0.5 - 0.5j,
%           -45 deg). The squarer output has lines only at 0 and 2*w0.
clear; close all; rng(1);

hs    = {[0.5 0.5], ones(1,3)/3, [1 -1], 0.2*0.8.^(0:79)};
names = {'2-pt average', '3-pt average', 'first difference', 'one-pole a=0.8 (80 taps)'};
ws    = [0.1 0.25 0.5 0.75] * pi;          % rad/sample
n     = 0:399;                              % sample index
Hof   = @(h, w) sum(h .* exp(-1j * w * (0:numel(h)-1)));

%% Cosine in: fit y = a cos(wn) + b sin(wn) on the steady-state part
fprintf('%-26s %6s %9s %9s %9s %9s %8s %8s\n', 'filter', 'w/pi', 'Re H', 'a(fit)', '-Im H', 'b(fit)', '|H|', 'angH');
maxerr = 0;
for i = 1:numel(hs)
    h = hs{i}; M = numel(h);
    for w = ws
        x = cos(w * n);
        y = filter(h, 1, x);                % causal FIR, zero initial state
        keep = n >= M;                      % drop the start-up transient
        A = [cos(w * n(keep)).', sin(w * n(keep)).'];
        ab = A \ y(keep).';
        H = Hof(h, w);
        maxerr = max([maxerr, abs(ab(1) - real(H)), abs(ab(2) + imag(H))]);
        fprintf('%-26s %6.2f %9.4f %9.4f %9.4f %9.4f %8.4f %7.1f\n', names{i}, w/pi, ...
                real(H), ab(1), -imag(H), ab(2), abs(H), angle(H)*180/pi);
    end
end
fprintf('max |fit - prediction| over all cases: %.1e\n', maxerr);

%% Complex phasor in: y./x is the constant H (the eigenvalue)
h = [0.5 0.5]; w = pi/2;
x = exp(1j * w * n); y = filter(h, 1, x);
ratio = y(3:end) ./ x(3:end);
fprintf('phasor in, [1/2 1/2], w = pi/2: mean(y/x) = %s, std(y/x) = %.1e, H = %s\n', ...
        num2str(mean(ratio), 4), std(ratio), num2str(Hof(h, w), 4));
yb = filter(h, 1, cos(w*(-1:6)));          % x(-1) = 0 is included, so y(0..6) is exact
fprintf('  cos(pi n/2) in, n = 0..6: y = %s\n', mat2str(round(yb(2:end)*1e4)/1e4));

%% Nonlinear system: squarer creates DC and 2*w0
N = 200; w0 = 0.3*pi;                       % w0 on DFT bin 30 of 200
xs = cos(w0 * (0:N-1));
S1 = abs(fft(xs)) / N; S2 = abs(fft(xs.^2)) / N;
fprintf('cos(w0 n), w0 = 0.3 pi: |DFT|/N at bins 0, 30, 60: input %.3f %.3f %.3f, squared %.3f %.3f %.3f\n', ...
        S1(1), S1(31), S1(61), S2(1), S2(31), S2(61));
fprintf('  (squared: 1/2 at DC and 1/4 at +-2w0, i.e. cos^2 = 1/2 + 1/2 cos 2w0 n)\n');

%% Figures
wg = linspace(0, pi, 400);
figure('Name', 'Re H and Im H');
subplot(2,1,1);
H1 = arrayfun(@(v) Hof([0.5 0.5], v), wg); H3 = arrayfun(@(v) Hof([1 -1], v), wg);
plot(wg/pi, real(H1), 'b', wg/pi, -imag(H1), 'b--', wg/pi, real(H3), 'r', wg/pi, -imag(H3), 'r--');
xlabel('\omega/\pi'); ylabel('gain');
title('Re H (solid) and -Im H (dashed): 2-pt average (b), first difference (r)'); grid on;
subplot(2,1,2);
h = [0.5 0.5]; w = pi/2; nn = 0:12; tt = linspace(0, 12, 400); H = Hof(h, w);
plot(tt, cos(w*tt), 'Color', [0.5 0.5 0.5]); hold on;
plot(tt, real(H)*cos(w*tt), 'b--', tt, -imag(H)*sin(w*tt), 'r--', tt, abs(H)*cos(w*tt + angle(H)), 'k');
stem(nn, abs(H)*cos(w*nn + angle(H)), 'k', 'filled'); hold off;
xlabel('n (samples)'); title('out (k) = Re H cos (b) - Im H sin (r); input gray');

What to look for: measured in-phase gain equals \(\mathrm{Re}\,H\) and measured quadrature gain equals \(-\mathrm{Im}\,H\) to about \(10^{-15}\); for \([\tfrac12, \tfrac12]\) at \(\omega = \pi/2\), 0.5 and 0.5 (\(H = 0.5 - 0.5j\), −45°). The ratio \(y/x\) for a complex phasor is constant in \(n\). The squarer's spectrum has lines at 0 and \(2\omega_0\) only.

m03_difference_equations.mRuns moving average, first difference and one-pole smoother by explicit loops and by filter, prints the by-hand samples, −3 dB points, group delays, noise gains, and an unstable pole · change a, L · runtime ≈ 1 s

%% m03_difference_equations.m
% Purpose : run the three basic recursions (moving average, first
%           difference, one-pole smoother) by explicit loops and with
%           filter, print the by-hand samples, and check the standard
%           numbers: -3 dB points, group delay, noise variance gain,
%           moving-average sidelobes, instability, pole-pair geometry.
% Change  : a (one-pole coefficient), L (moving-average length).
% Expect  : step response 0.2 0.36 0.488 0.5904 0.67232; a = 0.9: -3 dB at
%           0.1055 rad/sample, group delay 9.00 samples at DC, noise gain
%           0.0526; a = 1.05 doubles every 14.21 samples; the pole pair
%           0.9 exp(+-j pi/4) peaks 12.04 dB above DC.
clear; close all; rng(1);

%% One-pole smoother y(n) = a y(n-1) + (1-a) x(n), by loop and by filter
a = 0.8; N = 40;
x = ones(1, N);                              % step input
y = zeros(1, N); yprev = 0;
for n = 1:N
    y(n) = a * yprev + (1 - a) * x(n);
    yprev = y(n);
end
yf = filter(1 - a, [1 -a], x);
fprintf('one-pole a = 0.8, step response y(0..5): %s\n', mat2str(round(y(1:6)*1e5)/1e5));
fprintf('  closed form 1 - a^(n+1)            : %s\n', mat2str(round((1 - a.^(1:6))*1e5)/1e5));
fprintf('  loop vs filter max diff: %.1e; time constant -1/ln(a) = %.2f samples\n', max(abs(y - yf)), -1/log(a));
hi = filter(0.5, [1 -0.5], [1 zeros(1, 5)]);
fprintf('one-pole a = 0.5, impulse response: %s\n', mat2str(hi));

%% Moving average and first difference, by loop and by filter
xr = [zeros(1,3) ones(1,4) zeros(1,4)];      % a rect
ma = zeros(size(xr)); fd = zeros(size(xr));
xa = @(k) (k >= 1) * xr(max(k, 1));          % x with zeros before the start
for n = 1:numel(xr)
    ma(n) = (xa(n) + xa(n-1) + xa(n-2)) / 3;
    fd(n) = xa(n) - xa(n-1);
end
fprintf('rect                 : %s\n', mat2str(xr));
fprintf('3-pt average (loop)  : %s\n', mat2str(round(ma*1e4)/1e4));
fprintf('3-pt average (filter): %s\n', mat2str(round(filter(ones(1,3)/3, 1, xr)*1e4)/1e4));
fprintf('first difference     : %s  (edges only)\n', mat2str(fd));
fprintf('conv([1 -1], [1 1 1 1]) = %s\n', mat2str(conv([1 -1], [1 1 1 1])));

%% Standard numbers for the one-pole smoother
wg = linspace(0, pi, 200001);
for a = [0.8 0.9]
    H = (1 - a) ./ (1 - a * exp(-1j * wg));
    w3 = acos(1 - (1 - a)^2 / (2 * a));          % closed form
    [~, i3] = min(abs(abs(H) - 1/sqrt(2)));      % numerical check
    ph = unwrap(angle(H)); gd0 = -(ph(2) - ph(1)) / (wg(2) - wg(1));
    hh = (1 - a) * a.^(0:2000);
    fprintf('a = %.1f: -3 dB at %.4f (grid %.4f) rad/sample; |H(pi)| = %.4f (%.1f dB); gd(0) = %.2f = a/(1-a) = %.2f\n', ...
            a, w3, wg(i3), abs(H(end)), 20*log10(abs(H(end))), gd0, a/(1-a));
    fprintf('        noise variance gain sum(h^2) = %.4f = (1-a)/(1+a) = %.4f; tau = %.2f samples\n', ...
            sum(hh.^2), (1-a)/(1+a), -1/log(a));
end

%% Moving-average sidelobes and the first-difference approximation
for L = [3 5 8]
    Hm = abs(sin(L*wg/2) ./ (L*sin(wg/2))); Hm(1) = 1;
    fprintf('L = %d moving average: first sidelobe %.2f dB, -3 dB at %.3f rad/sample\n', ...
            L, 20*log10(max(Hm(wg > 2*pi/L))), wg(find(Hm < 1/sqrt(2), 1)));
end
fprintf('first difference at w = pi/2: |H| = %.4f, angle = %.1f deg, 2 sin(w/2)/w = %.3f\n', ...
        abs(1 - exp(-1j*pi/2)), angle(1 - exp(-1j*pi/2))*180/pi, 2*sin(pi/4)/(pi/2));

%% Instability: pole at 1.05
hu = filter(1, [1 -1.05], [1 zeros(1, 40)]);
fprintf('a = 1.05: h(40) = %.2f, doubling time ln2/ln1.05 = %.2f samples\n', hu(41), log(2)/log(1.05));

%% Pole pair 0.9 exp(+-j pi/4): |H| from distances
p = 0.9 * exp(1j*pi/4);
A = real(poly([p conj(p)]));                     % [1, -2 r cos(theta), r^2]
Hp = 1 ./ polyval(A, exp(1j*wg)) .* exp(1j*wg).^2;   % H(z) = 1/(1 + a1 z^-1 + a2 z^-2)
e = exp(1j*pi/4);
fprintf('pole pair: distances at w = pi/4: %.3f, %.3f -> |H| = %.3f; at DC: %.3f -> |H| = %.3f\n', ...
        abs(e - p), abs(e - conj(p)), 1/(abs(e-p)*abs(e-conj(p))), abs(1-p), 1/abs(1-p)^2);
[pk, ip] = max(abs(Hp));
fprintf('  peak %.3f at w = %.3f pi, %.2f dB above DC; sum|h| = %.2f\n', pk, wg(ip)/pi, 20*log10(pk/abs(Hp(1))), ...
        sum(abs(filter(1, A, [1 zeros(1, 2000)]))));

%% Figures
figure('Name', 'Difference equations');
subplot(2,2,1); stem(0:39, yf, 'filled'); xlabel('n (samples)'); title('one-pole a = 0.8: step response');
subplot(2,2,2); stem(0:40, hu, 'r', 'filled'); hold on; stem(0:40, 0.9.^(0:40), 'b'); hold off;
xlabel('n (samples)'); title('p^n: 0.9 (b) decays, 1.05 (r) grows');
subplot(2,2,3);
H8 = 0.2 ./ (1 - 0.8*exp(-1j*wg)); H9 = 0.1 ./ (1 - 0.9*exp(-1j*wg));
plot(wg, abs(H8), 'b', wg, abs(H9), 'r', [0 pi], [1 1]/sqrt(2), 'k:'); xlim([0 1]);
xlabel('\omega (rad/sample)'); title('|H|: a = 0.8 (b), 0.9 (r)');
subplot(2,2,4); plot(wg/pi, 20*log10(abs(Hp)/abs(Hp(1)))); grid on;
xlabel('\omega/\pi'); ylabel('dB re DC'); title('pole pair 0.9e^{\pm j\pi/4}');

What to look for: step response 0.2, 0.36, 0.488, 0.5904, 0.67232; one-pole \(a = 0.9\): −3 dB at 0.105 rad/sample, group delay 9.00 samples at DC, noise variance gain 0.0526; \(a = 1.05\) doubles every 14.2 samples; \(\lvert H\rvert\) of the pole pair \(0.9\,e^{\pm j\pi/4}\) peaks 12.0 dB above DC.

m03_rf_filter_deconv.mA 5 MHz RF line with five reflectors: band-pass filtering as a sum of shifted filtered pulses, SNR gain for three pass bands, then naive versus Wiener deconvolution over ε · change the noise, the band, ε · runtime ≈ 2 s

%% m03_rf_filter_deconv.m
% Purpose : 1D pulse-echo line r = p * rho + noise at fs = 40 MHz. (1) Band-
%           pass filter it and show that the output is a sum of shifted,
%           scaled copies of the filtered pulse; compute the SNR gain for
%           three pass bands. (2) Wiener-deconvolve it and tabulate the
%           trade-off between effective pulse width and noise.
% Change  : sigma (noise std re pulse peak), the pass bands, the eps list.
% Expect  : copies-vs-filtered difference ~1e-16; expected SNR gains +5.8,
%           +7.3, +8.5 dB for 2-8, 3-7, 4-6 MHz; Wiener: max gain 1/(2 sqrt eps),
%           width 0.294 us without deconvolution, 0.214 us at eps = 1e-2,
%           0.155 us at 1e-4, with peak-SNR changes -0.6 and -15.7 dB.
clear; close all; rng(1);

fs = 40;                                    % sampling rate (MHz), time in us
f0 = 5; B = 3;                              % centre frequency, -6 dB bandwidth (MHz)
st = sqrt(2*log(2)) / (pi*B) * fs;          % Gaussian envelope std (samples)
k  = -20:20;                                % pulse support (samples)
p  = exp(-k.^2 / (2*st^2)) .* cos(2*pi*f0/fs*k);
N  = 256;                                   % 6.4 us line
pos = [50 90 96 150 200]; amp = [1 -0.7 0.6 0.5 0.8];
rho = zeros(1, N); rho(pos+1) = amp;
fprintf('pulse: sigma_t = %.3f samples, -6 dB envelope width %.3f us (%.3f mm)\n', st, ...
        2*sqrt(2*log(2))*st/fs, 2*sqrt(2*log(2))*st/fs*1.54/2);

sfull = conv(rho, p); s = sfull(21:20+N);   % centred: no delay
sigma = 0.1; w = sigma * randn(1, N); r = s + w;

%% Band-pass FIR (41-tap Hamming-windowed difference of sincs)
hw = 0.54 - 0.46*cos(2*pi*(0:40)/40);
sincd = @(f, m) (m == 0) * 2*f/fs + (m ~= 0) .* sin(2*pi*f/fs*m) ./ (pi*m + (m == 0));
bpf = @(f1, f2) hw .* (sincd(f2, k) - sincd(f1, k));
bands = [2 8; 3 7; 4 6];
for b = 1:3
    hb = bpf(bands(b,1), bands(b,2));
    yf = conv(r, hb);  y  = yf(21:20+N);    % centred output (20-sample delay removed)
    sf = conv(s, hb);  sf = sf(21:20+N);
    nf = conv(w, hb);  nf = nf(21:20+N);
    gS = sum(sf.^2) / sum(s.^2); gN = sum(hb.^2);
    fprintf('band %d-%d MHz: signal %+.2f dB, noise %+.2f dB, expected SNR gain %+.2f dB, this realization %+.2f dB\n', ...
            bands(b,1), bands(b,2), 10*log10(gS), 10*log10(gN), 10*log10(gS/gN), ...
            10*log10(sum(sf.^2)/sum(nf.^2)) - 10*log10(sum(s.^2)/sum(w.^2)));
    if b == 2, y37 = y; s37 = sf; hb37 = hb; end
end
q = conv(p, hb37);                           % filtered pulse, 81 taps, centre at index 41
copies = zeros(1, N);
for i = 1:numel(pos)
    idx = (pos(i)+1) + (-40:40);
    ok = idx >= 1 & idx <= N;
    copies(idx(ok)) = copies(idx(ok)) + amp(i) * q(ok);
end
fprintf('3-7 MHz: max |noise-free output - sum of shifted copies| = %.1e\n', max(abs(s37 - copies)));
fprintf('linear-phase delay of the 41-tap FIR: 20 samples = %.2f us = %.3f mm\n', 20/fs, 20/fs*1.54/2);

%% Wiener deconvolution: G = conj(P) / (|P|^2 + eps |P|max^2)
NF = 512;
pc = zeros(1, NF); pc(mod(k, NF) + 1) = p;    % pulse centred at sample 0 (circular)
P = fft(pc); Pm = max(abs(P));
fq = (0:NF-1) * fs / NF;
Pat = @(f) abs(sum(p .* exp(-1j*2*pi*f/fs*k)));   % DTFT at exactly f (MHz), not the nearest bin
fprintf('|P| re peak: %.1f dB at 2 MHz, %.1f dB at 8 MHz, %.1f dB at 10 MHz\n', ...
        20*log10(Pat(2)/Pat(5)), 20*log10(Pat(8)/Pat(5)), 20*log10(Pat(10)/Pat(5)));
R = fft(r, NF); Wn = fft(w, NF);
eps_list = [1e-1 1e-2 1e-3 1e-4 1e-5 1e-7];
fprintf('%9s %10s %12s %12s %12s %14s\n', 'eps', 'max gain', 'width (us)', 'width (mm)', 'dSNR (dB)', 'noise std out');
est = zeros(numel(eps_list), N);
for i = 1:numel(eps_list)
    G = conj(P) ./ (abs(P).^2 + eps_list(i)*Pm^2);
    qe = real(ifft(P .* G));                  % effective pulse
    Q = fft(qe); Q(2:NF/2) = 2*Q(2:NF/2); Q(NF/2+2:end) = 0;
    e = fftshift(abs(ifft(Q))); e = e / max(e);
    c = find(e == 1, 1); il = c; ir = c;
    while e(il) >= 0.5, il = il - 1; end
    while e(ir) >= 0.5, ir = ir + 1; end
    wl = il + (0.5 - e(il)) / (e(il+1) - e(il));  wr = ir - (0.5 - e(ir)) / (e(ir-1) - e(ir));
    width = wr - wl;                           % samples
    rmsG = sqrt(mean(abs(G).^2));               % white-noise gain of G
    no = real(ifft(Wn .* G)); 
    xr = real(ifft(R .* G)); est(i, :) = xr(1:N);
    fprintf('%9.0e %+9.1f dB %12.3f %12.3f %+12.1f %14.3f\n', eps_list(i), 20*log10(max(abs(G))*Pm), ...
            width/fs, width/fs*1.54/2, 20*log10(qe(1)/rmsG), std(no(1:N)));
end

%% Figures
t = (0:N-1) / fs;
figure('Name', 'RF line: band-pass filtering');
subplot(3,1,1); stem(pos/fs, amp, 'k', 'filled'); xlim([0 N/fs]); ylim([-1.2 1.2]); title('reflectivity \rho');
subplot(3,1,2); plot(t, r, 'b'); xlim([0 N/fs]); title('RF line r = p * \rho + noise');
subplot(3,1,3); plot(t, copies, 'Color', [0.8 0.6 0.1]); hold on; plot(t, y37, 'Color', [0.1 0.5 0.4]); hold off;
xlim([0 N/fs]); xlabel('time (\mus)'); title('3-7 MHz output (green) and noise-free sum of copies (gold)');
figure('Name', 'Spectra and Wiener deconvolution');
subplot(2,1,1); half = 1:NF/2+1;
Hb = abs(fft([hb37(21:41) zeros(1, NF-41) hb37(1:20)]));
plot(fq(half), 20*log10(abs(P(half))/Pm), 'k', fq(half), 20*log10(Hb(half)), 'r', ...
     fq(half), 20*log10(abs(R(half))/max(abs(R))), 'b'); ylim([-80 5]); grid on;
xlabel('f (MHz)'); ylabel('dB'); title('|P| (k), |H_{bp}| (r), noisy |R| (b)');
subplot(2,1,2); plot(t, est(2,:), 'Color', [0.1 0.5 0.4]); hold on; plot(t, est(4,:) - 0.6, 'm');
stem(pos/fs, 0.3*amp, 'k'); hold off; xlim([0 N/fs]);
xlabel('time (\mus)'); title('Wiener: eps 1e-2 (green), 1e-4 (magenta, -0.6); 0.3 rho (k)');

What to look for: noise-free filtered line minus the sum of shifted copies ≈ \(2\times10^{-16}\); \(\lvert P\rvert\) = −24.1 dB at 2 and 8 MHz and −67.0 dB at 10 MHz; expected SNR gains +5.8, +7.3 and +8.5 dB for 2–8, 3–7 and 4–6 MHz (a single noise realization lands within about 0.8 dB of them); the Wiener table reproduces the one in 03.8 (0.214 µs and −0.6 dB at \(\varepsilon = 10^{-2}\), 0.155 µs and −15.7 dB at \(10^{-4}\)) and the measured output noise std grows by a factor of 2.4 to 2.8 per decade of \(\varepsilon\).

03.14Self-check

Q1Is y[n] = n·x[n] linear? Time-invariant?

Linear (a gain applied sample by sample), not time-invariant (the gain depends on \(n\)). It is a TGC with \(g[n] = n\).

Q2Compute [1, −1] ∗ [1, 1, 1, 1]. What does the result tell you about differencing a rect?

\([1, 0, 0, 0, -1]\): length \(2 + 4 - 1 = 5\). Differencing a rect leaves only its two edges, +1 at the rising edge and −1 at the falling one; \(\sum y = 0 = \sum x \cdot \sum h = 0 \cdot 4\).

Q3For h = [½, ½], find Re H and Im H at ω = π/3 and the output to cos(πn/3).

\(H = \tfrac12(1 + e^{-j\pi/3})\) \(= \tfrac12(1 + 0.5 - 0.866j)\) \(= 0.75 - 0.433j\), \(\lvert H\rvert = 0.866\), \(\angle H = -30^\circ\) (half a sample at 60° per sample). \(y = 0.75\cos(\pi n/3) + 0.433\sin(\pi n/3)\) \(= 0.866\cos(\pi n/3 - 30^\circ)\).

Q4Which of these has a purely imaginary H: (a) h = [1, −1] at n = 0, 1; (b) h[−1] = 1, h[1] = −1; (c) h = [1, 2, 1] at n = −1, 0, 1?

(b): it is odd about \(n = 0\), \(H = e^{j\omega} - e^{-j\omega} = 2j\sin\omega\). (c) is even, \(H = 2 + 2\cos\omega\), real. (a) is the causal first difference, \(2j\sin(\omega/2)\,e^{-j\omega/2}\): imaginary only after removing the half-sample delay.

Q5You convolve a 2000-sample RF line with a 49-tap filter by FFT. What is the smallest FFT length, and the smallest power of two?

\(N \ge 2000 + 49 - 1 = 2048\), which is already a power of two. One sample less and the last output would wrap onto the first.

Q6A causal filter has poles at 0.95e^±j0.3π and zeros at ±1. Where is |H| largest, where is it zero, and is it stable?

Peak near \(\omega = 0.3\pi\) (pole distance about 0.05), nulls at \(\omega = 0\) and \(\omega = \pi\) (zeros on the circle at \(z = 1\) and \(z = -1\)): a band-pass resonator. Poles have radius \(0.95 < 1\): stable, ringing time about \(1/(1 - 0.95) = 20\) samples.

Q7Why does a Wiener deconvolution filter look like a matched filter at frequencies where the pulse is weak?

\(G = P^{*}/(\lvert P\rvert^2 + \varepsilon) \approx P^{*}/\varepsilon\) when \(\lvert P\rvert^2 \ll \varepsilon\). \(P^{*}\) is the matched filter, so where the data are mostly noise the Wiener filter correlates with the pulse and attenuates instead of inverting.

03.15Figures from the literature

Published figures to hold this page against: two textbooks that draw the same constructions step by step, and one paper with measured echoes. Links only (no images reproduced); every figure number was checked against the published text.

Input-side convolution drawn copy by copy

S. W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing, California Technical Publishing, 1997, Chapter 6 "Convolution", Figs. 6-5 and 6-6 (also 6-7). dspguide.com/CH6.PDF

A 9-point input convolved with a 4-point impulse response gives 12 output points (Fig. 6-5); Fig. 6-6 shows the nine scaled, shifted impulse responses whose sum is the output, the same construction as the by-hand table of 03.2 and Demo 03.2. Fig. 6-7 swaps the roles of input and impulse response and gets the same output: commutativity.

Output-side view: the "convolution machine" and end effects

Same book, Chapter 6, Figs. 6-8, 6-9 and 6-10. dspguide.com/CH6.PDF

Fig. 6-8 is a flow diagram of the flipped impulse response sliding along the input; Fig. 6-9 sets it to compute \(y[0]\), \(y[3]\), \(y[8]\) and \(y[11]\). Compare with Figure 03.B and the left panel of Demo 03.3. Fig. 6-10 shows why the first and last \(M - 1\) output points of an M-point filter may be unusable, the transient part of the \(L + M - 1\) rule.

Single-pole recursive filters, and zero, linear and nonlinear phase

Same book, Chapter 19 "Recursive Filters", Figs. 19-2, 19-5, 19-7 and 19-8. dspguide.com/CH19.PDF

Fig. 19-2: a single-pole low-pass smoothing a step edge (compare Demo 03.5). Fig. 19-5: single-pole frequency responses. Fig. 19-7 contrasts zero, linear and nonlinear phase impulse and step responses (the symmetric step response of zero phase); Fig. 19-8 shows bidirectional filtering, which makes any recursive filter zero phase at the price of \(\lvert H\rvert^2\) (rule card on even and odd \(h\), interview I10).

Deconvolution and the Wiener filter on a pulse in noise

Same book, Chapter 17 "Custom Filters", Figs. 17-4, 17-5 (deconvolution), 17-6 (deconvolution of old phonograph recordings) and 17-7 to 17-9 (optimal filters). dspguide.com/CH17.PDF

Figs. 17-4 and 17-5 design a filter that narrows the pulses of a gamma-ray detector (Fig. 17-3), in time and in frequency, the operation of (03.27)–(03.28). Fig. 17-6 applies deconvolution to real recordings. Figs. 17-8 and 17-9 compare moving-average, matched and Wiener filters on a signal whose spectrum overlaps the noise: the Wiener filter is designed from the relative amounts of signal and noise at each frequency, which is the role of \(\varepsilon\) here.

Two-pole resonator at θ = π/4 for several pole radii

J. O. Smith III, Introduction to Digital Filters with Audio Applications, W3K Publishing, 2007 (online edition), Appendix B, "Resonator Bandwidth in Terms of Pole Radius", Fig. B.6. ccrma.stanford.edu/~jos/fp

Amplitude (a) and phase (b) of \(y(n) = x(n)\) \(+\, 2R\cos(\theta_c)\,y(n-1)\) \(-\, R^2 y(n-2)\) with \(\theta_c = \pi/4\), the same pole angle as the resonator rule card and Figure 03.D. The page gives the bandwidth rule \(B = -\ln R/(\pi T)\) Hz, i.e. a full −3 dB width of \(-2\ln R \approx 2(1 - R)\) rad/sample. Below: the same equation re-plotted here for radii chosen on this page (not necessarily the book's).

Exact formulaRe-plotted from the difference equation given with Fig. B.6 (θ_c = π/4); R = 0.5 (gray), 0.8 (blue), 0.9 (pink), 0.95 (green), chosen here.

Measured echoes: pulse estimate, spectrum and Wiener-type deconvolution

R. H. Herrera, R. Orozco, M. Rodríguez, "Wavelet-based deconvolution of ultrasonic signals in nondestructive evaluation," arXiv:1202.0609, 2012. Figs. 6, 7, 8 and 10. doi:10.48550/arXiv.1202.0609

Real data: a 2.25 MHz unfocused transducer on a steel test block (STBA1), sampled at 50 MHz. Fig. 6 is the measured signal and its power spectrum, Fig. 7 the pulse estimated from it, Fig. 8 the deconvolution of three back-wall echoes and Fig. 10 the reflectivity estimated from three closely spaced echoes. Read it as 03.8 on measured data: the pulse must be estimated, the inverse is regularized (Wiener with wavelet denoising), and the gain in axial resolution is limited by the transducer band. Steel (5931 m/s) rather than tissue, so attenuation-driven pulse change with depth is small.

03.16Go deeper

A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson, 2010. Chapters 2 (LTI systems, convolution, eigenfunctions), 3 (z-transform) and 5 (phase, group delay, minimum phase).
A. V. Oppenheim, A. S. Willsky and S. H. Nawab, Signals and Systems, 2nd ed., Prentice Hall, 1997. The superposition derivation in continuous and discrete time, side by side.
S. W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing, California Technical Publishing, 1997. Input-side and output-side views of convolution, drawn step by step.
R. G. Lyons, Understanding Digital Signal Processing, 3rd ed., Prentice Hall, 2010. Practical moving-average, differentiator and one-pole filters.
R. N. Bracewell, The Fourier Transform and Its Applications, 3rd ed., McGraw-Hill, 2000. Convolution theorem, even and odd parts, and the pictures behind them.
N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, MIT Press, 1949. The origin of the Wiener filter.
J. A. Jensen, "Deconvolution of ultrasound images," Ultrasonic Imaging, vol. 14, 1992. Pulse-echo deconvolution and its limits.

← 02 The DFT and FFT by hand 04 Sampling, Fourier, LTI →