How do I use the Lazzarini-Timoney DSVF as a sine generator?

Question

I've tested out the Chamberlin digital state variable filter (DSVF) as a sine generator as per this well-known approach. At $f_0=686\text{Hz}$ tone and $f_s=48000\text{Hz}$ it sounds fine and the FFT shows exactly two peaks at -220Hz and 220Hz; at $f_0=687\text{Hz}$ it buzzes.

The Chamberlin filter looks as below. Setting $Q=\infty$ makes $-q=0$ and deletes the bandpass feedback. From there the bandpass is cosine and the lowpass is sine.

The Chamberlin filter is known to be unstable above a certain fraction of the sample rate (generally $\dfrac{f_s}{6}$); apparently as an oscillator it's a bit more extreme. Lazzarini and Timoney proposed a slight adjustment that gives the filter exactly the same amplitude response as a biquadratic filter. I'd like to try that as a sine wave generator, but I have no idea what I'm doing!

That filter looks like this:

I encourage you to read the paper, it's pretty cool.

My attempts have produced noise.

This isn't really super-critical for my project. I'm using it as a cheap LFO generator and it's stable up to 500Hz; I'm just curious if it can be improved.

Edit: I have now posted code because I can't figure out why the filter implementation is broken.

Answer 1

I liked that paper, too, when I first saw it. It does add to some of the confusion in the AES about the definition of Q. I really think that we need to nail things down to a standard definition that makes sense for audio people.

Also, when you use any second-order IIR filter for a sinusoid generator, it eventually boils down to an input (what you're calling "$s_n$") being 0 and setting the two initial states to exactly what you need to get the amplitude and phase to be what you want. Ultimately the math comes out to be

$$ x[n] = 2 \cos(\omega_0) x[n-1] - x[n-2] $$

$$ x[n] - 2 \cos(\omega_0) x[n-1] + x[n-2] = 0 $$

$$ X(z) - 2 \cos(\omega_0) X(z)z^{-1} + X(z)z^{-2} = 0 $$

$$ X(z)z^{2} - 2 \cos(\omega_0) X(z)z + X(z) = 0 $$

$$ X(z)(z^{2} - 2 \cos(\omega_0)z + 1) = 0 $$

$$ z^{2} - 2 \cos(\omega_0)z + 1 = 0 $$

$$ z_\text{p} = \begin{cases} \cos(\omega_0) \pm \sqrt{\cos^2(\omega_0)-1} \quad & \cos^2(\omega_0) \ge 1 \\ \\ \cos(\omega_0) \pm j\sqrt{1-\cos^2(\omega_0)} \\ =\cos(\omega_0) \pm j\sin(\omega_0) \\ = e^{ \pm j\omega_0} \quad & \cos^2(\omega_0) \le 1 \\ \end{cases}$$

But only the second case counts for real $\omega_0$.

It doesn't matter what form of filter you make, it boils down to that and, if your word size is big enough (like IEEE float or IEEE double), that's just the simplest way to generate a sinusoid.

At one time I worried about this being stable for long enough to be useful (because there are two poles directly on the unit circle $|z|=1$, so this is "critically stable" or "marginally stable", not fully stable) but James McCartney set me straight on that some 20+ years ago.

Now, this inherent stability is not necessarily the case if you want to make a quadrature sine and cosine wave generator:

$$ \begin{align} x[n] &= \cos(\omega_0)x[n-1] - \sin(\omega_0)y[n-1] \\ y[n] &= \sin(\omega_0)x[n-1] + \cos(\omega_0)y[n-1] \\ \end{align} $$

Then you need to make sure that your $\cos(\omega_0)$ and $\sin(\omega_0)$ coefficients square and add to 1. If they add to something less than 1, your cosine and sine waves will slowly die out and if they add to more than 1, the waves will slowly grow indefinitely. So that quadrature generator needs some adaptive feedback to scale this coefficients based on the amplitude of the complex sinusoidal wave.