By now, even my mom knows classical wave-digital filter theory (WDF) as invented by Alfred Fettweis in 1971. What is less appreciated — even among people who use them — is that WDFs are the answer to a specific question the pure-DSP toolkit does not know how to phrase:

Which discrete-time realization of a passive continuous-time system remains passive at every sample, under composition, and under coefficient quantization?

Nothing in H(z)H(z)-land guarantees this. Poles inside the unit circle give BIBO stability, but stability is a spectral property — a statement about the impulse response integrated over all time — and it fails in exactly the ways a DSP engineer spends their career dodging: coefficient rounding drifts poles across the unit circle in narrow-band filters, feedback loops of two stable filters go unstable, nonlinear extensions of a stable linear system limit-cycle.

Passivity is the algebraic strengthening. Wave variables are the coordinate system that turn it into a pointwise identity you can check, not a global property you can only hope for. The rest of the theory is the machinery that lets you build networks in that coordinate system.

Classical

Waves and the reference impedance

At a port with voltage vv and current ii, pick a positive scalar R>0R > 0 and define

a=v+R i,b=v−R i.a = v + R\,i, \qquad b = v - R\,i.

RR is a reference impedance: a coordinate choice on the (v,i)(v,i) plane, not a property of any element. Changing RR changes what one calls "incident" vs "reflected" but not the physical trajectory. This freedom is the entire budget classical WDF has to spend.

Reflection coefficient and adaptation

For a linear one-port with Laplace impedance Z(s)Z(s), the wave-domain ratio is a Möbius transform of Z/RZ/R:

ρ(s)  =  b(s)a(s)  =  Z(s)−RZ(s)+R.\rho(s) \;=\; \frac{b(s)}{a(s)} \;=\; \frac{Z(s) - R}{Z(s) + R}.

The Möbius zero sits at Z=RZ = R. Adaptation spends the reference-impedance budget to place that zero exactly where ρ(z)\rho(z)'s constant-in-z−1z^{-1} term needs to vanish — after which b[n]b[n] carries no a[n]a[n] dependence and the port has no self-reflection at the current sample. Adaptation is not solving an algebraic loop; it is avoiding creating one.

Why bilinear

The location "where ρ(z)\rho(z)'s constant term needs to vanish" is determined by the discretization. The bilinear transform

s  =  2T 1−z−11+z−1s \;=\; \frac{2}{T}\,\frac{1 - z^{-1}}{1 + z^{-1}}

is the trapezoidal-rule discretization of d/dtd/dt, and it is the only rational s↔zs \leftrightarrow z map with all three of:

  1. LHP →\to open unit disk bijectively (stability preserved).
  2. Imaginary axis →\to unit circle bijectively (no frequency aliasing).
  3. On that boundary, ∣S(s)∣=1  ⟺  ∣S(z)∣=1|S(s)|=1 \iff |S(z)|=1 pointwise (losslessness preserved as an algebraic identity, not merely in integrated form).

Property 3 is the load-bearing one. It is why adaptation of a reactive element yields a pure delay in zz: trapezoidal integration is the unique choice under which the reflection rule of CC or LL becomes a finite polynomial in z−1z^{-1} instead of an infinite series.

Worked example: the capacitor

Constitutive law i=C v˙i = C\,\dot v. Trapezoidal integration between (n−1)T(n-1)T and nTnT:

v[n]−v[n−1]  =  T2C (i[n]+i[n−1]).v[n] - v[n-1] \;=\; \frac{T}{2C}\,\bigl(i[n] + i[n-1]\bigr).

Substitute v=(a+b)/2v = (a+b)/2, i=(a−b)/(2R)i = (a-b)/(2R):

(a[n]+b[n])−(a[n−1]+b[n−1])=  T2CR[(a[n]−b[n])+(a[n−1]−b[n−1])].\begin{aligned} &\bigl(a[n] + b[n]\bigr) - \bigl(a[n-1] + b[n-1]\bigr) \\ &\qquad =\; \frac{T}{2CR}\Bigl[\bigl(a[n] - b[n]\bigr) + \bigl(a[n-1] - b[n-1]\bigr)\Bigr]. \end{aligned}

Choose R=T/(2C)R = T/(2C) so the bracketed coefficient is 11. Every aa cancels and

  b[n]  =  a[n−1].  \boxed{\;b[n] \;=\; a[n-1].\;}

The capacitor, with its adapted reference impedance, is a unit-delay reflection: one memory cell, nothing else. For an inductor LL the analogous calculation gives R=2L/TR = 2L/T and b[n]=−a[n−1]b[n] = -a[n-1]. Both reactances become memory cells at the port; all dynamics get pushed into the interconnect.

Passivity vs. stability

At any port, at any sample, the instantaneous power reads

p[n]  =  v[n] i[n]  =  a[n]2−b[n]24R.p[n] \;=\; v[n]\,i[n] \;=\; \frac{a[n]^{2} - b[n]^{2}}{4R}.

Losslessness of an element is ∑portsp[n]≡0\sum_{\text{ports}} p[n] \equiv 0; passivity is ≤0\le 0. Both are algebraic identities in the wave variables, checkable per sample. Consequences the H(z)H(z) toolkit does not give:

  • BIBO stability is automatic. ∣b∣2≤∣a∣2|b|^{2} \le |a|^{2} per port implies ∑nb[n]2≤∑na[n]2\sum_n b[n]^{2} \le \sum_n a[n]^{2}; no output grows unboundedly.
  • Composition preserves it. Connecting two passive networks yields a passive network by matching waves at the shared port — no fresh stability proof.
  • Symmetric quantization preserves it. If wave amplitudes are rounded by a quantizer QQ with ∣Q(x)∣≤∣x∣|Q(x)| \le |x|, then ∣Q(b)∣2≤∣b∣2≤∣a∣2|Q(b)|^{2} \le |b|^{2} \le |a|^{2} — the standard failure mode of pole-inside-disk stability under coefficient rounding cannot occur.
  • Nonlinear extension preserves it. Any b=g(a)b = g(a) with ∣g(a)∣≤∣a∣|g(a)| \le |a| pointwise is passive and can be inserted at any port.

None of these hold for a Direct-Form-II biquad realization of the same transfer function.

What that actually buys over a biquad cascade

Two concrete gains:

  1. Coefficient sensitivity. A WDF derived from a doubly-terminated LC ladder inherits the ladder's coefficient sensitivity, which is minimal in the sense of Orchard (1966): sensitivity of the transfer function to component variations vanishes to first order in the passband. For narrow-band or high-order filters this is measurably 6–12 bits of fixed-point headroom over Direct Form II at the same target response.
  2. Nonlinear elements compose. A diode's i=Is(ev/VT−1)i = I_s(e^{v/V_T} - 1) becomes a wave map b=g(a)b = g(a) solvable at one port. Slotting it into an existing WDF is a local operation. Slotting the same diode into a biquad cascade is not a well-defined operation at all — biquads have no port to attach it to.

Delay-free loops

Adaptation zeros the z0z^0 term of ρ\rho for a single one-port. In a topology with multiple non-adapted ports connected through an adaptor, each port's zero condition wants a different reference impedance; the adaptor cannot satisfy all of them at once. The constant term reappears at the adaptor level — a delay-free loop across the network.

Newton or fixed-point iteration always solves it. The cost is threefold:

  • Constant time per sample dies. Iteration count is data-dependent, breaking real-time-scheduling analysis.
  • Convergence is not guaranteed for nonlinear elements. The contractivity condition depends on the port impedances, which are themselves the adaptation parameters — a coupled problem.
  • The passivity identity survives only at the fixed point. Capping iterations for real-time reasons loses the algebraic guarantee that motivated waves in the first place.

The classical framework treats DFLs as a topology failure to be fixed at model-build time — R-type multi-port adaptors, added delays where physically justified, reformulation — because every DFL chips at the load-bearing wall.


The three follow-up questions the framework leaves open are what the modern reformulations answer.

BWDF — decoupling numerical range from RR

The classical RR is fixed per element by adaptation. But numerical dynamic range of a,ba, b is a hardware constraint, not a physical one — adapted R=T/(2C)R = T/(2C) can be arbitrarily small or large, forcing wave amplitudes to occupy pathological ranges in fixed point.

Substitute (1, R)↦(ν, ζ)(1,\, R) \mapsto (\nu,\, \zeta) in the wave definition:

a=ν v+ζ i,b=ν v−ζ i.a = \nu\, v + \zeta\, i, \qquad b = \nu\, v - \zeta\, i.

Adaptation still cares only about the physical port resistance R=ζ/νR = \zeta/\nu; the pair (ν,ζ)(\nu, \zeta) has a genuine gauge symmetry

(ν,ζ)  ⟼  (λν,λζ),λ>0,(\nu, \zeta) \;\longmapsto\; (\lambda\nu, \lambda\zeta), \quad \lambda > 0,

that rescales wave amplitudes independently of RR.

The interesting fact is that two apparently distinct "good conditioning" requirements pick the same gauge. Demand:

  • (i) pseudopower (a2−b2)/(4νζ)(a^2 - b^2)/(4\nu\zeta) is RR-independent, and
  • (ii) under time-varying R(t)R(t), the geometric term in a˙\dot a contains only bb, not aa (no delay-free a→aa \to a coupling from modulation).

Both fix νζ≡c\nu\zeta \equiv c, giving

ν(R)=c/R,ζ(R)=cR\nu(R) = \sqrt{c/R}, \qquad \zeta(R) = \sqrt{cR}

— the power-wave section, uniquely. Time-varying RR, parametric filters, and state-dependent adaptation in nonlinear elements all live on this section for that reason.

VWDF — genuine multi-port coupling

Classical adaptors are series or parallel: KVL or KCL at a scalar junction. Circuits with intrinsic multi-port coupling — transformers, coupled inductors, gyrators, MIMO scattering structures — have no series/parallel decomposition. In classical WDF they surface as global DFLs of nn scalar ports each.

Substitute scalar →n\to n-vector, R→R∈SPDnR \to \mathbf R \in \mathrm{SPD}_n:

a=v+R i,b=v−R i.\mathbf a = \mathbf v + \mathbf R\,\mathbf i, \qquad \mathbf b = \mathbf v - \mathbf R\,\mathbf i.

Adaptation lifts to: choose R\mathbf R so that a designated block of the vector scattering matrix S\mathbf S vanishes. What was an nn-fold global DFL — resolved only by per-sample Newton iteration — becomes one build-time matrix inverse. The nonlinearity of the runtime problem is traded for the linear algebra of the model-build step.

Matrix WDF — both at once

The last substitution combines the previous two: (1,R)↦(A,B)(1, R) \mapsto (\mathbf A, \mathbf B) with A,B∈GLn\mathbf A, \mathbf B \in \mathrm{GL}_n,

a=A v+B i,b=A v−B i,R=B A−1.\begin{aligned} \mathbf a &= \mathbf A\,\mathbf v + \mathbf B\,\mathbf i, \\ \mathbf b &= \mathbf A\,\mathbf v - \mathbf B\,\mathbf i, \\ \mathbf R &= \mathbf B\,\mathbf A^{-1}. \end{aligned}

Gauge is GLn\mathrm{GL}_n acting by (A,B)↦(GA,GB)(\mathbf A, \mathbf B) \mapsto (\mathbf G\mathbf A, \mathbf G\mathbf B); the BWDF flat-section condition νζ≡c\nu\zeta \equiv c generalizes to A ⁣⊤ ⁣B=c I\mathbf A^{\!\top}\!\mathbf B = c\,\mathbf I, under which the pseudopower quadratic form diagonalizes as (∥a∥2−∥b∥2)/(4c)(\|\mathbf a\|^{2} - \|\mathbf b\|^{2})/(4c). The four levels collapse:

levelA\mathbf AB\mathbf Bgauge
Classical11RRnone
BWDFν\nuζ\zetaR>0\mathbb R_{>0}
VWDFI\mathbf IR\mathbf Rnone
MatrixA\mathbf AB\mathbf BGLn\mathrm{GL}_n

None of the four levels linearizes a genuinely nonlinear port, removes a delay-free cycle containing more than one non-adapted port, extends to active elements without indefinite-metric wave definitions, or handles distributed elements without prior discretization. The theory tells you exactly which problems it can and cannot algebraically eat.