Phase response of a biquad

[fuo]

Hello,

This is my first post here.
I found that site and forum thanks to a paper published by Robin Schmidt here:
http://www.rs-met.com/misc.html#DSPNotes
namely "Basic Digital Filters": http://www.rs-met.com/documents/dsp/Bas ... ilters.pdf

In this note Robin lists transfer functions as well as real/imaginary parts and magnitude/phase response for different filters.

Unfortunately for me the one I was looking for is missing: I need the phase response of the biquad (should be just after the magnitude response (18)).
Real and imaginary functions would also do of course.

I suppose I should be able to calculate them from the complex transfer function, but I have to admit my mathematical skills are lacking

So any help would be more than welcome here

Thanks
Fuo

[thevinn]

DSPFilters (in my signature) calculates the phase response, you can crib the formula from there:

https://github.com/vinniefalco/DSPFiltersDemo

[Music Engineer]

it should be straightforward to derive the biquad case from the general case - just compute ca, sa, cb, cb with N = M = 2 and plug it into the last formula. beware that the formula, as it stands, does not yet include a phase-unwrapping, which - if i'm not mistaken - is possible only when you compute the phase for all frequencies (and then choose as unwrapped phase the computed value plus a suitable multiple of 2*pi, such that the plot becomes continuous)

[fuo]

DSPFilters (in my signature) calculates the phase response, you can crib the formula from there:

Thanks!
I will try digging into the source code.

Nice software by the way!

[fuo]

it should be straightforward to derive the biquad case from the general case - just compute ca, sa, cb, cb with N = M = 2 and plug it into the last formula. beware that the formula, as it stands, does not yet include a phase-unwrapping, which - if i'm not mistaken - is possible only when you compute the phase for all frequencies (and then choose as unwrapped phase the computed value plus a suitable multiple of 2*pi, such that the plot becomes continuous)

Thanks Robin, that looks straightforward enough indeed
First tests look promising so I hope I am on the right track now.
Thanks again!

[fuo]

Robin,
The implementation based on the general direct form, as per your directives, works like a charm!
It looks like the phase is already unwrapped when using atan2.

On the other hand, I also tried the magnitude formula from the biquad section, just for the sake of it, but failed to make it work...

Thank you!

[Music Engineer]

On the other hand, I also tried the magnitude formula from the biquad section, just for the sake of it, but failed to make it work.

hmm - it's a while ago and i'm not sure anymore whether i verified the formula numerically. i tried to re-derive the magnitude formula and it seems bogus indeed. now, i arrived at something else. could you try the following:

Code: Select all

g  = 1 / (ca^2 + sa^2);  // scale factor
Hr = g*(cb*ca + sb*sa);  // real part of H
Hi = g*(cb*sa - sb*ca);  // imaginary part of H
M  = sqrt(Hr^2 + Hi^2);  // magnitude

[fuo]

Yes that formula works

It gives the same result as the magnitude formula from the General Direct Form section :

Code: Select all

M = sqrt((cb^2+sb^2)/(ca^2+sa^2))
P = atan2(sb,cb)-atan2(sa,ca))

Here is the ca, sa, cb, sb calculation from biquad parameters (normalized for a0=1):

Code: Select all

ca = 1 + a1*cos(w) + a2*cos(2*w)
sa = a1*sin(w) + a2*sin(2*w)
cb = b0 + b1*cos(w) + b2*cos(2*w)
sb = b1*sin(w) + b2*sin(2*w)

[Music Engineer]

Yes that formula works

It gives the same result as the magnitude formula from the General Direct Form section :

Code: Select all

M = sqrt((cb^2+sb^2)/(ca^2+sa^2))
P = atan2(sb,cb)-atan2(sa,ca))

oh - really? it was actually that magnitude formula (for the general case) which seemed bogus to me now. for the formula in the pdf, i did (i think so at leat):

Code: Select all

                          cb - j*sb     cb + j*sb
H^2(w) = H(w) * H'(w) =  ----------- * -----------
                          ca - j*sa     ca + j*sa

where H' should be the complex conjugate of H. i wasn't sure anymore, if i can obtain the complex conjugate of H by simply using the complex conjugates inside computation of H (which was apparently what i did back then). so what i did today, was to actually split H into its real and imaginary parts (by multiplying numerator and denominator by the complex conjugate of the denominator).

but you were refering to the biquad specific magnitude formula? that's actually a formula that i have in practical use and it never gave me any problems. this is how my implementation looks like:

Code: Select all

double biquadMagnitudeAt(double b0, double b1, double b2, double a1, double a2, double w)
{
  double c1  = 2*cos(w);
  double c2  = 2*cos(2*w);
  double num = b0*b0 + b1*b1 + b2*b2 + c1*(b0*b1 + b1*b2) + c2*b0*b2;
  double den = 1     + a1*a1 + a2*a2 + c1*(   a1 + a1*a2) + c2*   a2;
  double p   = num/den;
  return sqrt(p);
}

(sorry for the messed up formatting)

[Music Engineer]

i wasn't sure anymore, if i can obtain the complex conjugate of H by simply using the complex conjugates inside computation of H

OK - i just verified it. it works. we generally have (z/w)' = z' / w'. so both formulas should be equivalent.

that wasn't obvious to me right now and wikipedia doesn't list that identity either:
http://en.wikipedia.org/wiki/Complex_conjugate

[Music Engineer]

however, this "new" formula doesn't seem to be entirely useless because, if i'm not mistaken, it should now be possible to compute the phase as:

P = atan2(Hi, Hr)

which gets rid of one of the two calls to atan2.

as for the phase (un)wrapping - it's perhaps not relevant for your specific cases because the phase response remains inside the wrapping-interval anyway. for other (especially higher order) filters, it may exceed the wrapping interval and hence be "wrapped" into it.

[fuo]

Strange, it is indeed formula ( 18 ) that does not work for me...
Anyway, I now have a working solution thanks to your help

as for the phase (un)wrapping - it's perhaps not relevant for your specific cases because the phase response remains inside the wrapping-interval anyway. for other (especially higher order) filters, it may exceed the wrapping interval and hence be "wrapped" into it.

Yes, stupid me...
I summed several biquads to exceed 90° and could not see any wrapping... of course as it would be a per biquad wrapping...