Options (Affects News & Product results only):

 OS: Any Format: Any Include: ProductsDevelopersNewsForumsVideos
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.
KVRer

6 posts since 25 Jan, 2013
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
KVRian

767 posts since 30 Nov, 2008
DSPFilters (in my signature) calculates the phase response, you can crib the formula from there:

https://github.com/vinniefalco/DSPFiltersDemo

KVRAF

3648 posts since 8 Mar, 2004, from Berlin, Germany
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)
KVRer

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

6 posts since 25 Jan, 2013
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
First tests look promising so I hope I am on the right track now.
Thanks again!
KVRer

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

3648 posts since 8 Mar, 2004, from Berlin, Germany
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 factorHr = g*(cb*ca + sb*sa);  // real part of HHi = g*(cb*sa - sb*ca);  // imaginary part of HM  = sqrt(Hr^2 + Hi^2);  // magnitude`
KVRer

6 posts since 25 Jan, 2013
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)`
KVRAF

3648 posts since 8 Mar, 2004, from Berlin, Germany
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*sbH^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)
KVRAF

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

3648 posts since 8 Mar, 2004, from Berlin, Germany
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.
KVRer

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