# Bigverb

This is a `sndkit`

algorithm. A more up-to-date version
can be found here.

### Introduction

Bigverb is an implementation of a classic waveguide reverberator. It is composed up of 8 feedback delay lines running in parallel, in what is known as a Feedback Delay Network. It is largely inspired by the works Julius O. Smith, and the work he did in using waveguide networks. In particular, this algorithm is based on ideas in a 1985 paper called "A New Approach to Digital Reverberation using Waveguide Networks".

Bigverb, as the name implies, has a very nice warm sound when set to be a very large space, and is ideal for use in ambient style music. Bigverb is not designed to sound realistic, and is not an ideal choice for anything that needs room modelling.

### Distinct Characterstics

The "sound" of this algorithm comes from a few major things.

#### Use of Jitter Modulation

Every delay line has a bit of jitter applied to the delay time, which causes a little of pitch modulation. Several schroeder and FDN-style reverb designs use this kind of modulation, such as the Mutable instruments clouds reverb, and the tank reverb described in the Dattorro paper "Effect design, part 1: reverberator and other filters". Such modulation creates the illusion of the reverb being a higher-order system than it actual (aka: a more complex reverb algorithm with more delay lines). While these designs tend to use a sinusoidal LFO on only a few delay lines, the Costello topology applies non-correlated modulation to every single delay line. As a result, you get a very dense-sounding reverb output.

#### Parameter Tunings

The parameter tunings of this reverb add to the particular character of this sound. All units are set in such a way so that they can be stored as integer values.

The tuning consists of 8 parameter sets, corresponding to
each delay line "unit". A set has the following parameters: delay time
(`delay`

) (in units of samples with a samplerate of
44.1kHz), the amount of variation in delay time to apply
(`drift`

) (in units of tenth milliseconds, in order to use
integer values), the random variation frequency (`randfreq`

)
(in Hz multiplied by 1000 to avoid decimals),
as well as a starting seed for the internal RNG.

The credit for these parameter tunings (and this topology) go to Sean Costello, who also originally designed this reverb algorithm in 1999. This work would later pave the way for his famous commercial reverb plugins and algorithms published by his company ValhallaDSP.

*<<parameter_set>>=*```
struct bigverb_paramset {
int delay; /* in samples, 44.1 kHz */
int drift; /* 1/10 milliseconds */
int randfreq; /* Hertz * 1000 */
int seed;
};
static const struct bigverb_paramset params[8] = {
{0x09a9, 0x0a, 0xc1c, 0x07ae},
{0x0acf, 0x0b, 0xdac, 0x7333},
{0x0c91, 0x11, 0x456, 0x5999},
{0x0de5, 0x06, 0xf85, 0x2666},
{0x0f43, 0x0a, 0x925, 0x50a3},
{0x101f, 0x0b, 0x769, 0x5999},
{0x085f, 0x11, 0x37b, 0x7333},
{0x078d, 0x06, 0xc95, 0x3851}
};
```

#### Delay Line Configuration

Normally, this family of reverbs is a combination of allpass
filters and comb filters in various series/parallel
combinations. This implementation has a far simpler
topology: it is essentially (but not **quite**) 8 feedback delay
lines in parallel with a 1-pole lowpass IIR filter for tone
control. If it weren't for the massive amount delay
modulation, this would be a pretty metallic sounding reverb!

### Top Level Files and Structs

#### Files to Tangle To

This Bigverb document tanlges to two files:
`bigverb.c`

and `bigverb.h`

.

`bigverb.h`

contains typedefs, structs, and forward
function declarations. To enable structs, define the
macro `SK_BIGVERB_PRIV`

.

*<<bigverb.c>>=*```
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define SK_BIGVERB_PRIV
#include "bigverb.h"
<<parameter_set>>
<<fdelay_constants>>
<<static_funcdefs>>
<<funcs>>
```

*<<bigverb.h>>=*```
#ifndef SK_BIGVERB_H
#define SK_BIGVERB_H
#ifndef SKFLT
#define SKFLT float
#endif
<<typedefs>>
<<funcdefs>>
#ifdef SK_BIGVERB_PRIV
<<main_struct>>
#endif
#endif
```

#### The Bigverb Struct

An instance of Bigverb is contained inside of a struct
called `sk_bigverb`

.

*<<typedefs>>=*`typedef struct sk_bigverb sk_bigverb;`

*<<main_struct>>=*```
<<delay_struct>>
struct sk_bigverb {
int sr;
<<sk_bigverb>>
};
```

### Setup and Cleanup

A new instance of bigverb is created with
`sk_bigverb_new`

. The only argument required is the sampling
rate. If something goes wrong, this will return `NULL`

.

*<<funcdefs>>=*`sk_bigverb * sk_bigverb_new(int sr);`

*<<funcs>>=*```
sk_bigverb * sk_bigverb_new(int sr)
{
sk_bigverb *bv;
bv = calloc(1, sizeof(sk_bigverb));
bv->sr = sr;
<<init_variables>>
<<setup_delay_lines>>
return bv;
}
```

When it is done being used, bigverb must be cleanly freed
with `sk_bigverb_del`

.

*<<funcdefs>>=*`void sk_bigverb_del(sk_bigverb *bv);`

*<<funcs>>=*```
void sk_bigverb_del(sk_bigverb *bv)
{
<<cleanup>>
free(bv);
bv = NULL;
}
```

### High level parameters

High level parametric control of bigverb includes "size" and "cutoff". Set parameters before computing audio. These are are just floating point values contained in the struct that can be indirectly set with setters in situations where the struct is opaque.

#### Size

Set the reverb size with `sk_bigverb_size`

*<<funcdefs>>=*`void sk_bigverb_size(sk_bigverb *bv, SKFLT size);`

Size is a variable between 0-1, which controls the feedback level for the delay line.

The `size`

parameter is stored as a variable called
`size`

, and is set to be a pretty sounding value of
`0.93`

.

*<<sk_bigverb>>=*`SKFLT size;`

*<<init_variables>>=*`sk_bigverb_size(bv, 0.93);`

*<<funcs>>=*```
void sk_bigverb_size(sk_bigverb *bv, SKFLT size)
{
bv->size = size;
}
```

#### Cutoff

The tone of bigverb can be set with `sk_bigverb_cutoff`

.

*<<funcdefs>>=*`void sk_bigverb_cutoff(sk_bigverb *bv, SKFLT cutoff);`

`cutoff`

is a parameter in Hz that determines the overall
timbre of the reverb. This controls the cutoff frequency of
the one pole lowpass filter applied to the reverb.

It is set to be a default value of 10kHz, or 10,000 hz.

*<<init_variables>>=*`sk_bigverb_cutoff(bv, 10000.0);`

Cutoff uses caching in order to monitor if the parameter
has changed. It does this in order to prevent needing to
compute filter coefficients every sample. The main variable
to be set is `cutoff`

, and the cached variable is `pcutoff`

.
At the beginning, `pcutoff`

is set to be a negative value,
which will cause bigverb to calculate coefficients in the
first call to the tick function after initialization.

*<<sk_bigverb>>=*```
SKFLT cutoff;
SKFLT pcutoff;
```

*<<init_variables>>=*`bv->pcutoff = -1;`

*<<funcs>>=*```
void sk_bigverb_cutoff(sk_bigverb *bv, SKFLT cutoff)
{
bv->cutoff = cutoff;
}
```

### Filter

State in a constant called `filt`

.

*<<sk_bigverb>>=*`SKFLT filt;`

*<<init_variables>>=*`bv->filt = 1.0;`

### Computing Audio

After bigverb has been initialized, it is ready to
process audio. This implementation uses what is known
as a `tick`

function, or a function that computes audio
one sample at a time instead of one block at a time. This
simplifies the implementation at the cost of a little bit
of performance overhead, depending on the compiler and
optimization settings.

#### Top-Level Tick Function

The function to tick one sample unit of audio is done with
`sk_bigverb_tick`

. It takes in two stereo input
values, and returns two stereo output values.

*<<funcdefs>>=*```
void sk_bigverb_tick(sk_bigverb *bv,
SKFLT inL, SKFLT inR,
SKFLT *outL, SKFLT *outR);
```

*<<funcs>>=*```
void sk_bigverb_tick(sk_bigverb *bv,
SKFLT inL, SKFLT inR,
SKFLT *outL, SKFLT *outR)
{
/* TODO: implement */
SKFLT lsum, rsum;
lsum = 0;
rsum = 0;
<<update_filter_coefficients>>
<<calculate_junction_pressure>>
<<compute_delay_bank>>
*outL = lsum;
*outR = rsum;
}
```

#### Updating filter coefficients

Bigverb uses parameter caching for the `cutoff`

parameter in
order to save on computation time.

Any time `cutoff`

changes, the filter coefficients must be
updated. This happens in the tick function, before any
computation happens.

The filter is a simple 1-pole IIR lowpass filter whose difference equation been reduced to only require a single parameter. This in turn then gets used in each filter delay line.

*<<update_filter_coefficients>>=*```
if (bv->pcutoff != bv->cutoff) {
bv->pcutoff = bv->cutoff;
bv->filt = 2.0 - cos(bv->pcutoff * 2 * M_PI / bv->sr);
bv->filt = bv->filt - sqrt(bv->filt * bv->filt - 1.0);
}
```

#### Calculating Resultant Junction Pressure Amount

The resultant junction pressure amount is calculated from the delay bank, and then factored into the input signals.

Sum of all the delay line signals, and scaled by 0.25, or 2/N, where N is the number of delay lines (8).

*<<calculate_junction_pressure>>=*```
{
int i;
SKFLT jp;
jp = 0;
for (i = 0; i < 8; i++) {
jp += bv->delay[i].y;
}
jp *= 0.25;
inL = jp + inL;
inR = jp + inR;
}
```

#### Computing the delay bank

The delay bank is then computed. Each delay line is computed and summed with either the left or right input signal, and then sent to a corresponding left or right channel.

At the end, a final scaling out of the output happens. This is hard coded to be 35 percent.

*<<compute_delay_bank>>=*```
{
int i;
for (i = 0; i < 8; i++) {
if (i & 1) {
rsum += delay_compute(&bv->delay[i],
inR,
bv->size,
bv->filt,
bv->sr);
} else {
lsum += delay_compute(&bv->delay[i],
inL,
bv->size,
bv->filt,
bv->sr);
}
}
}
rsum *= 0.35f;
lsum *= 0.35f;
```

### The Feedback Delay Line Bank

8 delay units come together to make the delay line bank. Each is initialized using one of the parameter sets.

#### Memory Allocation + Setup

*<<sk_bigverb>>=*`SKFLT *buf;`

*<<init_variables>>=*`bv->buf = NULL;`

*<<setup_delay_lines>>=*```
{
unsigned long total_size;
int i;
SKFLT *buf;
total_size = 0;
buf = NULL;
<<calculate_pool_size>>
<<allocate_memory>>
<<initialize_delay_banks>>
}
```

The delay bank is the abstraction in charge of properly allocating all the memory needed for the buffers.

Memory is allocated in one giant chunk, and then divied up to each delay line.

The total memory size is obtained by summing all the delay times. These times are stored as fixed delay times in units of samples. These parameters assume a sampling rate of 44.1kHz. If this is not the case, this value must be scaled accordingly, and then truncated to be an integer. This value is used again to properly slice up the big memory chunk.

*<<static_funcdefs>>=*`static int get_delay_size(const struct bigverb_paramset *p, int sr);`

*<<funcs>>=*```
static int get_delay_size(const struct bigverb_paramset *p, int sr)
{
SKFLT sz;
sz = (SKFLT)p->delay/44100 + (p->drift * 0.0001) * 1.125;
return floor(16 + sz*sr);
}
```

*<<calculate_pool_size>>=*```
for (i = 0; i < 8; i++) {
total_size += get_delay_size(¶ms[i], sr);
}
```

Allocation is done with `calloc`

, which zeros out the memory
as well. This memory will eventually be freed in
`sk_bigverb_del`

.

*<<allocate_memory>>=*```
buf = calloc(1, sizeof(SKFLT) * total_size);
bv->buf = buf;
```

*<<cleanup>>=*`free(bv->buf);`

*<<sk_bigverb>>=*`sk_bigverb_delay delay[8];`

*<<initialize_delay_banks>>=*```
{
unsigned long bufpos;
bufpos = 0;
for (i = 0; i < 8; i++) {
unsigned int sz;
sz = get_delay_size(¶ms[i], sr);
delay_init(&bv->delay[i], ¶ms[i],
&buf[bufpos], sz, sr);
bufpos += sz;
}
}
```

#### A Single Delay Line Unit

A delay unit in a bank consists of variable delay line with cubic interpolation with a 1 pole low-pass filter for tone control, whose frequency is determined using a master parameter, as well as a jitter generator. Feedback as well.

##### Struct Declaration

A delay unit is known as a struct called `sk_bigverb_delay`

.

*<<typedefs>>=*`typedef struct sk_bigverb_delay sk_bigverb_delay;`

*<<delay_struct>>=*```
struct sk_bigverb_delay {
<<bigverb_delay>>
};
```

##### Initialization

*<<static_funcdefs>>=*```
static void delay_init(sk_bigverb_delay *d,
const struct bigverb_paramset *p,
SKFLT *buf,
size_t sz,
int sr);
```

*<<funcs>>=*```
static void delay_init(sk_bigverb_delay *d,
const struct bigverb_paramset *p,
SKFLT *buf,
size_t sz,
int sr)
{
SKFLT readpos;
<<delay_init>>
}
```

Set up buffer + sz

*<<bigverb_delay>>=*```
SKFLT *buf;
size_t sz;
```

*<<delay_init>>=*```
d->buf = buf;
d->sz = sz;
```

Initialize write position (0), abbreviated as `wpos`

.

*<<bigverb_delay>>=*`int wpos;`

*<<delay_init>>=*`d->wpos = 0;`

Initialize read position. Based on delay time, `drift`

and
initial seed. Read position has to components, an integer
read position, and a floating point read position. These
will be abbreviated `irpos`

and `frpos`

.

*<<bigverb_delay>>=*```
int irpos;
int frpos;
```

Seed value is multiplied by the initial drift value, and then divided by 32767.

*<<bigverb_delay>>=*`int rng;`

*<<delay_init>>=*```
d->rng = p->seed;
<<setup_readpos>>
```

The initial time is added to this.

`bufsize - (readpos * sr)`

<-- this puts the read position at
the end of the buffer.

Truncate (using integer cast).'

*<<setup_readpos>>=*```
readpos = ((SKFLT)p->delay / 44100);
readpos += d->rng * (p->drift * 0.0001) / 32768.0;
readpos = sz - (readpos * sr);
d->irpos = floor(readpos);
d->frpos = floor((readpos - d->irpos) * FRACSCALE);
```

Create first random segments.

*<<delay_init>>=*```
<<init_jitter>>
generate_next_line(d, sr);
```

##### Top-Level Compute

The delay line computation is done in a tick function. It takes in an input sample, returns an output sample. In addition to delay, filtering, feedback, and jittering happens as well.

Because feedback + filtering are global options, these are passed in as parameters on the stack. What is required is the feedback amount, and the calculated filter coeffecient used in the filter.

*<<static_funcdefs>>=*```
static SKFLT delay_compute(sk_bigverb_delay *d,
SKFLT in,
SKFLT fdbk,
SKFLT filt,
int sr);
```

*<<funcs>>=*```
static SKFLT delay_compute(sk_bigverb_delay *del,
SKFLT in,
SKFLT fdbk,
SKFLT filt,
int sr)
{
SKFLT out;
SKFLT frac_norm;
SKFLT a, b, c, d;
SKFLT s[4];
out = 0;
<<write_to_buffer>>
<<increment_write_position>>
<<update_fractional_read_position>>
<<update_integer_read_position>>
<<normalize_fractional_component>>
<<calculate_interpolation_coefficients>>
<<read_from_buffer>>
<<compute_interpolation>>
<<increment_fractional_read_position>>
<<apply_feedback_and_filter>>
<<update_jitter>>
return out;
}
```

The following things happen:

Write the to delay buffer and pre-filter the input by
subtracting the filter state `y`

.

*<<write_to_buffer>>=*`del->buf[del->wpos] = in - del->y;`

Increment the write position. If this is greater than the buffer size, wrap around.

*<<increment_write_position>>=*```
del->wpos++;
if (del->wpos >= del->sz) del->wpos -= del->sz;
```

Update the fractional read position. If the read position exceeds the maximum fractional scale amount, it means it has bits that must carry over to the integer read position. After these bits have been carried over, mask out the upper bits to keep the range in bounds.

*<<update_fractional_read_position>>=*```
if (del->frpos >= FRACSCALE) {
del->irpos += del->frpos >> FRACNBITS;
del->frpos &= FRACMASK;
}
```

If needed, update the read position with wrap-around.

*<<update_integer_read_position>>=*`if (del->irpos >= del->sz) del->irpos -= del->sz;`

Normalize the fractional component so that it is in range
0 and 1. This is done by dividing the amount by fractional
scaling factor `FRACSCALE`

.

*<<normalize_fractional_component>>=*`frac_norm = del->frpos / (SKFLT)FRACSCALE;`

Calculate interpolation coefficients. These are 4
pre-derived coefficents used to compute third-order
lagrangian interpolation.
Derivation of these is currently beyond the scope of this
document. These will be called `a`

, `b`

, `c`

, and `d`

,
respectively, and will correspond to `x(n - 1)`

, `x(n)`

,
`x(n + 1)`

, and `x(n + 2)`

, respectively.

*<<calculate_interpolation_coefficients>>=*```
{
SKFLT tmp[2];
d = ((frac_norm * frac_norm) - 1) / 6.0;
tmp[0] = ((frac_norm + 1.0) * 0.5);
tmp[1] = 3.0 * d;
a = tmp[0] - 1.0 - d;
c = tmp[0] - tmp[1];
b = tmp[1] - frac_norm;
}
```

Read the samples needed, based on the current playhead position. When the read position is in regular bounds, this means reading the previous, current, two next samples. Otherwise, this means the same thing, but with wrap-around and bounds checks.

*<<read_from_buffer>>=*```
{
int n;
SKFLT *x;
n = del->irpos;
x = del->buf;
if (n > 0 && n < (del->sz - 2)) {
s[0] = x[n - 1];
s[1] = x[n];
s[2] = x[n + 1];
s[3] = x[n + 2];
} else {
int k;
n--;
if (n < 0) n += del->sz;
s[0] = x[n];
for (k = 0; k < 3; k++) {
n++;
if (n >= del->sz) n -= del->sz;
s[k + 1] = x[n];
}
}
}
```

Calculate interpolation. Using the coefficents
described above and the fractional component `f`

, one can compute
cubic interpolation with the following expression:

$$ y(n) = (a x(n - 1) + b x(n) + c x(n + 1) + d x(n + 2)) \cdot f + x(n) $$

*<<compute_interpolation>>=*`out = (a*s[0] + b*s[1] + c*s[2] + d*s[3]) * frac_norm + s[1];`

Increment fractional read position, as determined by the jitter.

*<<increment_fractional_read_position>>=*`del->frpos += del->inc;`

Apply feedback and filter. The feedback will scale the delay output. The filtering is a difference equation, optimized and factored to use a minimum number of arithmetic operations.

*<<apply_feedback_and_filter>>=*```
out *= fdbk;
out += (del->y - out) * filt;
del->y = out;
```

# # #

Update jitter, if needed. When the counter zeros out (or worse), it is time to find a new random target to lerp to.

*<<update_jitter>>=*```
del->counter--;
if (del->counter <= 0) {
generate_next_line(del, sr);
}
```

##### Feedback Fractional Delay Line

A delay line is initialized with a pre-allocated zeroed buffer its size. Memory will be managed outside of this abstraction.

Being a fractional delay line means the read position has two components: an integer component and a fractional component. The integer component is the current position in the delay buffer. The fractional component tells how much it goes over into the next discrete sample position. In a way, interpolation can be thought of as the process of using these two values to make a really good guess of what lies in between the samples.

The fractional delay component is maintained as a 28 bit integer. This is done to avoid some of the weirdness found in floating point operations. The remaining upper bits are "carry-over" samples, that get accumulated in integer component of the read position.

A few constants are used to conveniently work with this fractional delay component.

`FRACSCALE`

is the fractional scaling amount,
which is $1^28$, or `0x10000000`

. Multiplied with a uniform
scalar, this is used to calculate the increment.

*<<fdelay_constants>>=*`#define FRACSCALE 0x10000000`

`FRACMASK`

is the bitmask used to keep the fractional
position in 28-bit range. It is $1^28 - 1$, or `0xFFFFFFF`

.
This is in particular is used to filter out upper bits that
get carried over to the integer read position.

*<<fdelay_constants>>=*`#define FRACMASK 0xFFFFFFF`

`FRACNBITS`

is the number of bits in the number. set to be
28.

*<<fdelay_constants>>=*`#define FRACNBITS 28`

##### Jitter

Jitter in this context, is a random line segment generator. It linearly interpolates between random values in a given range, using random durations in a given range.

A line generator stores a counter and increment amount.

*<<bigverb_delay>>=*```
int inc;
int counter;
```

*<<init_jitter>>=*```
d->inc = 0;
d->counter = 0;
```

The most significant thing to happen in the jitter signal is
calculating the next random segment. This is done in a
static function called `generate_next_line`

.

*<<static_funcdefs>>=*`static void generate_next_line(sk_bigverb_delay *d, int sr);`

To begin, another random value is created based on the previous random value.

The RNG algorithm used is quite simple, and is used to produce a 16-bit value.

$$x(n) = (5^{6}x(n - 1) + 1) \& 0xFFFF$$

Before and after this equation, the value is balanced so that is a 16-bit bipolar signal. Before, it adds 0x10000 if the value is less than 0. After, it substracts 0x10000 if the seed value is greater than 0x8000.

*<<generate_random_number>>=*```
if (d->rng < 0) d->rng += 0x10000;
/* 5^6 = 15625 */
d->rng = (1 + d->rng * 0x3d09);
d->rng &= 0xFFFF;
if (d->rng >= 0x8000) d->rng -= 0x10000;
```

This new random value is used to produce the next random value in seconds.

The line counter is reset. This value comes from the high-level parameter.

*<<bigverb_delay>>=*`int maxcount;`

*<<init_jitter>>=*`d->maxcount = round((sr / ((SKFLT)p->randfreq * 0.001)));`

*<<reset_counter>>=*`d->counter = d->maxcount;`

Compute delay time values. The current delay time,
`curdel`

, is obtained by subtracting the write + integer
read positions, then adding in the fractional component.
Wraparound is applied.

*<<compute_delay_values>>=*```
curdel = d->wpos -
(d->irpos + (d->frpos/(SKFLT)FRACSCALE));
while (curdel < 0) curdel += d->sz;
curdel /= sr;
```

The next delay time to lerp to is derived from the RNG and drift amount.

*<<compute_delay_values>>=*`nxtdel = (d->rng * (d->drift * 0.0001) / 32768.0) + d->dels;`

The delay time, in seconds (`dels`

).

*<<bigverb_delay>>=*`SKFLT dels;`

*<<init_jitter>>=*`d->dels = p->delay / 44100.0;`

*<<bigverb_delay>>=*`SKFLT drift;`

*<<init_jitter>>=*`d->drift = p->drift;`

The linear increment value is the difference between the current
and next delay times, divided by the number of steps needed
to draw a line between them (`counter`

). This value is then
converted into samples. An extra sample is tacked on to
prevent nil values.

*<<compute_increment>>=*```
inc = ((curdel - nxtdel) / (SKFLT)d->counter)*sr;
inc += 1;
```

This increment value is truncated and converted to the fractional read.

*<<set_fractional_read>>=*`d->inc = floor(inc * FRACSCALE);`

*<<funcs>>=*```
static void generate_next_line(sk_bigverb_delay *d, int sr)
{
SKFLT curdel;
SKFLT nxtdel;
SKFLT inc;
<<generate_random_number>>
<<reset_counter>>
<<compute_delay_values>>
<<compute_increment>>
<<set_fractional_read>>
}
```

##### Filter Memory

A one-pole lowpass filter such as the one used in the delay
line requires one sample of memory, which stores the output
of the previous filter. In a difference equation, this would
be known as $y(n - 1)$. In C code, we abbreviate this as
`y`

.

*<<bigverb_delay>>=*`SKFLT y;`

*<<delay_init>>=*`d->y = 0.0;`