Audio DSP

Converts between the decibel values used by meters, compressors and gain controls and the linear values the audio maths works in.

Use 20 for amplitude. The equivalent for power values uses 10.

The lower limit prevents log10(0). At 1e-12 the floor sits at -240 dB.

gain = pow(10.0, gainDb / 20.0);

levelDb = 20.0 * log10(max(level, 1e-12));

Saturation with level compensation, so raising the amount does not just make it louder.

x = input sample, p = amount (0 to 1).

The x * x term leaves a DC offset behind, so follow it with a DC blocker.

drive  = 1.0 + p * 3.0;
warmth = p * 0.25;

sat  = tanh(x * drive * 0.9);   // symmetric, odd harmonics
harm = x * x * warmth;          // asymmetric, even harmonics

output = (sat + harm) * (1.0 / (1.0 + p * 0.8));

Nonlinear processing creates harmonics above the input frequency. When those cross Nyquist they fold back as aliasing, usually heard as a brittle high end. Oversample around the shaper when strong drive or bright material makes it audible.

More on aliasing

Removes the offset any asymmetric shaper leaves behind.

cornerHz around 5 is a safe choice. Work R out when the sample rate changes, not per sample.

x1 and y1 are the previous input and output. One set per channel.

R = exp(-2.0 * pi * cornerHz / sampleRate);

y  = x - x1 + R * y1;
x1 = x;
y1 = y;

output = y;

One set of formulas covering every band in an EQ. The standard reference, from Robert Bristow-Johnson's cookbook.

f0 = frequency, Q = resonance or bandwidth, dBgain for peaking and shelving only.

Work the coefficients out when a parameter changes, then run the last two lines per sample.

w0    = 2.0 * pi * f0 / sampleRate;
cosw  = cos(w0);
alpha = sin(w0) / (2.0 * Q);
A     = pow(10.0, dBgain / 40.0);   // peaking and shelving only
t     = 2.0 * sqrt(A) * alpha;      // shelving only

// low-pass
b0 = (1.0 - cosw) * 0.5;  b1 = 1.0 - cosw;     b2 = b0;
a0 = 1.0 + alpha;         a1 = -2.0 * cosw;    a2 = 1.0 - alpha;

// high-pass
b0 = (1.0 + cosw) * 0.5;  b1 = -(1.0 + cosw);  b2 = b0;
a0 = 1.0 + alpha;         a1 = -2.0 * cosw;    a2 = 1.0 - alpha;

// band-pass, 0 dB peak
b0 = alpha;               b1 = 0.0;            b2 = -alpha;
a0 = 1.0 + alpha;         a1 = -2.0 * cosw;    a2 = 1.0 - alpha;

// notch
b0 = 1.0;                 b1 = -2.0 * cosw;    b2 = 1.0;
a0 = 1.0 + alpha;         a1 = -2.0 * cosw;    a2 = 1.0 - alpha;

// peaking
b0 = 1.0 + alpha * A;     b1 = -2.0 * cosw;    b2 = 1.0 - alpha * A;
a0 = 1.0 + alpha / A;     a1 = -2.0 * cosw;    a2 = 1.0 - alpha / A;

// low shelf
b0 =      A * ((A+1) - (A-1)*cosw + t);
b1 =  2.0*A * ((A-1) - (A+1)*cosw);
b2 =      A * ((A+1) - (A-1)*cosw - t);
a0 =           (A+1) + (A-1)*cosw + t;
a1 = -2.0   * ((A-1) + (A+1)*cosw);
a2 =           (A+1) + (A-1)*cosw - t;

// high shelf
b0 =      A * ((A+1) + (A-1)*cosw + t);
b1 = -2.0*A * ((A-1) + (A+1)*cosw);
b2 =      A * ((A+1) + (A-1)*cosw - t);
a0 =           (A+1) - (A-1)*cosw + t;
a1 =  2.0   * ((A-1) - (A+1)*cosw);
a2 =           (A+1) - (A-1)*cosw - t;

// per sample, direct form 1
y  = (b0/a0)*x + (b1/a0)*x1 + (b2/a0)*x2 - (a1/a0)*y1 - (a2/a0)*y2;
x2 = x1; x1 = x; y2 = y1; y1 = y;

Q above about 10 at a low f0 puts the poles very close to the unit circle, where the filter rings and becomes sensitive to rounding. Direct form 1 holds up better there than the transposed forms.

Audio EQ Cookbook

The cheapest useful filter: 6 dB per octave, no resonance, one sample of state. Good for tone controls, smoothing and taming a feedback path.

Both share the same coefficient, so a pair costs one exp. Work it out when the cutoff changes.

Measured -3 dB at the cutoff, as expected.

x = exp(-2.0 * pi * cutoffHz / sampleRate);

// low-pass
y1 = (1.0 - x) * in + x * y1;
lowOut = y1;

// high-pass, same x
b0 = 0.5 * (1.0 + x);
hpY = b0 * in - b0 * hpX1 + x * hpY;
hpX1 = in;
highOut = hpY;

Four waveforms plus sample and hold, all read from one phase running 0 to 1.

Advance with phase += freqHz / sampleRate, and wrap by subtracting 1.0 when it passes 1.

Trigger the sample and hold on that wrap. Use your own random source, not rand(), inside the audio callback.

sine     = sin(2.0 * pi * phase);
triangle = 2.0 * fabs(2.0 * (phase - floor(phase + 0.5))) - 1.0;
square   = (phase < 0.5) ? 1.0 : -1.0;
saw      = 2.0 * phase - 1.0;

// sample and hold: refresh only when the phase wraps
if (wrapped) shValue = nextRandom();   // -1 to +1
sampleHold = shValue;

Tracks the level of a signal. The core of every compressor, gate, ducker and envelope-driven filter.

Rises fast and falls slow by using a different coefficient in each direction.

coeff = 1 - exp(-1 / (timeMs * sampleRate / 1000)), worked out once per time value.

level = fabs(x);

if (level > env)
    env += attackCoeff * (level - env);
else
    env += releaseCoeff * (level - env);

output = env;

This defines the stated time as one time constant, the point where 63.2% of the movement is done. Some processors instead quote the time to settle within -60 dB, which is about seven time constants, so coefficients copied from elsewhere may look different without being wrong.

Stops a knob move or automation step from clicking. Apply to any parameter read per block.

coeff = exp(-1 / (timeMs * sampleRate / 1000)). Around 20 ms suits most controls.

Smooth the value, not the audio: run this on the parameter before it reaches the maths.

current = current * coeff + target * (1.0 - coeff);

output = current;

Same time-constant convention as the envelope follower: the stated time is 63.2% of the way there, not fully settled.

Maps an input level in dB to an output level in dB, bending gradually into the ratio instead of switching at the threshold.

x, threshold and knee in dB, ratio as N to 1. Feed it the envelope, not the sample.

Continuous at both knee edges, so there is no step where the regions meet.

if (x <= threshold - knee * 0.5)
    y = x;                                  // below the knee
else if (x >= threshold + knee * 0.5)
    y = threshold + (x - threshold) / ratio;  // fully compressed
else {
    d = x - threshold + knee * 0.5;
    t = d / knee;
    y = x + (1.0 / ratio - 1.0) * d * t * 0.5;
}

gainDb = y - x;   // apply this, not y

One stage of an ADSR. Curves towards the target instead of sliding linearly, which is what makes it sound natural.

N = length of the segment in samples. reach = how close to the target it gets by the end, e.g. 0.01.

mul is computed once when the time changes, not per sample.

mul = pow(reach, 1.0 / N);      // once, when the time changes

value = target + (value - target) * mul;   // per sample

Turns stereo into a centre channel and a difference channel. The basis of mid/side EQ, width controls and processing centred and wide material separately.

This halves while encoding and uses unity gain while decoding. Others use 1 / sqrt(2) in both directions; either works, but do not mix them.

Width above 1 widens the result and can raise the decoded peak, so watch the output meter.

mid  = (left + right) * 0.5;
side = (left - right) * 0.5;

side *= width;    // 0 = mono, 1 = unchanged

left  = mid + side;
right = mid - side;

What Mid and Side actually contain

Moves a mono source across the stereo field without it dipping in the middle.

pan runs -1 (hard left) to +1 (hard right). Centre gives 0.707 on both sides, not 0.5.

The two gains satisfy L squared plus R squared equal to 1 at every position, which is what keeps the perceived level steady.

theta = ((pan + 1.0) * 0.5) * (pi / 2.0);

leftGain  = cos(theta);
rightGain = sin(theta);

outL = in * leftGain;
outR = in * rightGain;

Noise falling at 3 dB per octave, which is closer to how the ear weights a spectrum than white noise.

Paul Kellet's economy method. b0 to b6 are persistent state, white is a fresh random value in -1 to +1.

The final multiply keeps it roughly inside -1 to +1.

b0 = 0.99886 * b0 + white * 0.0555179;
b1 = 0.99332 * b1 + white * 0.0750759;
b2 = 0.96900 * b2 + white * 0.1538520;
b3 = 0.86650 * b3 + white * 0.3104856;
b4 = 0.55000 * b4 + white * 0.5329522;
b5 = -0.7616 * b5 - white * 0.0168980;

pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362;
b6 = white * 0.115926;

output = pink * 0.11;

Delay and LFO times that land on the grid.

noteValue: 1.0 = quarter, 0.5 = eighth, 0.25 = sixteenth.

Multiply by 1.5 for dotted, by 2/3 for a triplet. At 120 BPM a dotted eighth is 375 ms.

ms = (60000.0 / bpm) * noteValue;

Note 69 is A4 at 440 Hz. Change the 440.0 to retune the whole instrument.

Pitch bend arrives as 14 bits, 0 to 16383, centred at 8192.

range is the bend width in semitones, usually 2 or 12.

freq = 440.0 * pow(2.0, (note - 69) / 12.0);

// with pitch bend
semitones = ((bend - 8192.0) / 8192.0) * range;
freq = freq * pow(2.0, semitones / 12.0);

Clears inaudibly small persistent state before it eventually reaches the subnormal range, where some processors fall into a slow path.

Run it on state that decays towards silence, such as filter memory or a feedback buffer, not on every sample of the audio.

1e-15 is roughly -300 dB, so this is a practical silence cutoff rather than the actual IEEE denormal boundary, which is far lower.

if (fabs(x) < 1e-15)
    x = 0.0;