Phase response of a biquad
-
- KVRer
- Topic Starter
- 21 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/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
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
- KVRian
- 775 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
My Open Source:
Beast, rippled, DSPFilters, LayerEffects, SimpleDJ
Beast, rippled, DSPFilters, LayerEffects, SimpleDJ
-
Music Engineer Music Engineer https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=15959
- KVRAF
- 4265 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
- Topic Starter
- 21 posts since 25 Jan, 2013
Thanks!thevinn wrote:DSPFilters (in my signature) calculates the phase response, you can crib the formula from there:
I will try digging into the source code.
Nice software by the way!
-
- KVRer
- Topic Starter
- 21 posts since 25 Jan, 2013
Thanks Robin, that looks straightforward enough indeedRobin 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)
First tests look promising so I hope I am on the right track now.
Thanks again!
-
- KVRer
- Topic Starter
- 21 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!
-
Music Engineer Music Engineer https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=15959
- KVRAF
- 4265 posts since 8 Mar, 2004 from Berlin, Germany
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: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.
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
- Topic Starter
- 21 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))
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 Music Engineer https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=15959
- KVRAF
- 4265 posts since 8 Mar, 2004 from Berlin, Germany
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):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))
Code: Select all
cb - j*sb cb + j*sb
H^2(w) = H(w) * H'(w) = ----------- * -----------
ca - j*sa ca + j*sa
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);
}
-
Music Engineer Music Engineer https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=15959
- KVRAF
- 4265 posts since 8 Mar, 2004 from Berlin, Germany
OK - i just verified it. it works. we generally have (z/w)' = z' / w'. so both formulas should be equivalent.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
that wasn't obvious to me right now and wikipedia doesn't list that identity either:
http://en.wikipedia.org/wiki/Complex_conjugate
-
Music Engineer Music Engineer https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=15959
- KVRAF
- 4265 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
- Topic Starter
- 21 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
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
Yes, stupid me...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.
I summed several biquads to exceed 90° and could not see any wrapping... of course as it would be a per biquad wrapping...