done. and - lo and behold - it worked on first go. well, at least, plotting the impulse responses looked reasonable and LP+2*R*BP+HP reconstruct the original signal, as it should be. here it is:Robin from www.rs-met.com wrote: i think, now i'll finally get my hands dirty with implementing the TPT/ZDF-SVF (which is also the first time, that i implement an SVF at all - as strange as that may sound)
Code: Select all
// Example code for a TPT/ZDF-SVF (obviously, this could be optimized for production code)
// written by Robin Schmidt, licensing/copyright: none (public domain)
// Parameters:
// fc: cutoff frequency
// fs: sample-rate
// R: damping coefficient
// N: number of samples
// x: input signal
// yL: lowpass output
// yB: bandpass output
// yH: highpass output
void svf(double fc, double fs, double R, int N, double x[], double yL[], double yB[], double yH[])
{
double wd = 2*PI*fc; // target radian frequency
double T = 1/fs; // sampling period
double wa = (2/T)*tan(wd*T/2); // prewarped radian frequency for analog filter (Eq. 3.7)
double g = wa*T/2; // embedded integrator gain (Fig 3.11), wc == wa
// states of the 2 integrators, static so we can call the function block-wise while maintaining
// the states from block to block:
static double s1 = 0.0;
static double s2 = 0.0;
// loop over the samples:
for(int n = 0; n < N; n++)
{
// compute highpass output via Eq. 5.1:
yH[n] = (x[n] - 2*R*s1 - g*s1 - s2) / (1 + 2*R*g + g*g);
// compute bandpass output by applying 1st integrator to highpass output:
yB[n] = g*yH[n] + s1;
s1 = g*yH[n] + yB[n]; // state update in 1st integrator
// compute lowpass output by applying 2nd integrator to bandpass output:
yL[n] = g*yB[n] + s2;
s2 = g*yB[n] + yL[n]; // state update in 2nd integrator
// Remark: we have used two TDF2 integrators (Fig. 3.11) where one of them would be in code:
// y = g*x + s; // output computation
// s = g*x + y; // state update
}
}
yeah - so, as we see, the feedback is indeed crucial for this DC reconstruction. i guess, i'm already way too spoiled with this strictly sequential thinking in DSP.