2
$\begingroup$

Using formulae from Audio EQ Cookbook from here it was easy to implement biquad IIR filters in C++, i.e. to calculate coefficients b0, b1, b2, a0, a1, a2.

Now, I want to calculate the phase response for a given IIR filter. I have read a number of papers about IIR design, read the websites (DSP Guide) but could not find the way to do it. How I can calculate and specify the phase of an IIR filter during the design?

Improve this question
$\endgroup$
3
  • $\begingroup$ You are asking two different questions here. Computing the phase response of a given filter is straightforward: just take the argument of the complex frequency response. Matlab can help you. The other question how to design a filter which approximates a given desired phase response is much more difficult. You can always use some optimization routine but the problem is to find a reasonable phase specification which can realistically be satisfied by an IIR filter of a given order. $\endgroup$ – Matt L. Jun 18 '15 at 13:11
  • $\begingroup$ Thanks Matt. I guess I am looking to computing phase response of a given filter. I have to do it without MATLAB, i.e. to implement it in C++. Is it possible? Any guide or example how to do it? I assume I can use transfer function for biquad from the Audio Cookbook and try to find the argument. Is that the right approach? Thank you. $\endgroup$ – Nebojsa Jun 18 '15 at 13:16
  • $\begingroup$ I've posted an answer, hope it helps. $\endgroup$ – Matt L. Jun 18 '15 at 13:25
4
$\begingroup$

Finding the phase response of a biquad at a specific frequency is simple. Recall the transfer function of a biquad:

$$ H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}} $$

The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where $\omega$ is a normalized frequency in the range $[-\pi, \pi)$. SO, it would look like this:

$$ H(e^{j\omega}) = \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}} $$

Because of the complex exponentials, the value of $H(e^{j\omega})$ will be complex. The phase response at the frequency $\omega$ is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

The only other detail you might need is how to arrive at $\omega$: given a signal sampled at sample rate $f_s$ Hz, if you want to know the frequency response at a given frequency $f$ Hz, you can use the above equation, and let:

$$ \omega = \frac{2 \pi f}{f_s} $$

Improve this answer
$\endgroup$
8
  • $\begingroup$ It should be $e^{-j\omega}$ in the formula for the frequency response $H(\omega)$. And note that in the cookbook formulas usually $a_0\neq 1$. $\endgroup$ – Matt L. Jun 18 '15 at 13:25
  • $\begingroup$ Thanks for pointing out the missing signs. I haven't read through the cookbook, so I'll take your word for that. I find it universally easier to normalize the coefficients so that $a_0 = 1$. $\endgroup$ – Jason R Jun 18 '15 at 13:30
  • $\begingroup$ I agree that normalizing such that $a_0=1$ makes a lot of sense. But unfortunately it's not always done, and it will cause errors in implementations if $a_0$ is considered to be $1$ if it isn't. $\endgroup$ – Matt L. Jun 18 '15 at 13:33
  • $\begingroup$ @JasonR, you need not read through any cookbook to know that when $e^{j\omega}$ is substituted for $z$ in $z^{-1}$, the result is $e^{-j\omega}$. the only thing i might add to one of the answers is that if this is implemented in a limited precision context (fixed-point or single-precision floating-point), my experience is that with either magnitude or phase response is that you run into the "cosine problem" when $|\omega| \ll 1$. it's because $$ \cos(\omega) = 1 - \sin^2 \left( \frac{\omega}{2} \right) $$. all of the information regarding $\omega$ is in the difference from 1. $\endgroup$ – robert bristow-johnson Jun 18 '15 at 15:44
  • $\begingroup$ so even with floating-point, all of the significant bits of $ \sin^2 \left( \frac{\omega}{2} \right) $ get added to 1 (which becomes the "hidden 1 bit" in IEEE floating point) and when $|\omega| \ll 1$ most of those significant bits fall offa the edge. so, if you want to calculate and display the phase response (or, better yet, either the group or phase delay, i think that would be more meaningful to look at) in an embedded context, i would replace all instances of $\cos(\cdot)$ with the above identity and re-jigger your equations. otherwise it will look like shit for small $\omega$. $\endgroup$ – robert bristow-johnson Jun 18 '15 at 15:51
3
$\begingroup$

You need the argument of the filter's frequency response. The latter is given by

$$H(\omega)=\frac{b_0+b_1e^{-j\omega}+b_2e^{-2j\omega}}{a_0+a_1e^{-j\omega}+a_2e^{-2j\omega}}\tag{1}$$

This is a complex function and its argument can be computed by using the function atan2(y,x), where y is the imaginary part of $H(\omega)$ and x is the real part of $H(\omega)$.

Improve this answer
$\endgroup$
0
$\begingroup$

If you factor the transfer function of an IIR into the locations of its poles and zeros, you can graphically solve (or even visually inspect for a small number of poles/zeros) the phase response. Walk around the unit circle in the complex plane, and add/subtract the angles from all the transfer functions poles and zeros to any point on the unit circle for the phase response at the frequency represented by that point.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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