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/BasicDigitalFilters.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 (1).
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

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

https://github.com/vinniefalco/DSPFiltersDemo

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!
I will try digging into the source code.

Nice software by the way!

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

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!

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:

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

Yes that formula works

It gives the same result as the magnitude formula from the General Direct Form section :
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):

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)

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):

                          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:

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)

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

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.

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


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