0
$\begingroup$

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?

$\endgroup$
10
  • $\begingroup$ 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. $\endgroup$ Commented Jul 3, 2023 at 16:31
  • $\begingroup$ @robertbristow-johnson but I need response of this musicdsp.org/en/latest/Filters/29-resonant-filter.html, not the biquad (which I already have the code). $\endgroup$
    – markzzz
    Commented Jul 3, 2023 at 17:24
  • $\begingroup$ 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. $\endgroup$ Commented Jul 3, 2023 at 18:02
  • 1
    $\begingroup$ It can be evaluated if I work on it a little. $\endgroup$ Commented Jul 4, 2023 at 19:31
  • 1
    $\begingroup$ Sure, put that Paul Kellet filter in a question. I'll work it out. $\endgroup$ Commented Jul 4, 2023 at 20:05

1 Answer 1

1
$\begingroup$

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)$.

Filter Example

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Commented Jul 4, 2023 at 19:33
  • $\begingroup$ i just might modify your answer a little. $\endgroup$ Commented Jul 4, 2023 at 19:37
  • $\begingroup$ @robertbristow-johnson Yes thank you! That has always bothered me but I never had a good idea of how to fix that. $\endgroup$ Commented Jul 4, 2023 at 19:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.