🎯 Why naive discretization fails

In the prior article on synthesizer filters, we explored analog VCF topologies and mentioned that "direct discretization often fails." This article explains exactly why—and presents the solution that modern virtual analog synths use: zero-delay feedback (ZDF), also called the topology-preserving transform (TPT).

Consider a simple one-pole lowpass filter. In the analog domain, it's described by the differential equation:

where ω_c = 2πf_c is the angular cutoff frequency. The output y follows the input x with a time constant determined by ω_c. The naive approach is to discretize this using the bilinear transform or simple backward/forward Euler differences. This gives you working filter coefficients—but with problems that become severe in synthesizer applications:

These problems are manageable for EQ and mixing tasks where you set parameters once and leave them. But in a synthesizer, the filter is a dynamic instrument—cutoff sweeps across the full range, resonance goes to self-oscillation, and modulation sources run at audio rates. Standard IIR filters fall apart under these demands.

🔄 Zero-delay feedback: the core idea

Zero-delay feedback (ZDF) solves the delay problem by solving for the output algebraically within each sample, rather than relying on the previous sample's output. The key insight comes from Vadim Zavalishin's work on the topology-preserving transform: instead of discretizing the transfer function, we discretize the circuit topology directly, preserving the instantaneous feedback relationships.

The mathematical foundation uses the trapezoidal integrator. An analog integrator has the transfer function H(s) = 1/s. The trapezoidal rule discretizes this as:

This can be rewritten using a state variable that captures the "memory" of the integrator:

where g = tan(π·f_c/f_s) is the prewarped coefficient and s is the integrator state. The tan() function provides automatic frequency prewarping—it compensates for the frequency compression that the bilinear transform introduces, so the filter tracks the specified cutoff accurately up to Nyquist.

📐 TPT state variable filter: step by step

Let's derive the topology-preserving state variable filter from scratch. The analog SVF consists of two integrators in series with feedback. In the s-domain, it has lowpass, bandpass, and highpass outputs simultaneously. We'll build the digital version that preserves this topology.

Step 1: The analog block diagram

An analog SVF has this signal flow:

• Input x feeds into a summer
• The summer also receives the negated lowpass output (scaled by 2R for resonance) and the negated bandpass output
• The summer output is the highpass: hp = x - 2R·lp - bp
• Highpass feeds integrator 1, producing bandpass: bp = ∫(ω_c·hp)
• Bandpass feeds integrator 2, producing lowpass: lp = ∫(ω_c·bp)

Step 2: Substitute trapezoidal integrators

Replace each analog integrator with a trapezoidal integrator. Each integrator has a state variable (s1 for bandpass, s2 for lowpass). The relationships become:

where g = tan(π·f_c/f_s). This is a system of three equations with three unknowns (hp, bp, lp). The trick is that we know x (the input) and the states s1, s2 (from the previous sample), so we can solve algebraically.

Step 3: Solve the algebraic loop

Substitute bp into the lp equation:

Substitute both into the hp equation:

Collect hp terms:

Define the denominator D = 1 + 2R·g + g²:

Now we can compute all outputs directly:

Step 4: Update the states

After computing outputs, update the integrator states for the next sample:

Equivalently: s1 += 2·g·hp, s2 += 2·g·bp. This form is sometimes more numerically stable.

Complete TPT SVF implementation

// TPT State Variable Filter - per-sample processing
struct TPTSVF {
    float s1 = 0;  // integrator 1 state (BP)
    float s2 = 0;  // integrator 2 state (LP)
    float g;       // tan(pi * fc / fs)
    float R;       // damping: R = 1/(2*Q)
    float D;       // precomputed denominator

    void setParams(float cutoffHz, float Q, float fs) {
        g = tanf(PI * cutoffHz / fs);
        R = 1.0f / (2.0f * Q);
        D = 1.0f / (1.0f + 2.0f*R*g + g*g);
    }

    void process(float x, float& lp, float& bp, float& hp) {
        // Solve the algebraic loop
        hp = (x - (2.0f*R + g)*s1 - s2) * D;

        // Compute outputs from integrators
        float v1 = g * hp;
        bp = v1 + s1;
        float v2 = g * bp;
        lp = v2 + s2;

        // Update states (equiv: s = 2*out - s)
        s1 = bp + v1;
        s2 = lp + v2;
    }
};

🌊 Resonance and stability in ZDF filters

The parameter R (or equivalently Q) controls resonance. In the TPT SVF:

• R = 1 (Q = 0.5): Critically damped—no resonant peak
• R = 0.5 (Q = 1): Mild resonance, common default
• R → 0 (Q → ∞): Self-oscillation—the filter becomes a sine oscillator

Unlike naive IIR filters, the TPT SVF remains stable even at R = 0. The self-oscillation is clean and tracks the cutoff frequency precisely. This is exactly what analog VCFs do—and what makes them musically useful as pitched elements.

📈 Adding nonlinearities: tanh and saturation

Analog filters aren't linear. Transistors and op-amps saturate, diodes conduct asymmetrically, and capacitors have voltage-dependent characteristics. These nonlinearities give analog filters their "warmth"—they compress dynamics, add harmonics, and prevent runaway oscillation at high resonance. To capture this character digitally, we need to add nonlinear elements—but carefully, to maintain stability.

Where to add nonlinearity

In a state variable filter, the natural places for nonlinearity are:

• Integrator inputs: Apply tanh() to the signal entering each integrator, modeling op-amp slew limiting
• Feedback path: Saturate the signal being fed back, preventing runaway at high resonance
• Output stage: Add soft clipping to the final output for overall warmth

The challenge: nonlinear algebraic loops

When you add nonlinearity inside the feedback loop, you can no longer solve for the output with simple algebra. The equation becomes implicit:

This requires iterative solving (Newton-Raphson) or careful approximation. Two practical approaches:

void processNonlinear(float x, float& lp) {
    float hp_est = 0, bp_est = 0;

    for (int i = 0; i < 4; i++) {
        float hp = (x - 2.0f * R * tanhf(lp) - bp_est) * D;
        float v1 = tanhf(g * hp);
        bp_est = v1 + s1;
        float v2 = tanhf(g * bp_est);
        lp = v2 + s2;
    }
    // Update states after convergence
    s1 = bp_est + tanhf(g * hp);
    s2 = lp + tanhf(g * bp_est);
}

void processLinearWithSaturation(float x, float& lp) {
    // Linear ZDF core
    float hp = (x - (2.0f*R + g)*s1 - s2) * D;
    float v1 = g * hp;
    float bp = v1 + s1;
    float v2 = g * bp;
    lp = v2 + s2;

    // Saturate outputs before state update
    float bp_sat = tanhf(bp);
    float lp_sat = tanhf(lp);

    // Update states with saturated values
    s1 = bp_sat + v1;
    s2 = lp_sat + v2;

    // Return saturated output
    lp = lp_sat;
}

Fast tanh approximations

The standard library tanh() is expensive. For real-time audio, use polynomial approximations:

// Fast tanh (Pade approximant, ~0.1% error)
inline float fastTanh(float x) {
    float x2 = x * x;
    return x * (27.0f + x2) / (27.0f + 9.0f*x2);
}

// Even faster: cubic soft clip
inline float softClip(float x) {
    if (x > 1.0f) return 1.0f;
    if (x < -1.0f) return -1.0f;
    return x - x*x*x / 3.0f;
}

🎹 Case study: TPT ladder filter

The classic Moog ladder filter is a cascade of four one-pole lowpass stages with global feedback for resonance. Let's build it using the TPT approach with nonlinear stages.

Ladder topology

Each stage is a one-pole TPT lowpass. The feedback is taken from the fourth stage output and subtracted from the input, scaled by the resonance factor k (where k = 4 gives self-oscillation in the linear case).

The algebraic loop solution

With four cascaded stages and feedback, we need to solve for the final output algebraically. For the linear case, define g = tan(π·fc/fs) and G = g/(1+g) for each stage. The closed-form solution is:

struct TPTLadder {
    float s[4] = {0};  // four integrator states
    float g;           // tan(pi * fc / fs)
    float k;           // resonance 0-4
    float G;           // g / (1 + g)

    void setParams(float cutoffHz, float res, float fs) {
        g = tanf(PI * cutoffHz / fs);
        G = g / (1.0f + g);
        k = res * 4.0f;  // 0-1 maps to 0-4
    }

    float process(float x) {
        // Feedback contribution from states
        float G2 = G * G;
        float G3 = G2 * G;
        float G4 = G3 * G;
        float S = G3*s[0] + G2*s[1] + G*s[2] + s[3];

        // Solve for output (linear case)
        float y4 = (G4*x + S) / (1.0f + k*G4);

        // Back-substitute: input after feedback
        float u = x - k * y4;

        // Process through stages
        float y = u;
        for (int i = 0; i < 4; i++) {
            float v = (y - s[i]) * G;
            y = v + s[i];
            s[i] = y + v;
        }
        return y;
    }
};

Adding tanh saturation to each stage

Real ladder filters have transistor saturation in each stage. Adding tanh to each stage creates an implicit nonlinear system that requires iteration:

float processNonlinear(float x) {
    // Initial estimate from previous output
    float y4_est = s[3];

    // Newton-Raphson iterations
    for (int iter = 0; iter < 3; iter++) {
        float u = tanhf(x - k * y4_est);

        float y = u;
        float s_temp[4];
        for (int i = 0; i < 4; i++) {
            float v = tanhf((y - s[i]) * G);
            y = v + s[i];
            s_temp[i] = y + v;
        }
        y4_est = y;
    }

    // Final pass with state update
    float u = tanhf(x - k * y4_est);
    float y = u;
    for (int i = 0; i < 4; i++) {
        float v = tanhf((y - s[i]) * G);
        y = v + s[i];
        s[i] = y + v;
    }

    return y;
}

⚙️ Implementation patterns

Per-sample vs. per-block coefficient updates

Computing g = tan(π·fc/fs) for every sample is expensive. Options:

// Fast tan approximation for per-sample updates
// Valid for |x| < pi/4 (fc < fs/4)
inline float fastTan(float x) {
    float x2 = x * x;
    return x * (1.0f + x2*(0.333333f + x2*0.1333333f));
}

// Lookup table with linear interpolation
// Precompute: tanTable[i] = tan(pi * i / 8192)
float tanTable[4096];

// normalizedFreq = fc/fs, range 0..0.5
float tableTan(float normalizedFreq) {
    float idx = normalizedFreq * 8192.0f;
    int i = (int)idx;
    float frac = idx - i;
    return tanTable[i] + frac*(tanTable[i+1] - tanTable[i]);
}

Oversampling considerations

ZDF filters are more stable than naive IIR at high frequencies, but nonlinearities still alias. Consider oversampling when:

• Heavy saturation (tanh in every stage)
• Self-oscillation at high frequencies (the sine wave needs bandwidth)
• Processing already-bright material through high resonance

2x oversampling often suffices. Use minimum-phase IIR decimation filters to preserve transient response. The oversampled section should include the entire nonlinear filter, not just the saturation functions.

Modulation response

ZDF filters handle modulation better than naive IIR, but there are still considerations:

For audio-rate modulation (filter FM), coefficient smoothing defeats the purpose. Instead:

• Update coefficients every sample (use fast tan approximation)
• Accept that very fast modulation will create sidebands (that's the point of filter FM)
• Keep the modulation depth reasonable; extreme sweeps still cause transients

🧪 Testing and validation

Digital filter implementations need rigorous testing. Here's a practical checklist:

// Simple stability test
void testStability(Filter& f) {
    float maxSeen = 0;

    for (int i = 0; i < 100000; i++) {
        // Random modulation
        float fc = 20.0f + rand01() * 19980.0f;
        float res = rand01();
        f.setParams(fc, res, 44100);

        // Process noise
        float x = rand01() * 2.0f - 1.0f;
        float y = f.process(x);

        maxSeen = fmaxf(maxSeen, fabsf(y));
        if (fabsf(y) > 10.0f) {
            printf("UNSTABLE at sample %d!\n", i);
            return;
        }
    }
    printf("Stable. Max output: %.2f\n", maxSeen);
}

⚠️ Common pitfalls and solutions

// Denormal prevention - add tiny DC offset
const float DENORMAL_GUARD = 1e-18f;

void process(float x) {
    // Add to input or states
    x += DENORMAL_GUARD;
    // ... filter processing ...
    s1 += DENORMAL_GUARD;
    s2 += DENORMAL_GUARD;
}

// Alternative: set FTZ at plugin init (x86)
#include <xmmintrin.h>
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);

📚 Further resources

For deeper understanding of ZDF and virtual analog filter design:

• The Art of VA Filter Design by Vadim Zavalishin — the definitive reference, covers TPT, ladder, and advanced nonlinear techniques
• DAFx Conference Proceedings — academic papers on virtual analog, many with source code
• Physical Audio Signal Processing by Julius O. Smith III — covers waveguide and physical modeling relevant to filter resonance
• KVR DSP Forum — active community discussing implementation details and sharing code
• Demystifying Synthesizer Filters — our companion article covering analog topologies and filter fundamentals

🎯 Summary

Zero-delay feedback filters solve the fundamental problem of digital filter resonance: the unit-sample delay in feedback that causes instability and frequency errors. By solving the filter equations algebraically within each sample, ZDF/TPT filters achieve:

• Stable self-oscillation that tracks the cutoff frequency accurately
• Correct phase response throughout the audio range
• Graceful behavior under fast modulation
• Natural integration of nonlinear elements (saturation, tanh) for analog character

The trade-off is complexity: you can't just plug in biquad coefficients from a cookbook. You need to understand the topology of the filter you're modeling and derive the algebraic solution for its feedback loops. For the TPT SVF and ladder filters presented here, the math is tractable and the results are worth it—filters that behave musically across their entire parameter range, just like their analog inspirations.