0

Using formula from Audio EQ Cookbook I implemented biquad IIR filters in Excel.

I can now plot every transfer function and total of 8 band EQ with shelvings and peakings.

My problem is, I can not implement any of the phase response formula in excel. I'm not very good at these kinds of functions and especially implementing complex number functions in excel.

Here is the transfer function.

enter image description here

I have $\omega_0, \alpha, A_x, \omega$ of the frequency, $\phi$ of the frequency $( = 4 * \sin{(\omega/2)}^2$, I don't know what this is), $b_0,b_1,b_2,a_1,a_2$ coefficients of every filter and dB values at each frequency ($H(s)$).

How can I write the phase response function in excel?

lennon310
  • 3,590
  • 19
  • 24
  • 27
BugraKezan
  • 87
  • 7

1 Answers1

2

Excel is an awkward language for this type of thing so we can make easier if we do a little math upfront. Doing complex math in Excel is just about as much fun as slamming your fingers in the car door, so we shall try to avoid it.

We have

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

If we want to evaluate this at normalized frequency $z = e^{j\omega}$ we get

$$B(z) = b_0 + b_1\cos(\omega)+b_2\cos(2\omega)-j\left[b_1\sin(\omega)+b_2\sin(2\omega)\right] = B_r+jB_i$$ with $$B_r = b_0 + b_1\cos(\omega)+b_2\cos(2\omega) \\ B_i = -b_1\sin(\omega)-b_2\sin(2\omega) $$

So the whole thing becomes

$$H(\omega) = \frac{B_r+jB_i}{A_r+jA_i}$$

With this we can calculate both magnitude and phase

$$|H(\omega)| = \sqrt{\frac{B_r^2+B_i^2}{A_r^2+A_i^2}} $$

$$\angle{H(\omega)} = \tan^{-1}\left(\frac{B_i}{B_r}\right)-\tan^{-1}\left(\frac{A_i}{A_r}\right)$$

We can implement this in excel by first calculating $B_r$, $B_i$, $A_r$, and $A_i$ and then using the last two formulas to calculate magnitude and phase. The phase can be done using the Excel function atan2() which in typical Excel fashion uses non-standard conventions and reverses the argument order.

You can find an example with a second order Butterworth lowpass at https://docs.google.com/spreadsheets/d/1ye7-WzmwNRpTglJsGk6Dtvb8kge9PICc/edit?usp=sharing&ouid=110367350423053204811&rtpof=true&sd=true

EDIT based on comments

It's been rightfully pointed out that the subtraction of the two phases can run into numerical problems and wrapping problems.

A better way to calculate the phase is the following one

$$H(\omega) = \frac{B_r+jB_i}{A_r+jA_i} = \frac{B_r+jB_i}{A_r+jA_i} \cdot \frac{A_r-jA_i}{A_r-jA_i} = \cdot \frac{B_rA_r+B_iA_i+j(B_iA_r-B_rA_i)}{X}, x \in \mathbb{R}$$

So we can get the phase of a single biquad as

$$\angle{H(\omega)} = \text{atan2}(B_iA_r-B_rA_i, B_rA_r+B_iA_i) $$ where $\text{atan2}()$ is the quadrant correct version of the inverse tangent. The spreadsheet has an added column for the new phase calculation.

Hilmar
  • 44,604
  • 1
  • 32
  • 63
  • That's amazing, you're very helpful. It works. Thanks a lot. I was trying to do this with =COMPLEX(), IMEXP() and =IMARGUMENT() functions but i guess more things get wrong. Just one thing, there is a very little shift at high frequencies, maybe about 0,2°. It's not much important but i wonder why? Image link: [link]https://im.ge/i/1cTVFG – BugraKezan Sep 30 '22 at 14:12
  • Hard to tell. Could be numerical: Sigma Studio generally works for fixed point processors. Could also be coefficient quantization (number of digits). If you want to know which one is correct, you need to use a known-good reference (Matlab, Octave or Python perhaps) or use filter where you can calculate the result by hand. – Hilmar Sep 30 '22 at 14:45
  • @BugraKezan if this helped you please also check this off as the right answer to close the open question, thank you! – Dan Boschen Sep 30 '22 at 16:10
  • 1
    I know the two arctans are split up, but I am concerned that there still will be a disparity between $\angle{H(\omega)} = \tan^{-1}\left(\frac{B_i}{B_r}\right)-\tan^{-1}\left(\frac{A_i}{A_r}\right)$ and the general complex argument. Also, if the OP wants the phase unwrapped, we should be computing the delta of phase between adjacent points of the phase response. – robert bristow-johnson Sep 30 '22 at 19:53
  • Not sure I understand your concern. As long as you use atan2() to calculate the angle, this should work fine. Unwrapping should not be necessary for a single biquad, I think. – Hilmar Sep 30 '22 at 21:41
  • Hello @Hilmar, I've tried an all-pass filter for phase response, but result is corrupted after -360 degree. As an example, a0 a1 a2 b0 b1 b2 omega/rad = 1,001022471 -1,998929175 0,998977529 0,998977529 -1,998929175 1,001022471 0,081812309

    IT2s 1kHz Q=16 All pass filter. It's being weird after 1750Hz. How can I resolve this?

     – BugraKezan Oct 06 '22 at 17:35
  • @robertbristow-johnson ; how can I unwrap the phase? – BugraKezan Oct 06 '22 at 17:57
  • 1
    @BugraKezan , please take a look at this answer to a related question. – robert bristow-johnson Oct 09 '22 at 21:40