## Phase response of a biquad

12 posts
• Page

**1**of**1**- 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

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

https://github.com/vinniefalco/DSPFiltersDemo

Open Source: DSPFilters, LayerEffects, VFLib, SimpleDJ

- 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

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!

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

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

Here is the ca, sa, cb, sb calculation from biquad parameters (normalized for a0=1):

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

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

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

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

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