What is KVR Audio? | Submit News | Advertise | Developer Account

Options (Affects News & Product results only):

OS:
Format:
Include:
Quick Search KVR

"Quick Search" KVR Audio's Product Database, News Items, Developer Listings, Forum Topics and videos here. For advanced Product Database searching please use the full product search. For the forum you can use the phpBB forum search.

To utilize the power of Google you can use the integrated Google Site Search.

Products 0

Developers 0

News 0

Forum 0

Videos 0

Search  

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
 
767 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
Robin from www.rs-met.com
KVRAF
 
3648 posts since 8 Mar, 2004, from Berlin, Germany

Postby Robin from www.rs-met.com; 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!
Robin from www.rs-met.com
KVRAF
 
3648 posts since 8 Mar, 2004, from Berlin, Germany

Postby Robin from www.rs-met.com; 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)
Robin from www.rs-met.com
KVRAF
 
3648 posts since 8 Mar, 2004, from Berlin, Germany

Postby Robin from www.rs-met.com; 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
Robin from www.rs-met.com
KVRAF
 
3648 posts since 8 Mar, 2004, from Berlin, Germany

Postby Robin from www.rs-met.com; 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
Robin from www.rs-met.com
KVRAF
 
3648 posts since 8 Mar, 2004, from Berlin, Germany

Postby Robin from www.rs-met.com; 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