can I calculate freq response generically from filters?

Question

I'd like to display the freq response of a Sallen Key filter type (in the specific case, this simple resonant filter).

on biquad, for example, I calculate magnitude and phase from zeros and poles, such as (in cpp):

// zeros
std::complex<double> resZeros(0.0, 0.0);
resZeros += mCoefficients.mA0 * std::exp(std::complex<double>(0.0, -0 * freq * k2PI));
resZeros += mCoefficients.mA1 * std::exp(std::complex<double>(0.0, -1 * freq * k2PI));
resZeros += mCoefficients.mA2 * std::exp(std::complex<double>(0.0, -2 * freq * k2PI));
// poles
std::complex<double> resPoles(0.0, 0.0);
resPoles += 1.0 * std::exp(std::complex<double>(0.0, -0 * freq * k2PI));
resPoles += mCoefficients.mB1 * std::exp(std::complex<double>(0.0, -1 * freq * k2PI));
resPoles += mCoefficients.mB2 * std::exp(std::complex<double>(0.0, -2 * freq * k2PI));
std::complex<double> res = resZeros / resPoles;
// magnitude
float magnitude = (float)std::abs(res);
if (std::isnan(magnitude)) {
    magnitude = 0.0f;
}
// phase
float phase = (float)std::arg(res); // i.e. math.atan2(res.im, res.re)
if (std::isnan(phase)) {
    phase = 0.0f;
}
return std::polar(magnitude, phase);

this code is taken from the web :) can I do the same with the filter above? generally, can I do it with different kind of filters?

For magnitude response of a biquad (which is one thing you're doing), you can sorta rearrange some of the mathematics to make it numerically a little better behaved. For phase, there is something else to work out, but what you should do is work out the phase change from the previous plotted point to the current and simply accumulate that phase change. That will automatically calculate group delay and nicely unwrap the phase. — robert bristow-johnson, Jul 03 '23 at 16:31
@robertbristow-johnson but I need response of this https://www.musicdsp.org/en/latest/Filters/29-resonant-filter.html, not the biquad (which I already have the code). — markzzz, Jul 03 '23 at 17:24
That example has two states (´buf0´ and ´buf1´) and it's LTI. So it's a biquad. With the exception of ´in´, attach a $z^{-1}$ to every variable and solve for the transfer function. Then you'll know how to apply the biquad frequency response. — robert bristow-johnson, Jul 03 '23 at 18:02
@robertbristow-johnson so its not a Sallen Key filter type in your opinion? — markzzz, Jul 03 '23 at 18:11
Sallen-Key is an analog circuit for second-order alalog filter. It gives you a "biquad" $H(s)$. But, even if you emulate the analog filter, you will end up with an $H(z)$ that may also be second order, and that's what I usually mean by a biquad and that's all your code can run. But you're right, you can run the frequency response of an analog filter on a computer, and your code looks like it. With $j\omega$ instead of $e^{j\omega}$. — robert bristow-johnson, Jul 03 '23 at 20:55
But that resonant filter from Paul Kellet that you pointed to, that's a digital filter with an $H(z)$. buf0 and buf1 get defined in terms of their previous values. — robert bristow-johnson, Jul 03 '23 at 20:59
@robert bristow-johnson: I'll do per step :) First of all, how would you convert the normalized freq used on that filter (0 to 1) to hz? Can't find a correct map function :O It seems 1000 hz is 0.14? Mmm... any clue? Can be evaluated? — markzzz, Jul 04 '23 at 19:29
@robert bristow-johnson maybe make sense open a separated question for this? — markzzz, Jul 04 '23 at 20:02
Sure, put that Paul Kellet filter in a question. I'll work it out. — robert bristow-johnson, Jul 04 '23 at 20:05

Dan Boschen · Answer 1 · 2023-07-03T17:28:42.313

To calculate the frequency response from the poles and zeros in continuous time (Laplace), write out the transfer function, $H(s)$, as the numerator (zeros) divided by the denominator (poles) and simply replace $s$ with $j\omega$, since the frequency response is the transfer function with $s$ restricted to be the $j\omega$ axis:

$$H(s) = \frac{(s-z_1)(s-z_2)(s-z_3)\ldots}{(s-p_1)(s-p_2)(s-p_3)\ldots} $$

With $z_1$, $z_2$, etc referring to the zero locations on the complex s-plane, and $p_1$, $p_2$, etc referring to the poles.

$$H(j\omega) = \frac{(j\omega-z_1)(j\omega-z_2)(j\omega-z_3)\ldots}{(j\omega-p_1)(j\omega-p_2)(j\omega-p_3)\ldots} $$

This is plug and chug with complex numbers, so multiply it out and then calclulate the magnitude and phase as a function of frequency, which is the frequency response.

For discrete time (Z-transform), it is the same process, except the frequency response is the transfer function in z, $H(z)$, with $z$ replaced with $e^{j\omega}$ with $\omega$ extending from $0$ to $2 \pi$ (the unit circle):

$$H(z) = \frac{(z-z_1)(z-z_2)(z-z_3)\ldots}{(z-p_1)(z-p_2)(z-p_3)\ldots} $$

$$H(e^{j\omega}) = \frac{(e^{j\omega}-z_1)(e^{j\omega}-z_2)(e^{j\omega}-z_3)\ldots}{(e^{j\omega}-p_1)(e^{j\omega}-p_2)(e^{j\omega}-p_3)\ldots} $$

Here's a very simple and graphical demonstration for a first order single pole continuous time RC low pass filter:

The RC filter results in a single real pole in the left half plane at $s=-\frac{1}{RC}$.

The transfer function for this is: $$H(s) = \frac{1}{s-p_1} = \frac{1}{s+1/(RC)}$$

Graphically on the s-plane we have a ratio of phasors where the numerator is always $1/angle{0}$ and the denominator $s-p_1$ with $s=j\omega$ starts off at DC when \omega=0 to be a phasor with magnitude $1/(RC)$ and angle $0$. This is our starting magnitude referenced by 0 dB on the Bode plot showing the relative magnitude and phase over frequency. As we increase $\omega$ up to the magnitude of the pole (when $s=j (1/(RC))$), the length of the denominator phasor will be $\sqrt{2}$ longer and at an angle of 45°. Given this is in the denominator and the numerator in this case was always 1, this is consistent with a -3 dB magnitude and -45° and is noted as the cutoff frequency. As we continue to sweep the frequency toward $\infty$ we see that the magnitude of the denominator approaches $\infty$ and the phase approaches +90°, consistent with the magnitude approaching $0$ and phase approaching -90° for $H(j\omega)$.

What I would do is replace the symbol for zeros from $z_n$ to $q_n$ so that when you go to digital and $H(z)$, you don't have symbol confusion. — robert bristow-johnson, Jul 04 '23 at 19:33
@robertbristow-johnson Yes thank you! That has always bothered me but I never had a good idea of how to fix that. — Dan Boschen, Jul 04 '23 at 19:45

can I calculate freq response generically from filters?

1 Answers1