DSP filter. Is prewarping performed when discretizing using Forward/Backward?

Question

I am trying to derive the coefficients used for a IIR implementation for the lowpass portion of a SVF filter. I've seen finished expressions for the coefficients (smith (p.8) and victor), their discretizized transfer function looks something like (not super important for my question)

$$ H(z) = \frac{k_f^2 \, z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} \, , $$ with $k_q = \frac{1}{Q}$ and (next is central to my question) $k_f = 2 \, \sin (\pi f_c T) $ .

I would like to derive these parameters myself, which I have attempted, but I am having trouble with the parameters involving the cutoff frequency $f_c$!

This is my starting point, (from analysing SVF circuit), I get the following transfer function for the lowpass portion,

$$ H(s) = \frac{\Omega_0^2}{s^2 + \frac{\Omega_0}{Q}s + \Omega_0^2} \, . $$

Next step is performing the discretization. I did this using the bilinear transform ($z = e^{sT}$), using from what I gather the Tustin/trapezoidal approximation ($s = \frac{1}{T} \log(z) \approx \frac{2}{T} \frac{z-1}{z+1})$.

Prewarping

This transform requires prewarping and was done as follows,

$$ \Omega_0 = \tfrac{2}{T} \tan \left( \tfrac{\omega_0}{2} \right) \,, $$

where $\Omega_0 = 2 \pi f_c$ and $f_c$ is my wanted cutoff frequency.

Now, already I can see that my derivation is going to result in a different expression than what the two links above are getting. Their frequency parameters shows up in expressions that look like $ f = 2 \sin(\pi f_c T) $

My Question

Basically, where is the $\sin()$ term coming from?

Could it be that they use Euler forward (or backward) as an approximation?

$$ s = \frac{1}{T} \log(z) \approx \frac{z - 1}{T} $$

And then done some prewarping for that case? I have looked around the web and can not find any mention of prewarping done for euler backward/froward. Is that usually done?

That would give me $$ \begin{split} s = j \Omega &= \frac{z-1}{T} \\ j \Omega &= \frac{e^{j \omega}-1}{T} \\ j \Omega &= \frac{e^{j \omega/2}-e^{-j \omega/2}}{T e^{-j \omega/2}} \\ j \Omega &= \frac{j 2 \sin(\omega/2)}{T e^{-j \omega /2}} \\ \Omega &= \frac{2 \sin(\omega/2)}{T e^{-j \omega /2}} \\ \end{split} $$

Not sure what I am missing.

--

Thank you very much for any help regarding this!

Answer 1

(After some helpful comments, I think I see it a bit more clearly. Here is an attempt at answering my own question, please let me know if I am completely off track.)

Answer:

In my question I must have misunderstood from where the top expression $H(z)$ comes from. I thought that such an expression was the result of a direct derivation from an analog circuit using the bilinear transform.

As Robert comments to the original question, $H(z)$ as I've written it, is the transfer function for the Hal Chamberlin SVF. And is not the result of a direct derivation from an analog transfer function using the bilinear transform (also pointed out by Matt).

Also, to answer my other question as far as I know there is no pre-warping done for Euler forward/backward and the term $\sin(\pi f_c T)$ does not come up because of this.

Where does the $\sin(\pi f_c T)$ term come from?

Bluntly, I think from an approximation. An approximation that is baked into to what we call the Hal Chamberlin SVF (presented in "Musical Applications of Microprocessors").

By using Roberts fantastic EQ cookbook, which has transfer functions and expressions for the coefficients, but which are derived from analog prototypes, I can compare and see where the approximation is being made.

(I received a nice guide that walks through the Hal Chamberlin SVF its transfer function and how it compares to transfer functions derived using the Bilinear transform and later the Impulse Invariant mapping. Link here: "Exact Bilinear Transform or Impulse Invariant mappings of Q and resonant frequency from analog prototype to Hal Chamberlinʼs State Variable Filter (SVF)" by Robert Bristow-Johnson. I'll also paraphrase the part which explains the $\sin{()}$ term.)

Starting from (note: not derived via bilinear transform)

$$ H(z) = \frac{k_f^2 \, z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} \, , $$

I can compare the coefficients in the denominator with the coefficients that were derived from an analog lowpass biquad which are presented in the EQ cookbook. From this we see that

$$ \frac{a_2}{a_0} = \frac{1 - \frac{\sin(\omega_c)}{2 Q}}{1+\frac{\sin(\omega_c)}{2 Q}} = 1 - k_f k_q $$ $$ \frac{a_1}{a_0} = \frac{ - 2 \cos(\omega_c)}{1+\frac{\sin(\omega_c)}{2 Q}} = k_f^2 + k_f k_q - 2 $$

where $\omega_c = 2 \pi f_c T$.

Using these two and solving for $k_f$ gives us,

$$ k_f = \frac{2 \sin(\frac{\omega_c}{2})}{\sqrt{1+\frac{\sin(\omega_c)}{2 Q}}} \, , $$

If we assume $Q > \frac{1}{2}$ and that $\omega_c \ll \pi$ (cutoff frequency is small compared to Nyquist) we can approximate,

$$ k_f \approx 2 \sin\left(\frac{\omega_c}{2} \right) = 2 \sin(\pi f_c T) \, , $$

and we have the $\sin(\cdot)$ term I was looking for.