Phase response of a biquad

DSP, Plugin and Host development discussion.
Post Reply New Topic
RELATED
PRODUCTS

Post

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

Post

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

https://github.com/vinniefalco/DSPFiltersDemo

Image

Post

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)
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

thevinn wrote: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!

Post

Robin from www.rs-met.com wrote: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 :D
First tests look promising so I hope I am on the right track now.
Thanks again!

Post

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!

Post

fuo wrote: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
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

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)

Post

fuo wrote: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)
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Robin from www.rs-met.com wrote: 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
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

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.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Strange, it is indeed formula ( 18 ) that does not work for me...
Anyway, I now have a working solution thanks to your help :)
Robin from www.rs-met.com wrote: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...

Post Reply

Return to “DSP and Plugin Development”