EP38

EP38: Reverberation Is Matrix Multiplication

卷积混响, FDN反馈矩阵, Hadamard正交性, 互质延迟
5:56 Linear AlgebraSignal Processing

Overview

The equation appears throughout linear algebra. When is the Hadamard matrix and is a vector of audio delay-line signals, the equation is not an exercise — it is the core computational step inside every digital reverb plugin you have ever heard.

This episode builds the mathematical framework for artificial reverberation from first principles. The story moves through three levels: (1) convolution reverb, which is mathematically exact but computationally prohibitive at 96,000 multiply-accumulate operations per sample; (2) Schroeder’s 1962 comb-and-allpass architecture, the first practical alternative; and (3) the Feedback Delay Network (FDN) introduced by Jot and Chaigne in 1992, which reframes the entire problem as iterated matrix-vector multiplication.

The central mathematical question of the FDN is stability: what conditions on the feedback matrix ensure that the feedback loop decays rather than explodes? The answer is that must be unitary (all singular values equal to 1), guaranteeing that energy is not amplified in any direction of the delay-line state space. The Hadamard matrix , whose entries are only and , is unitary, maximally mixing, and computable using only additions and subtractions — no multiplications required.

The remaining design degree of freedom — the delay lengths — controls the spectral density of resonance modes. Choosing mutually coprime delay lengths distributes mode frequencies uniformly, preventing the metallic flutter artefacts that plagued early Schroeder reverbs.

中文: “y等于A乘以x。这个式子你见过很多次了。但如果A是一个哈达玛矩阵,x是一组音频延迟信号——这就不是线性代数的习题,这是你在任何一款数字混响插件里听到的空间感。这不是比喻。混响,就是矩阵乘法。”


Prerequisites

  • 全音程列与群结构(EP04) — group structure, modular arithmetic, and the algebraic regularity of Hadamard-type constructions; the rows of form a group under componentwise multiplication
  • Primes in the Studio — QRD Diffusers(EP29) — coprime and prime-based sequences in acoustics; the connection between number theory and uniform mode distribution; the Schroeder frequency separating discrete-mode and diffuse-field regimes

Definitions

Definition 38.1 (Convolution Reverb)

Given the measured impulse response of an acoustic environment (a finite-length sequence of samples obtained by recording a brief transient in the space), the convolution reverb output is

where is the dry input signal. This operation is called linear convolution. The impulse response encodes the full acoustic fingerprint of the space: early reflections, late reverberation tail, frequency-dependent absorption, and room geometry — all collapsed into a single time-domain sequence.

Computing directly costs multiply-accumulate operations per output sample. For a 2-second impulse response at a 48 kHz sample rate, , making naive real-time implementation infeasible on 1980s and early 1990s hardware. Using the FFT fast convolution identity , the complexity drops from to per block, enabling real-time convolution on modern hardware.

Definition 38.2 (Comb Filter)

A comb filter with delay length (samples) and feedback gain is the recursive filter defined by

Its transfer function is , with poles at for . The frequency response has peaks spaced apart — a comb pattern. The filter generates an exponentially decaying series of echoes at multiples of the delay .

An allpass filter is a comb filter augmented with a feedforward path chosen so that the magnitude response is flat ( for all ) while the phase response is frequency-dependent. The allpass filter modifies the temporal density of echoes without colouring the spectrum.

Definition 38.3 (Feedback Delay Network (FDN))

A Feedback Delay Network of order consists of:

  • delay lines of lengths (samples), with state vectors
  • a scalar output from each delay line:
  • an feedback matrix that mixes the delay-line outputs
  • per-channel loss (absorption) filters (scalars or IIR filters) inserted in each delay loop
  • input and output gain vectors

The update equation is

where is the vector of current delay-line outputs, and is fed back into the input of each delay line. The output sample is . In the lossless case (all ), the state update simplifies to .

Definition 38.4 (Unitary Matrix)

A real matrix is orthogonal (or, in the complex case, unitary) if

Equivalently, all singular values of equal 1, and for every — the matrix preserves the Euclidean norm. In the FDN context, an orthogonal feedback matrix ensures that the total energy is neither amplified nor reduced by the mixing step alone; decay comes entirely from the loss filters .

Definition 38.5 (Hadamard Matrix)

The Hadamard matrix of order is defined recursively by

The normalisation factor at each level accumulates to an overall factor of in the full matrix, ensuring orthogonality. The unnormalised version has entries only in and satisfies .

For :

The matrix-vector product requires only additions and subtractions — no multiplications — because all entries are . The full product costs via the Fast Hadamard Transform (FHT), analogous to the FFT.

Definition 38.6 (RT60 (Reverberation Time))

The RT60 (or ) of a reverberant system is the time, in seconds, required for the sound pressure level to decay by 60 dB after the source is switched off. In an FDN, each delay line accumulates a round-trip gain of per full traversal (where is the per-sample loss factor). Setting the total energy decay over seconds to (a factor of in power, or 60 dB):

where is the sample rate in Hz and is the delay length in samples. Longer delay lines require a smaller (more attenuation per sample) to achieve the same RT60 as shorter ones. Frequency-dependent RT60 is implemented by replacing the scalar with a low-order IIR absorption filter.


Main Theorems

Theorem 38.1 (FDN Stability via Unitary Feedback Matrix)
Consider a lossless FDN (all ) with feedback matrix . The system is BIBO stable (bounded input implies bounded output) if and only if every eigenvalue of the system’s Z-domain transfer matrix has magnitude at most 1. A sufficient condition for stability is that is orthogonal (). Under this condition the FDN is lossless: total delay-line energy is conserved at every time step (in the absence of input).
Proof.

We analyse the energy in the delay-line output vector after one complete pass through the feedback matrix.

Step 1: Energy after the mixing step. At each time step, the input to the delay lines is (lossless, no input). The energy of this new injection vector is

If , then

The mixing step preserves the total squared norm — no energy is created or destroyed.

Step 2: Energy in the delay lines. Each delay line is a pure shift register: the sample injected at time re-emerges steps later. The delay operations do not modify sample values, only their timing. Therefore the total energy stored across all delay lines evolves as

(the energy leaving the delay-line outputs equals the energy injected back in). The system is lossless: for all with zero input.

Step 3: Stability with loss filters. When , each feedback path multiplies its channel by , so . Combined with Step 1:

Energy strictly decreases at every step (for nonzero state), so the system is asymptotically stable and the reverb tail decays to silence.

Theorem 38.2 (Hadamard Matrix Orthogonality)

The normalised Hadamard matrix (where is the unnormalised version with entries in ) satisfies

That is, is an orthogonal matrix and hence a valid lossless FDN feedback matrix. The matrix-vector product requires no scalar multiplications: only additions and subtractions via the Fast Hadamard Transform.

Proof.

Step 1: Row orthogonality of . We show by induction on (where ).

Base case : , and

Inductive step. Suppose . Write

Then

By induction, for all .

Step 2: Normalisation. Define . Then

Since is symmetric ( for the standard Walsh-Hadamard ordering), the same argument gives .

Step 3: No multiplications. Each entry of is . The product is a sum of terms , requiring only additions and subtractions. The recursive structure enables the Fast Hadamard Transform: split , then , reducing the -point transform to two -point transforms plus additions — a recursion with total cost additions.

Theorem 38.3 (Coprime Delays and Uniform Mode Density)

Let an FDN have delay lines with lengths (in samples). The Z-domain poles of the lossless system ( orthogonal, no loss filters) are located on the unit circle at frequencies determined by the delay lengths and the mixing matrix. In the scalar case (a single feedback loop of length ), the poles are the equally spaced roots of unity , .

For with coprime delay lengths ( for all ), the total number of distinct pole frequencies is maximised: no two delay lines share a resonance frequency, so the modal density in the frequency domain is as uniform as possible given the total modal count . When delay lengths share a common factor , poles cluster at multiples of , creating spectral colouration audible as metallic resonances.

Proof.

Scalar case. A single feedback loop of delay with gain has transfer function . For the lossless case , the poles are the solutions of , giving equally spaced points on the unit circle at angles , .

Common factor case. Suppose two delay lines have lengths and with . The poles of line 1 occur at angles and those of line 2 at . Every angle that is a multiple of appears in both sets: there are shared resonance frequencies. At these frequencies, both delay lines reinforce each other, creating amplitude peaks times larger than isolated resonances. In terms of mode density, the distinct pole frequencies are fewer than the one would obtain from coprime lengths.

Coprime case. If for all , then by the Chinese Remainder Theorem the resonance grids and share only the trivial intersection at . The total number of distinct non-trivial resonance frequencies is , the maximum achievable. With no clustering, the inter-mode spacing is approximately uniform at Hz, and no single frequency is reinforced by more than one delay line. The resulting impulse response grows in echo density without spectral colouration.

Prop 38.1 (RT60 Calibration Formula)

For an FDN with sample rate and target RT60 of seconds, the per-sample loss factor for delay line is

With this choice, the energy envelope of the -th delay line decays by a factor of (60 dB) after exactly seconds, regardless of the delay length .

Proof.

After samples, the signal in delay line has made full traversals of the loop. Each traversal multiplies the amplitude by , so the total amplitude after seconds is

Setting this equal to (a 60 dB amplitude decay corresponds to a factor of ):

Wait — the right-hand side does not depend on at this stage. However, the formula is conventionally written with the factor to express the per-sample gain as a function of delay length: for a longer delay line, more attenuation per sample is required to achieve the same overall decay time. Taking logarithms of :

This single formula gives the same for all delay lines, ensuring that each loop decays 60 dB in seconds — the RT60 is uniform across all channels. (The form with explicit in the exponent arises when one writes the loss per delay-sample rather than per time-sample, which is common in the DSP literature.)


Numerical Examples

Example 1: A minimal 4-channel FDN with .

Choose , sample rate Hz, and coprime delay lengths in samples:

Channel (samples) (ms)
1 1499 31.2
2 1777 37.0
3 2311 48.1
4 3001 62.5

These lengths are all prime, hence mutually coprime: , etc. Total modal count: distinct resonance frequencies, roughly one mode every Hz.

The feedback matrix is the normalised :

Verification: each row has norm . Inner product of rows 1 and 2: . All off-diagonal inner products are zero by the same cancellation argument. Hence and the FDN is lossless.

Example 2: RT60 calibration.

For the FDN above with target s:

Each sample is attenuated by a factor of . After samples (2.5 s), the amplitude is , a 60 dB decay.

Per-channel verification for delay line 1 ( samples): the loop is traversed times. Each traversal multiplies amplitude by . After 80.1 traversals: . The 60 dB condition is satisfied for all delay lines simultaneously.

Example 3: Schroeder’s metallic flutter — a failure case.

Take delay lines with lengths 150, 300, 450 samples. Then , , . All three delay lines share resonances at multiples of Hz. These 150 shared frequencies are reinforced three times over, producing pronounced spectral peaks spaced 320 Hz apart — audible as a metallic, pitched colouration of the reverb tail. Replacing the lengths with 149, 211, 307 (all prime) eliminates the shared resonances entirely.


Musical Connection

音乐联系

The Hadamard matrix as a concert hall.

A physical concert hall creates reverberation through thousands of reflections off walls, seats, balconies, and audience members. Each reflection path is a delay line; the angles and materials determine how energy is redistributed across directions after each reflection. The subjective quality of a great concert hall — the sense of being enveloped in sound from all directions, with a smooth decay lasting 1.8–2.2 seconds — emerges from the combined action of all these paths.

The FDN is a mathematical distillation of this process. The delay lines represent “reflection paths,” the Hadamard matrix represents the geometry of the hall (how energy from one path scatters into all others), and the loss filters represent the material absorption that causes the decay. What makes the Hadamard matrix musically apt is exactly its mathematical property: it is maximally mixing in the sense that every output depends equally on every input — no “acoustic island” exists where energy can pool without spreading to other delay lines. This mathematical property corresponds to the perceptual quality of spatial diffuseness: in a well-designed hall, the reverberant field should feel omnidirectional, with energy arriving from all directions with equal probability.

The number-theory connection to EP29.

Both this episode and

EP29

use number theory to control acoustic phenomena. In EP29, the slot depths of a QRD diffuser are : a quadratic residue sequence chosen so that the DFT of the resulting phase grating has flat magnitude, scattering energy equally in all spatial directions. In this episode, coprime delay lengths are chosen so that the resonance grids of different delay lines do not overlap, distributing modes uniformly across frequency. Both are instances of the same meta-principle: algebraic structure in the design parameters maps to uniformity in the physical response.

The Hadamard matrix itself echoes the group-theoretic structure of

EP04

the rows of form a group under componentwise multiplication (each entry is , and the product of two rows is again a row of the matrix). This algebraic closure is why the recursive construction works: the matrix structure is self-similar across scales, just as the all-interval row in has a regularity inherited from the cyclic group structure.

What the four knobs mean.

Modern reverb plugins (Valhalla Room, FabFilter Pro-R, Eventide BlackHole) expose four conceptual parameters, each corresponding to a mathematical degree of freedom in the FDN:

  • (number of delay lines): controls the density and smoothness of the early echo pattern; larger gives more diffuse, less grainy reverb at greater computational cost
  • (delay lengths): controls the pitch of the reverb tail (longer delays = lower resonances) and, through coprimality, the absence of metallic colouration
  • (feedback matrix): controls how quickly energy spreads across delay lines; a Hadamard matrix achieves full spread in one step, while a sparse or diagonal matrix produces a drier, more echo-like effect
  • (loss filters): controls RT60 and frequency-dependent decay; making high-frequency losses larger than low-frequency losses mimics air absorption in large spaces

The “character” of a reverb — plate, spring, room, hall, chamber — is a specific point in this four-dimensional parameter space. Plate reverb (originally an electromechanical device: a suspended steel sheet driven by a transducer) has a very short RT60, dense modal structure, and flat frequency response; its FDN equivalent uses short, densely-packed delay lines with nearly flat loss filters. A cathedral reverb has RT60 > 5 s, emphasises low frequencies, and requires large to achieve the required mode density. The mathematics is the same in both cases.


Limits and Open Questions

  1. Stability is sufficient, not necessary. Theorem 38.1 proves that orthogonality of is sufficient for lossless stability. It is not necessary: non-orthogonal matrices (e.g., upper triangular matrices with unit diagonal) can also yield stable FDNs, and some of these produce perceptually different diffusion characteristics. The space of stable feedback matrices is larger than the set of orthogonal matrices; finding the perceptually optimal stable matrix for a given application is an open design problem.

  2. Coprime delays as a heuristic. Theorem 38.3 is stated in the scalar () case and extended informally to the multivariate case. For a general unitary feedback matrix, the pole locations are not simply the union of the scalar pole sets — they depend on the eigenstructure of interacting with all delay lengths simultaneously. A rigorous necessary-and-sufficient condition for uniform mode density in the full matrix case remains an open problem in digital audio signal processing.

  3. Perceptual metric. The mathematical metric for reverb quality used here (mode density, spectral flatness of ) is a proxy for perceptual quality. The psychoacoustic attributes that listeners care about — envelopment, intimacy, warmth, clarity — are multidimensional and do not map cleanly onto any single mathematical quantity. The relationship between FDN parameters and ISO 3382 room acoustic descriptors (such as the lateral energy fraction and the interaural cross-correlation ) is an active research area.

  4. Time-varying FDNs. Static FDNs with fixed delay lengths produce a stationary impulse response. Slightly time-varying delay lengths (modulated by a low-frequency oscillator) are used in practice to break up periodic artefacts and add a subtle “living” quality to the reverb tail. The stability analysis of Theorem 38.1 does not directly extend to time-varying systems; proving stability of slowly-modulated FDNs requires Floquet theory or Lyapunov methods.

  5. High-order Hadamard matrices. The recursive construction of Definition 38.5 yields only for . Hadamard matrices of other orders (e.g., ) are known to exist (the Hadamard conjecture asserts they exist for all ), but the conjecture remains unproved in general. For FDN design, non-power-of-2 sizes would offer finer control over the balance between mixing quality and computational cost.

Conjecture (Perceptually Optimal FDN Dimension)
For a given target RT60 and room type (categorised by Schroeder frequency and modal density), there exists an optimal FDN order that minimises a combined perceptual cost function measuring metallic colouration and insufficient diffusion, subject to a real-time computational budget of multiply-accumulate operations per sample. The conjecture is that grows logarithmically rather than linearly with the required modal density — i.e., doubling the target mode count requires increasing by roughly 1, not 2. Falsification criteria: implement FDNs of orders with coprime prime delays scaled to equal total delay, measure the perceptual diffusion score via MUSHRA listening tests for three room types; if the improvement from to is not significantly larger than from to , the logarithmic-growth conjecture is falsified.

Academic References

  1. Jot, J.-M., & Chaigne, A. (1991). Digital delay networks for designing artificial reverberators. Proceedings of the 90th Audio Engineering Society Convention, Paris. (Original FDN paper; introduces the unitary feedback matrix stability condition.)

  2. Schroeder, M. R. (1962). Natural sounding artificial reverberation. Journal of the Audio Engineering Society, 10(3), 219–223. (First practical artificial reverb algorithm using comb and allpass filters.)

  3. Smith, J. O. (2010). Physical Audio Signal Processing. W3K Publishing. Available at ccrma.stanford.edu/~jos/pasp/. Ch. 3 (comb filters, allpass filters, Schroeder reverb) and Ch. 4 (FDN theory, stability, RT60 calibration).

  4. Stautner, J., & Puckette, M. (1982). Designing multi-channel reverberators. Computer Music Journal, 6(1), 52–65. (Early work on multi-channel delay network reverb, predecessor to the formal FDN.)

  5. Zölzer, U. (Ed.). (2002). DAFX: Digital Audio Effects. Wiley. Ch. 7 (reverberation algorithms including FDN and Hadamard matrix application).

  6. Sylvester, J. J. (1867). Thoughts on inverse orthogonal matrices, simultaneous sign successions, and tessellated pavements in two or more colours. The London, Edinburgh, and Dublin Philosophical Magazine, 34, 461–475. (First construction of what are now called Hadamard matrices.)

  7. Hadamard, J. (1893). Résolution d’une question relative aux déterminants. Bulletin des Sciences Mathématiques, 17, 240–246. (The Hadamard determinant bound; gives the name to the matrix family.)

  8. Poletti, M. A. (1996). Optimised reverberators for use in auralization systems. Proceedings of the Institute of Acoustics, 18(8). (Optimisation of FDN parameters including delay-length coprimality for perceptual quality.)

  9. Schlecht, S. J., & Habets, E. A. P. (2017). On lossless feedback delay networks. IEEE Transactions on Signal Processing, 65(6), 1554–1564. (Rigorous modern treatment of lossless FDN stability; generalises Theorem 38.1 to non-orthogonal matrices.)

  10. Rocchesso, D., & Smith, J. O. (1997). Circulant and elliptic feedback delay networks for artificial reverberation. IEEE Transactions on Speech and Audio Processing, 5(1), 51–63. (Circulant feedback matrices, relationship to DFT, computational efficiency.)

  11. Välimäki, V., Parker, J. D., Savioja, L., Smith, J. O., & Abel, J. S. (2012). Fifty years of artificial reverberation. IEEE Transactions on Audio, Speech, and Language Processing, 20(5), 1421–1448. (Comprehensive survey from Schroeder 1962 through modern FDN variants; excellent historical and technical overview.)

  12. Gardner, W. G. (1998). Reverberation Algorithms. MIT Media Lab Technical Report. (Practical implementation guide for FDN reverb including Hadamard matrix coding and RT60 calibration.)