Login / Register  0 items | $0.00New @ What is KVR? Submit News Advertise

Phase response of a biquad

DSP, Plug-in and Host development discussion.

Moderator: Moderators (Main)

fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Fri Jan 25, 2013 8:16 am Phase response of a biquad

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 (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
KVRian
 
771 posts since 30 Nov, 2008

Postby thevinn; Fri Jan 25, 2013 8:30 am

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

https://github.com/vinniefalco/DSPFiltersDemo

Image
Music Engineer
KVRAF
 
3651 posts since 8 Mar, 2004, from Berlin, Germany

Postby Music Engineer; Fri Jan 25, 2013 8:37 am

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)
Image
fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Fri Jan 25, 2013 10:15 am

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!
fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Fri Jan 25, 2013 10:18 am

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!
fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Mon Jan 28, 2013 4:04 am

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
KVRAF
 
3651 posts since 8 Mar, 2004, from Berlin, Germany

Postby Music Engineer; Mon Jan 28, 2013 5:27 am

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
Image
fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Mon Jan 28, 2013 6:33 am

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
KVRAF
 
3651 posts since 8 Mar, 2004, from Berlin, Germany

Postby Music Engineer; Mon Jan 28, 2013 6:50 am

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)
Image
Music Engineer
KVRAF
 
3651 posts since 8 Mar, 2004, from Berlin, Germany

Postby Music Engineer; Mon Jan 28, 2013 7:07 am

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
Image
Music Engineer
KVRAF
 
3651 posts since 8 Mar, 2004, from Berlin, Germany

Postby Music Engineer; Mon Jan 28, 2013 7:18 am

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.
Image
fuo
KVRer
 
6 posts since 25 Jan, 2013

Postby fuo; Mon Jan 28, 2013 8:09 am

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

Moderator: Moderators (Main)

Return to DSP and Plug-in Development