Building a morphing reverb

From what a room does to sound, through delay lines and feedback delay networks, to the moving mixing matrix SD-Reverb runs and its Airwindows roots.

44 minutes read

9252 words

A reverb turns one sound into thousands of quieter copies of itself, spread over time. A digital reverb builds those copies from delay lines wired into a feedback loop. The table of numbers that wires them decides how the copies pile up. SD-Reverb is an audio plugin whose table does not hold still. It travels a closed path over seconds to minutes, fifty seconds at the default setting. The resonances of the room therefore drift while its decay time stays where the control put it. The architecture and the delay lengths come from the Airwindows reverbs of Chris Johnson. Moving the table is a published technique, reached here by its own route.

Everything below builds from one physical fact about rooms up to the six numbers SD-Reverb recomputes every 64 samples. No prior signal processing is assumed. A plugin is an effect loaded into a host program such as a recording application. It receives audio from that host in blocks of a few hundred samples at a time. SD-Reverb’s source is not published. The file paths below name where each mechanism is implemented, and the few short extracts are quoted from those files.

Reverberation in a room

Sound leaves a source in every direction. One path reaches the listener without touching anything, and it arrives first and loudest. Every other path reflects off a surface before arriving, so it travels further, arrives later, and loses energy at each reflection.

The first reflections arrive as separate events, spaced by the distance between the surfaces that produced them. Later arrivals have bounced many times and there are far more of them, because every reflection spawns further reflections. Past a few tens of milliseconds they overlap into a decaying wash that no listener can resolve into separate echoes. That wash is the tail.

The response of a room to one short click is called its impulse response . It records everything the room does to any sound, because any sound is a sum of clicks and the room treats each one the same way.

Density is the count of separate arrivals per second. It starts low in the early reflections and rises as the bounces multiply. A reverb that is too sparse is heard as a stutter of distinct echoes. Building density fast enough is half of what the structures below are for.

A decibel , written dB, is a ratio on a logarithmic scale. A drop of 20 dB is a drop to one tenth of the amplitude, and a drop of 60 dB is a drop to one thousandth. Reverberation time measures how long the tail lasts, and $T_{60}$ is the time the sound takes to fall by 60 dB after the source stops. SD-Reverb’s Decay control sets $T_{60}$ from 0.2 to 50 seconds on a logarithmic curve, with a default of 2.4 seconds.

Kinds of digital reverb

What goes between the input and the output is the whole design question, and the answers fall into families.

Convolution takes a recorded impulse response and replays it for every moment of the input, adding the results. What comes out is that room, heard from the position the recording was made at, as faithfully as the recording allows. The room cannot then be altered, only swapped for another recording or cut short. The arithmetic is heavy. A response of two seconds is tens of thousands of numbers long, and every one of them is used at every moment.

A delay network builds a structure whose own response has the right shape instead. Nothing is measured and no particular room is reproduced. What is aimed at is the density, the decay and the colour. It costs a fraction of the convolution, and each of those can be moved while the audio runs. Schroeder’s combs and allpasses are of this kind, and so are the feedback delay networks this page is about.

Between the two sit the physically motivated networks. A scattering delay network puts a node at each surface of a modelled room. The delay lines between the nodes carry the propagation, so the geometry is a parameter rather than a recording. It costs more than a plain delay network and much less than simulating the wave field across the whole room.¹

SD-Reverb is a delay network. The rest of this page is how one is built.

Samples and delay lines

An audio signal inside a computer is a list of numbers taken at a fixed rate. At a sample rate $f_s$ of 48 kHz there are 48000 numbers per second per channel. The number at position $n$ is written $x[n]$ , and the output the effect produces at that position is $y[n]$ .

A delay line stores the recent past of a signal. It is an array plus a write position. Each sample the input is written at that position and the position advances by one. A read taken $M$ places behind the write position returns what came in $M$ samples ago.

When the write position runs off the end of the array it wraps to the front, which is why the structure is called a circular buffer. SD-Reverb rounds every buffer up to a power of two, so the wrap costs one bitwise operation and no division. A delay of 961 samples at 48 kHz is 20.02 ms.

float read()  const noexcept { return buffer_[(writeIdx_ - delay_) & mask_]; }
void  write(float x) noexcept { buffer_[writeIdx_ & mask_] = x; }
void  advance()      noexcept { writeIdx_ = (writeIdx_ + 1) & mask_; }

That is the whole of src/dsp/DelayLine.hpp that the audio path uses. The mask is the buffer’s length minus one, and the AND is what a power of two buys.

A read from a delay line is called a tap. Taps in the feedback loop are whole numbers of samples. The pre-delay line is the exception, because a control that glides has to read between two stored samples, and it reads with four-point cubic Hermite interpolation (src/dsp/InterpolatedDelayLine.hpp). Feedback taps are kept whole so that changing a length cannot bend the pitch of the tail.

Feedback comb filters

One delay line with its output added back to its own input repeats the sound. The fraction sent back is the feedback gain $g$ .

Feed a single click into equation (2). The output is a click at position 0, another at $M$ , another at $2M$ , each one $g$ times the height of the one before. At $g = 0.9$ and $M = 961$ the repeats are 20.02 ms apart and the result is an audible flutter.

Each round trip multiplies the amplitude by $g$ . Sixty decibels down means the amplitude has reached $10^{-3}$ , which takes $-3/\log_{10} g$ round trips, and each round trip takes $M/f_s$ seconds. Multiplying those two gives $T_{60} = -3M/(f_s \log_{10} g)$ . Solving for the gain gives the line every reverberator needs.

Equation (3) is the sixty-decibel time read backwards, and it is what src/dsp/MultibandDamping.cpp evaluates to turn a decay time in seconds into a per-line gain.

Repeating at a fixed spacing also shapes the frequency response, which is the gain the structure applies at each frequency. Feed the loop a steady tone. The copy that returns after $M$ samples is late by $M/f_s$ seconds. When that lateness is a whole number of the tone’s periods the copy comes back in step with the input. The two add and the loop builds up. That happens whenever $fM/f_s$ is a whole number, which is at every multiple of $f_s/M$ . Between those frequencies the returning copy is out of step and partly cancels.

A frequency the structure builds up is a resonance. For $M = 961$ at 48 kHz the resonances stand 49.9 Hz apart, and the response takes the shape the filter is named for.

Six lines of different lengths put their resonances in six different places, which is the lower panel of the figure. Two whole numbers are coprime when no whole number above 1 divides both of them. Lengths that are not coprime stack resonances on top of each other and reinforce those frequencies. SD-Reverb’s six late lengths, as tabled for 48 kHz at the reference Size, are coprime in every pair. The first coincidence is therefore pushed as far out as the arithmetic allows. The Size control scales and rounds those lengths, and rounding does not preserve the property. At Size 0 the Hall’s second and fourth lines are 630 and 1552 samples, both even.

Allpass sections

A filter changes the level of each frequency it passes, and it changes the timing. The comb section measured levels and said nothing about timing. An allpass section leaves every level alone and moves only the timing. The signal $w$ below is internal to the section.

Both halves of equation (4) are doing work. The first, $w[n] = x[n] + g\,w[n-M]$ , is the feedback comb of equation (2), and on its own it would ring at every multiple of $f_s/M$ . The second subtracts a scaled copy of the same signal, and the two ripples cancel exactly. What survives the cancellation is the delay each frequency picked up on the way through, and that delay is different at every frequency.

Showing that the level is unchanged takes one piece of notation. Write $z^{-1}$ for a delay of one sample, so $z^{-M}$ is a delay of $M$ samples. Equation (4) then rearranges into a single ratio of output to input, which is called the transfer function.

Substituting $z = e^{j\omega}$ reads that ratio at one frequency. Here $\omega$ is $2\pi$ times the frequency divided by $f_s$ , and $j$ is the square root of minus one. Both are built up under Mathematics, complex numbers first and then the form $e^{j\omega}$ under exponential form . Those pages write the square root of minus one as $i$ .

Taking $e^{-j\omega M}$ out of the numerator of equation (5) leaves $1 - g e^{j\omega M}$ behind it. That is the complex conjugate of the denominator, meaning the same number with the sign of its imaginary part flipped. A number and its conjugate have equal distance from the origin, and $e^{-j\omega M}$ has distance 1, so $|H| = 1$ at every frequency.

One section is seven lines, and both terms of equation (4) are visible in it.

const float d = s.buf[(s.w - s.delay) & s.mask];   // w[n-M]
const float v = x + g * d;                         // the feedback half
s.buf[s.w & s.mask] = v;
s.w = (s.w + 1) & s.mask;
return d - g * v;                                  // the feed-forward half

A click is every frequency arriving at the same instant. After a section they no longer arrive together, which is the whole of what the structure does.

One section turns the click into 25 arrivals and drops the tallest from 1 to 0.70. Eight in series turn it into 1612 within the same window, and the tallest is 0.084. Chaining multiplies the count, because each arrival leaving one section becomes a burst in the next.

None of that is energy removed. Measured over 400 ms the chain has produced 10911 arrivals and its total energy is 1.0000, the same as the click that went in. That is $|H| = 1$ read in the time domain, and it is why the chain cannot be dropped from the path. It destroys the transient and leaves the spectrum where it was.

SD-Reverb chains eight sections per channel, with delays from 1.7 to 9.3 ms on the left and 1.9 to 9.7 ms on the right, at a gain of 0.70 (src/dsp/AllpassDiffuser.hpp). The chain is placed between the early stage and the late field, both described below.

Comb and allpass reverberators

Schroeder’s 1962 design put a parallel bank of feedback combs in front of a series chain of allpass sections. The combs supply the decay, the allpasses raise the density, and the whole thing is cheap enough to have run on the computers of the day. Julius O. Smith’s account gives the working numbers. Three or four parallel combs are needed for an adequate frequency response, and five series allpasses give roughly 810 echoes per second.²

Two properties bound that structure. The echo count from a comb bank grows in proportion to the number of combs, so density is bought one delay line at a time. And each comb repeats at one fixed interval forever, which colours the tail with its own resonances for as long as the tail lasts.

Later designs worked on both. Moorer put a lowpass filter inside each comb’s feedback path. A lowpass passes frequencies below a corner frequency and reduces those above it. High frequencies then decay faster than low ones, which is what air and soft surfaces do. Dattorro arranged allpass sections inside a single feedback loop, giving a dense tail from a small number of delay lines. The feedback delay network is the structure SD-Reverb uses.

Feedback delay networks

The construction generalises the comb bank. $N$ delay lines run in parallel, and each one receives a weighted sum of what all of them produced. A single input therefore reaches every line within one pass. The count of distinct paths multiplies by $N$ at each pass afterwards.

Gerzon proposed the structure for reverberation in 1971. His observation was that one feedback comb sounds poor while several cross-coupled ones sound good. Stautner and Puckette arrived at it again in 1982, with a four-channel version and the conditions under which it stays stable. Jot and Chaigne made a design method of it in 1991, by scaling each line’s damping to its own length. That is what lets the decay time be set per frequency band rather than per line.³

Building one means choosing the delay lengths, the matrix and the gains, and each governs something different. The lengths decide where the resonances fall and how sparse the first arrivals are. The matrix decides how fast energy spreads between the lines. The gains decide how long it takes to leave. The sections below take them in that order.

Write $y_i[n]$ for the output of line $i$ , $M_i$ for its length and $g_i$ for its gain. The value written into that line is $d_i[n]$ . Its share of the input is $b_i$ and its share of the output is $c_i$ .

The $\sum$ is shorthand for adding up a run of terms, and nothing in equation (6) is worse than long addition once it is written out. Take the first line of six. What it hands back is what it stored $M_1$ samples ago. What goes into it is its share of the input, plus every line’s output scaled by that line’s gain and by one weight each.

The other five lines are the same with their own row of weights. The output is one more sum of that shape. It reads $z[n] = c_1 y_1[n] + \dots + c_6 y_6[n]$ .

A matrix is a rectangular table holding all those weights at once, one row per line. Multiplying a list of numbers by it produces a new list, and each entry of the new list is a weighted sum of the old ones. The numbers $a_{ip}$ form the mixing matrix $A$ , and the choice of $A$ decides whether the network is stable.

The energy in the six lines is the sum of their squares. That is how the length of a line on paper comes from the squares of its two sides. Call the sum the length of the list, written $\|\mathbf{y}\|$ . Mixing is only safe if it leaves that total alone.

Orthogonal matrices are the ones that do. The transpose $A^{\mathsf{T}}$ is the same table with rows and columns swapped. The matrix $I$ has ones down the diagonal and zeros elsewhere, so it leaves any list exactly as it found it. A matrix is orthogonal when $A^{\mathsf{T}}A = I$ , which says that swapping the table round and multiplying by the original undoes it.

Read left to right, each step is small. The energy after mixing is the mixed list multiplied into itself. Multiplying a pair like that lets the two $A$ terms come together as $A^{\mathsf{T}}A$ . That product is $I$ by the condition just given. What is left is the list multiplied into itself, which is the energy before mixing. A worked instance with six actual numbers is at the end of the next section.

A length-preserving matrix moves energy between lines without creating or destroying any. Decay is then entirely the work of the gains $g_i$ , and every $g_i < 1$ makes the loop shrink whatever passes round it. Setting each $g_i$ from equation (3) with that line’s own length gives every line the same $T_{60}$ . That is Jot’s construction, and it is why such a network can hold one decay time across six lines of unequal length.⁴

Any length-preserving matrix serves, and the choice is a real one. A Hadamard matrix mixes every line into every other at equal weight, spreading energy fastest. A circulant matrix gives each line the same pattern of neighbours. The Householder form of the next section is the one both reverbs on this page use. It can be applied without ever being built, which is the reason.

Collecting the weights of equation (6) gives two more matrices. The input shares $b_i$ form $B$ , one row per line. The output shares $c_i$ form $C$ , one row per output channel.

Householder mixing matrices

The mixing matrix $A$ of equation (6) is the one thing in the network that connects the delay lines to each other. It is $N$ by $N$ , so 6 by 6 in the late field here and 3 by 3 in the early stage. On every sample it takes the $N$ values just read from the lines and returns $N$ new ones, each a weighted sum of all of them. That is a matrix times a vector , with the matrix holding the weights. Without it, six delay lines are six separate echoes.

Equation (8) already said what $A$ has to be, which is length-preserving. That leaves a large family to choose from, and the Householder form is one member of it.

The member SD-Reverb starts from needs no matrix at all to describe. Add the six numbers together. Subtract a third of that total from each of them. That is the whole mixing step.

Try it on one line carrying everything and the other five silent, so the six numbers are 1, 0, 0, 0, 0, 0. The total is 1, so a third comes off each. What leaves is $2/3$ on the line that had the signal and $-1/3$ on each of the other five.

Energy is the sum of the squares, and here that is $4/9 + 5 \times 1/9$ , which is 1 again. The signal has been spread across all six lines and none of it has been gained or lost. That is what equation (8) claimed in general, done in arithmetic.

Written as a table, that subtraction is a 6 by 6 grid. It holds $2/3$ down the diagonal and $-1/3$ in every other cell, which is the left panel above. Multiplying by it is the same arithmetic in a shorter notation. Each output is one row of the grid, multiplied entry by entry into the six inputs and summed. The Airwindows code described in the section after next computes $2x_i$ minus the sum of the other five, which is that same grid scaled by three.

That grid is one member of a family, and the family is written with a single vector. A Householder matrix is built from one vector $v$ of length 1, called a unit vector. The outer product $vv^{\mathsf{T}}$ is the grid whose cell in row $i$ and column $j$ holds $v_i v_j$ .

It is symmetric, so $A^{\mathsf{T}} = A$ , and multiplying it by itself gives

because $v^{\mathsf{T}}v = \|v\|^2 = 1$ . Every unit vector therefore yields an exactly length-preserving matrix. Geometrically $A$ reflects a list of $N$ numbers through the flat surface at right angles to $v$ . A reflection turns a direction round and leaves the length alone, which is what equation (8) asked for. Energy moves between the lines and none of it is made or lost.

Carrying $v$ instead of the matrix is what the rest of the page rests on. The matrix never has to be stored or multiplied out. The dot product $v \cdot x$ multiplies matching entries of two lists and adds the results, and applying equation (9) to a list $x$ is

which is one dot product and one scaled subtraction, and it is the regrouping of applying a reflection .

float dot = 0.0f;
for (std::size_t i = 0; i < kN; ++i) dot += v_[i] * x[i];
const float d = 2.0f * dot;
for (std::size_t i = 0; i < kN; ++i) x[i] -= v_[i] * d;

At $N = 6$ that is 12 multiplications where a general matrix would need 36. Six numbers hold the whole thing, and the matrix is never built.

Any unit vector at all gives a length-preserving matrix, so $v$ can be moved anywhere on the unit sphere and the result is still one. That is what makes the trajectory of a later section possible, and it is the reason this form was chosen over the rest of the family.

The grid the section opened with is this family’s even member. Every entry of $v$ equal to $1/\sqrt{N}$ gives $A = I - (2/N)J$ , with $J$ the grid of all ones. At $N = 6$ that is the $2/3$ and $-1/3$ already seen. It is the most even mixing a Householder matrix can do, and it is the matrix the Airwindows reverbs use.

Airwindows reverbs

Chris Johnson has published a decade of reverb research as MIT-licensed source at Airwindows. The individual plugin sources for VerbSixes, kRockstar, kCyberCity, kCosmos and kCathedral5 show one architecture with variations.⁵

The mixing matrix is the even Householder, written out by hand. Every stage of the 6 by 6 field computes 2*x_i minus the sum of the other five taps, which is $3(I - (2/6)J)$ , the even Householder times three. The 3 by 3 stages compute x_i minus twice the sum of the other two, which is $3(I - (2/3)J)$ on three entries. The vector $v$ never appears, because with a fixed even $v$ the sums can be written out directly.

The late field is a cascade. Six stages sit in series, each holding six delay lines. The field therefore has 36 delay lines per channel and applies the same matrix six times per pass. The output of the sixth stage returns to the first. Stereo is cross-coupled, with the left bank’s final stage producing the right channel’s feedback vector.

Four further properties matter for anything derived from that work.

The reverb core runs at an integer subdivision of the host sample rate, and one step of that slower clock is called a tick. Input samples are averaged between ticks, and one reverb iteration runs on the average. The sparse output stream is interpolated back up with a quadratic Bezier, a curve fitted through the three most recent tick outputs. In kCosmos the Bezier is replaced by a windowed-sinc filter, a standard reconstruction filter, built each block. The delay integers are therefore counted in ticks, which is why they are small.

Damping is done with running averages and nothing more elaborate. The kRockstar and kCyberCity fields apply a two-sample running average to the first feedback tap. VerbSixes and kCathedral5 apply a three-to-one weighted average of the same tap with its own past value, taken without its sign. In kRockstar and kCyberCity a cascade of running averages on the output doubles as the reconstruction filter for the undersampling and as the top-end damping.

Room tone is the technique Chris Johnson calls alive silence. The kRockstar and kCyberCity sources run two bounded random walks, each a value nudged up or down by a small random step every tick. The first is added to the tap-A feedback of each channel, so the noise circulates through the whole field and takes on the room’s colour. The second is the regeneration coefficient itself, the fraction of the loop signal fed back round. It wanders slowly around the decay setting, so the decay time breathes.

The feedback path contains no waveshaper, meaning no curve that each sample is mapped through. What kRockstar calls the nonlinearity is a level-dependent decay. Its applyNonlin term rises from 0 at silence towards 1 as the input approaches full scale. The only use of that term is to lift the ceiling of the wandering regeneration by a fraction of a percent. The kCosmos sources run the opposite sign, subtracting an energy measure from the regeneration so that a loud loop damps itself.

The delay tables were found by search, with years of genetic-algorithm compute behind them. A genetic algorithm scores many candidate sets, keeps the best, and recombines them into the next generation. The header comments name the room each set was aimed at, such as kRockstar’s “2 to 193 ms, 1137 seat hall”. Those millisecond ranges are path sums through the cascade, not single tap times. For kCathedral5 at full rate the six per-stage maxima sum to 11439 ticks, which at 44.1 kHz is 259 ms and matches the header exactly.

SD-Reverb signal chain

SD-Reverb keeps the two-stage shape, a 3 by 3 early field feeding a 6 by 6 late field, and runs both at the host sample rate with one stage each. The dry signal is the input untouched, and the wet signal is what the reverb produces from it.

A block arrives and is processed in chunks of 64 samples. The dry signal is copied first. The wet path then takes a highpass set by Low Cut, which passes frequencies above its corner and reduces those below. That control runs from 20 to 500 Hz and defaults to 80 Hz. The pre-delay line follows, adjustable from 0 to 200 ms and defaulting to 20 ms. An envelope follower reads the signal at that point and nowhere else, tracking the size of the signal by rising quickly and falling slowly.

The early stage is three delay lines of 6 to 26 ms at the default Size, at a fixed feedback of 0.4. It produces a handful of separate reflections and no tail. One copy of its output is tapped straight to the wet mix under the Early Level control. Another copy passes through the allpass chain and drives the late field. Between the two ends of the Diffusion control the late field receives a blend of the chain’s output with the signal that bypassed it. That blend is not magnitude-flat, although either end of it is. At a setting $\delta$ its response lies between $|1 - 2\delta|$ and 1 across frequency.

The late field runs seven steps per sample (src/dsp/FDN.cpp).

Reading the six lines comes first. The output taps of $C$ are taken second, straight off those values, before damping or mixing touches them. Damping and the decay gains come third, the saturation fourth, the room tone fifth, the mixing matrix sixth. Seventh, the input is injected through $B$ and the six results are written back.

The order of the first two steps is what keeps the wet signal brighter than the signal circulating in the loop. Every pass through the loop damps the highs further, and the output tap reads before that pass.

Wet and dry are combined with an equal-power law that holds the total power constant as the control moves. The dry gain is $\cos(\tfrac{\pi}{2}\mu)$ and the wet gain is $\sin(\tfrac{\pi}{2}\mu)$ for a Mix setting $\mu$ between 0 and 1. Width acts before that sum, on the mid and side components of the wet signal. The mid is half the sum of the two channels and the side is half their difference.

Mode and Size change the delay lengths, which cannot be done to a running network without a click. SD-Reverb keeps two complete engines for that reason. A change configures the standby engine and copies the sounding engine’s delay buffers into it. The same audio is then re-tapped at the new lengths, and the two engines crossfade over 40 ms, one gain rising while the other falls. When the wet signal is already below $-100$ dBFS there is nothing to carry across, and the change is made in place at once. Size scales every length by $0.5 + \text{Size}$ , and a change smaller than 0.03 in that factor does not trigger a rebuild.

The lengths are counted in samples, as tabled for 48 kHz below, and the same counts are used at every host rate. At 96 kHz each delay therefore lasts half as long, the whole echo pattern compresses by two and the room reads as smaller. The decay time is unaffected, because equation (3) uses the rate.

Density from six lines

The Airwindows late field holds 36 delay lines and this one holds six. Two measurements bear on what that costs, and they point opposite ways.

The first is modal density, the count of resonances per hertz. A delay of $M$ samples contributes $M$ resonances, so a network’s total is the sum of its lengths over the sample rate. That is the same number as its total delay in seconds.

Six times the lines gives 1.8 times the resonances. Each of his is short, because one line is a fraction of a circuit through the cascade, and because they are counted in undersampled ticks. The row above reads his in the 24 kHz ticks his core runs at when it undersamples by two.

The second measurement is echo density, the count of separate arrivals per second, and here the cascade is far ahead. Paths multiply by six at every stage. One circuit through six stages therefore offers up to $6^6 = 46656$ of them, where one stage offers six.

That gap is what the allpass chain closes, and it closes it before the late field sees the signal.

His cascade arrives first. It passes 21000 arrivals per second by 15 ms and reaches one per tick by 25 ms, which is what 36 lines in series buy. Six lines on their own take a quarter of a second to reach 6500. The same six behind the eight allpass sections are at 5600 within 30 ms and at 47800 by 50 ms. The chain costs nothing in the spectrum, because each section has magnitude 1 at every frequency.

The two ceilings differ because the two fields run at different rates. His core ticks at 24 kHz, so it cannot place more than 24000 arrivals in a second. This one runs at 48 kHz and saturates at twice that. The cascade gets there sooner and stops lower.

Six lines and eight allpass sections also cost less than 36 lines. The table below counts the structures rather than timing them, over the late field and the allpass chain. The early stage and the pre-delay are outside the count on both sides.

His 36 lines are per channel and tick at 24 kHz, which is what makes them affordable. The six here are shared by both channels at 48 kHz, and the 16 allpass sections are two channels of eight. Damping, room tone and his reconstruction filter sit outside the count, because they are not what is being compared.

Two mechanisms carry what is left. The mixing matrix moves, so a sparse set of resonances drifts rather than standing at fixed frequencies, and the next section is the path it takes. Room tone circulates under the tail and fills the gaps between resonances.

The lower order is a choice as much as a cost, and the design records it as a bet. Sparser mixing scrambles phase less, and that is what it credits for keeping the high end present.

Matrix trajectory

The whole state of a Householder matrix is one unit vector. The set of matrices reachable at $N = 6$ is therefore the set of unit vectors with six entries. That set is the unit sphere in six dimensions, written $S^5$ . Any point on it gives an exactly length-preserving matrix by equation (10). A path across the sphere is a path through length-preserving matrices. No correction step is needed, and no midpoint on the path stops conserving energy.

Movement along the sphere uses spherical linear interpolation , written slerp. A straight line between two unit vectors passes inside the sphere, so its midpoint is shorter than the radius. The matrix there would not preserve length. Normalising the line puts the point back on the surface, but it then moves unevenly. For two vectors 60 degrees apart, the first fifth of the path covers 10.9 degrees and the middle fifth covers 13.2. Slerp holds one angular speed the whole way, along the shortest arc. For unit vectors $a$ and $b$ separated by an angle $\theta = \arccos(a \cdot b)$ , and a fraction $\lambda$ between 0 and 1,

Every point equation (12) produces is again a unit vector, so no step of the path leaves the sphere.

The trajectory is built in three parts (src/dsp/MatrixMorph.hpp). A base vector $r$ comes from the Diffusion control. $K$ anchor directions $u_k$ come from a fixed formula. There are six of them for the late field and four for the early stage. The mean of each is written $\overline{u_k}$ and is subtracted, so an anchor displaces $r$ without pulling it towards any one line.

Here $m$ is the mode number and $s$ is the loop width, which is 0.35 for Hall, 0.50 for Chamber, 0.25 for Cathedral and 0.80 for Drift. The $K$ corners $L_k$ define a closed loop. A number between 0 and 1 picks a point on that loop, by equation (12) between the two nearest corners. Morph Depth applies equation (12) a second time, between $r$ and that point. Depth 0 therefore returns $r$ exactly, and the matrix is static to the last bit.

One corner of every loop is the base vector itself. Take the last anchor, $k = K - 1$ , with $K$ equal to $N$ . The fraction inside the sine of equation (13) is then a whole number of turns for every entry $i$ . The same holds for $k = 2$ when $N = 3$ and $K = 4$ . The sine is the same for every entry, and the subtraction of the mean leaves nothing. So $L_5$ is $r$ in the late field and $L_2$ is $r$ in the early stage. Once per turn, five sixths of the way round, the late vector passes through its base position exactly. At full depth the loop reaches 31 degrees from the base in Hall, 50 in Chamber, 29 in Cathedral and 54 in Drift.

Moving the feedback matrix is not a new idea. Schlecht and Habets published it in 2015, modulating an FDN’s feedback matrix nearly continuously. A listening test in that paper judged the tail better for it.⁶ Their stability argument is the one equation (8) makes here. A unitary matrix conserves energy at every orientation, so modulating it cannot destabilise the loop. Their paper names earlier time-varying allpass networks as its own prior art.

What is a choice rather than a discovery is the parameterisation. SD-Reverb carries the reflection vector and moves that. Every point on the path is then length-preserving by construction, and none has to be corrected back. The design considered a rotation driven by a matrix exponential and set it aside for this.

A picture of the sphere is not available, and a picture of the path costs something measurable. Every point on the loop is an interpolation between two corners, so the whole path lies in the span of the six. The five non-zero anchors are independent, so the loop uses every one of the five directions available at $r$ .

How much a flat drawing would lose follows from the spread across those five. In Hall two directions hold 68 per cent of the movement and three hold 87. The four modes read 68 and 87 for Hall, 82 and 92 for Chamber, 72 and 85 for Cathedral, and 79 and 91 for Drift. A flat drawing would therefore fold between a fifth and a third of the movement onto itself. A closed loop drawn that way reads as an arc doubling back. The figure above plots the components and the angle instead, and both of those are exact.

Morph Rate spans 0.002 to 0.5 Hz on a logarithmic curve. The bottom two per cent of the control is reserved for exactly zero, so a frozen matrix is reachable. The default is 0.02 Hz, one turn of the loop in 50 seconds. The early stage runs the same construction at three entries and four corners, at one third of the rate and 0.7 of the depth. The vector is recomputed once per 64-sample chunk and handed to the engine that is sounding, and during a crossfade to both.

The matrix preserves length at every point on the path. Decay time therefore does not move, because equation (8) holds whatever the orientation. Nothing shifts in pitch, because no delay length changes and the movement is a redistribution of energy across lines. What does change is which resonances an incoming sound excites, and how strongly. The colour of the tail drifts while its length stays put.

Diffusion as vector balance

The base vector $r$ is a single control on how evenly the network mixes. For a balance $\beta$ between 0 and 1, the construction sets the first entry to $1 + \beta(N^{-1/2} - 1)$ and the others to $\beta N^{-1/2}$ , then scales the result to length 1.

At $\beta = 1$ every entry equals $N^{-1/2}$ and the matrix is the even Householder, the one the Airwindows fields use. At $\beta = 0$ the vector is $(1,0,0,0,0,0)$ and equation (9) gives a matrix that is zero everywhere except its diagonal, which reads $(-1,1,1,1,1,1)$ . That matrix inverts line 0 and leaves the other five untouched, so nothing mixes and the network is six separate comb filters. The Diffusion control moves between those two, and low settings are grainy for exactly that reason.

Long decays and low diffusion together produce audible ringing, so a floor is applied. Below a 4 second decay the floor is zero. Above it the floor rises as $(T_{60} - 4)/8$ , capped at 0.5. The same floored value is what blends the allpass chain into the late field’s input.

Frequency-dependent decay

Air absorbs high frequencies, soft surfaces absorb them further, and a large room applies more of both because sound travels further between reflections. A single lowpass filter in the feedback path reproduces the direction of that effect and little else, since it also flattens everything above its corner.

A shelf is the alternative. It changes the level of everything on one side of a corner frequency by a fixed amount and leaves the other side alone. Two shelves and the untouched middle then give three bands with three separate decay times. SD-Reverb gives each delay line a low shelf and a high shelf, with corners at 250 Hz and 4000 Hz (src/dsp/MultibandDamping.cpp). The mid band takes $T_{60}$ from the Decay control, and the low band takes 1.3 times it at Tone Tilt 0. The high band takes a multiplier that falls from 0.9 at Size 0 to 0.55 at Size 1. Equation (3) then turns each target into a gain.

At Size 0 the resulting decay is 3.11 s at 20 Hz, 2.43 s at 1 kHz and 2.23 s at 20 kHz. At Size 1 the top end falls to 1.51 s. Neither curve reaches its band target exactly, because a shelf built from one feedback term approaches its limit gradually. Tone Tilt moves the two multipliers in opposite directions. At a tilt $\tau$ the low one is scaled by $1 - 0.25\tau$ and the high one by $1 + 0.6\tau$ . The curve therefore leans dark at negative values and bright at positive ones.

Room tone and saturation

Each of the six delay lines has its own noise generator, built from an xorshift source, a pink filter and a rolloff above 8 kHz (src/dsp/RoomTone.hpp). Xorshift is a cheap random number generator. Pink noise has equal energy in every octave, which is closer to air movement than the flat hiss it starts as. The technique comes from the Airwindows reverbs, where it is called alive silence.

The level at the default Air Level of 0.30 is $10^{-4}$ , around 80 dB below full scale, which is the largest value a sample can hold. One shared 0.15 Hz oscillator scales the amplitude of all six, which is what the Air Motion control sets. Air Level 0 injects exactly zero.

The injection point is step five of the loop, after damping and before the matrix. The noise therefore circulates and picks up the room’s resonances, which is what makes it read as air in a space.

The saturation at step four maps each sample through a curve with a threshold $T$ of 3.

Here $\operatorname{sign}(x)$ is 1 for positive $x$ and $-1$ for negative, so the second branch holds the size and keeps the direction. The slope of the first branch is $1 - (x/T)^2$ , which falls from 1 at the origin to 0 at the threshold and never goes below zero. The gain therefore never exceeds 1 anywhere, which is what stops the curve from raising loop gain and destabilising the network. Tail values sit far below 3, so decay time is untouched at normal levels.

Entry Drift and Morph Adapt

The input matrix $B$ is 6 by 2. Left enters lines 0, 2 and 4 and right enters lines 1, 3 and 5, with the middle line of each triple inverted, at a weight of 0.5. The output matrix $C$ taps lines 0, 1 and 2 to the left channel and 3, 4 and 5 to the right, at $1/\sqrt{3}$ . Two channels share one network and neither output can be predicted from the other, which is what true stereo means.

Entry Drift moves $B$ on its own trajectory. The input matrix affects neither decay nor density. What moves is which lines an arriving sound excites, so two identical hits seconds apart ring different subsets of the same structure. A general 6 by 6 transformation would mix the left-fed and right-fed lines and drift the stereo image. What is applied is therefore a 3 by 3 Householder reflection, the same one to each triple. Each channel’s injection keeps its length exactly, because equation (10) holds on three entries as well as six.

A reflection is never the identity, so switching the feature on does not start from the fixed pattern and nudge it. At the base of its loop the reflection vector is the even three-entry vector. It turns the left triple’s weights from $(\tfrac12, -\tfrac12, \tfrac12)$ into $(\tfrac16, -\tfrac56, \tfrac16)$ , a different door of the same width. The loop then wanders around that pattern, at a depth scaled by 0.4. It always runs at the Hall width whatever the Mode. The change between the fixed and the reflected injection is crossfaded over 30 ms when the control leaves or returns to zero. The rate runs from 0.08 to 0.33 Hz and is set independently of Morph Rate. Whether the two loops ever return to a joint starting point depends on the ratio the user has set, and nothing in the plugin fixes it.

Morph Adapt gates the speed of the matrix trajectory from the input level. The envelope follower rises towards a new level over 50 ms and falls back over 1.2 seconds. The factor multiplying the trajectory speed is $1 - \min(4E, 1)\,\alpha$ , for envelope $E$ and Morph Adapt amount $\alpha$ . Loud input slows the matrix and silence lets it drift. The slow fall is what keeps the movement from becoming rhythmic, since a fast one would resume the drift in every gap between notes. The input trajectory is never gated this way, because the two features want opposite relationships to input activity.

Delay sets from Airwindows

The delay lengths are derived from Chris Johnson’s tuned tables and not copied from them. His integers are in undersampled cascade ticks and describe six stages in series. One SD-Reverb delay line therefore stands in for a whole round trip through his cascade. Dropping his numbers in at 48 kHz would give a metallic box under 10 ms across.

What was taken is the relative pattern of each source plugin’s six first-stage input taps. That pattern is rescaled log-affinely into the target millisecond window for each mode and rounded to integers. Log-affine rescaling puts the logarithms of the source taps onto the logarithms of the window with a straight line. The ratios between taps survive the move. The six are then adjusted until they are coprime in every pair, with the shortest late line at 15 ms or more. The provenance is recorded per mode in src/dsp/tables/DelaySets.hpp.

The Airwindows licence is MIT. It asks that the copyright notice travel with the work, and places no constraint on the licence of anything built from it. The notice is retained in SD-Reverb’s LICENSE, in CREDITS.md, and in the header of every file carrying a derived constant.

Deviations from the Airwindows design

Two mechanisms therefore keep the tail moving at once, and neither began here. Room tone perturbs the input to a system that is otherwise fixed, which is Chris Johnson’s technique. Matrix morphing varies the system itself, which is Schlecht and Habets’. Running the two together, on delay sets derived from his, is what is left.

Credit

The architecture is Chris Johnson’s. The two-stage shape comes from the Airwindows reverbs, as does the Householder matrix that mixes the late field. So do the room tone in the regeneration path and the practice of searching delay sets. They are MIT licensed and freely available. His notice travels with every file here that carries a derived constant. The delay lengths are rescaled from his tables.⁵

Moving that matrix over time is Schlecht and Habets’. They published it in 2015 with the stability argument and a listening test. Earlier time-varying allpass networks are the prior art they name in turn.⁶

What both of those rest on is older again. Gerzon’s cross-coupled combs, Stautner and Puckette’s stability conditions and Jot and Chaigne’s damping proportional to delay are the network itself. Schroeder’s combs in front of allpasses, Moorer’s lowpass inside the comb and Dattorro’s allpasses inside a single loop are what came before it. The figures quoted for those designs are from Julius O. Smith’s account, as is the history above.²

The rest is assembly, described in the sections above. The reflection vector is what moves, the damping scales with size, and three controls govern the motion.

Figure sources

Every image is drawn by a script kept beside this page.

figures.py writes the five computed plots. The comb response comes from the transfer function of equation (2). The damping curve evaluates the shelf pair of MultibandDamping.cpp at the constants FDN.cpp passes it, and converts the resulting per-round-trip gain back into a decay time. The density plot runs three fields. Those are the eight allpass sections of AllpassDiffuser.hpp, the six-line late field at a fixed gain of 0.985, and a six-stage cascade on kRockstar’s transcribed delays. The late field and the cascade both use the even Householder. Each is the structure alone, driven from one channel, and it then counts the arrivals in what each produces. The allpass plot feeds one click through one section and through eight. The impulse response is the one drawing with invented numbers.

diagrams.py writes the three block diagrams, the mixing grids and the trajectory plot. The trajectory is computed, running the anchor construction of equation (13) and both interpolations at six entries, mode 0 and six corners. The published version of the script used four corners (trajectory(..., k=4)) and must be run with k=6 to draw the loop the plugin runs.

figtheme.py holds the palette both of them draw with, and writes every file. It replaces each colour with a custom property and puts the light and dark values inside the file. A figure is therefore dark ink on a light page and light ink on a dark one.

All three are published with this page, so any of them can be read or run without asking for anything: figures.py, diagrams.py and figtheme.py.

uv run --with numpy --with matplotlib python3 figures.py
uv run --with numpy --with matplotlib python3 diagrams.py

References

• Audio
• Dsp
• C++