Magnetude designed biquad problem !

DSP, Plugin and Host development discussion.
RELATED
PRODUCTS

Post

Hello everybody. Please be patient with my english witch is poor.

A few days back I was asking myself if it's possible that there is a better method than orfanidis.

I fall on a paper by Robin Schmidt the name is : "Digital Biquad Design by Magnitude Requirements" .http//www.rs-met.com/documents/dsp/Digi ... ements.pdf

The point of this paper is :

g2(ω)=ˆ|H(ejω)|2 =
b0^2 +b1^2 +b2^2 + 2*b1(b0 + b2)cos(ω)+2b0b2 cos(2ω)
------------------------------------------------------
a0^2 +a1^2 +a2^2 + 2*a1(a0 + a2)cos(ω)+2a0a2 cos(2ω)


we define B0 = b0^2+b1^2+b2^2
B1 = b1(b0 +b2)
B2 = b0*b2

(idem fo A0 A1 & A2)

and U(w) = 2* cos(w)
and V(w) = 2*cos(2*w)

we can make pairs of frequencies & magnetude : let's define :

=> p0 = |H(w0)^2|
=> p1 = |H(w1)^2|
=> p2 = |H(w2)^2|
=> p3 = |H(w3)^2|
=> p4 = |H(w4)^2|

in fact with our defining of B0 B1 B2 A0 A1 & A2 we can write :

|h(w)^2|
= B0+B1*U(w)+B2*V(w)
------------------
A0+A1*U(w)+A2*V(w)
<===>
p0 =
B0+B1*U(w0)+B2*V(w0)
------------------
A0+A1*U(w0)+A2*V(w0)

lets define A0 = 1 and multiply left and right part of the equation by the denum

(1+A1*U(w0)+A2*V(w0) ) * p0 = B0+B1*U(w0)+B2*V(w0)
=> p0 + A1*U(w0)*p0 + A2*V(w0)*p0 = B0+B1*U(w0)+B2*V(w0)
=> p0 = B0+B1*U(w0)+B2*V(w0) - A1*U(w0)*p0 - A2*V(w0)*p0

we can write the same for p1, p2, p3 and p4 so we got that system :

B0 + B1*U(w0) + B2*V(w0) - A1*U(w0)*p0 - A2*V(w0)*p0 = p0
B0 + B1*U(w1) + B2*V(w1) - A1*U(w1)*p1 - A2*V(w1)*p1 = p1
...
B0 + B1*U(w4) + B2*V(w4) - A1*U(w4)*p4 - A2*V(w4)*p4 = p4

we can solve this system and we obtain B0,B1,B2,A1 and A2

now we can express b0, b1 and b2 from B0,B1 and B2
we got B1 = b0^2+b1^2+b2^2
B2 = b0*b2 <==> b2 = B2/b0
B1 = b1(b0 +b2) <==> b1 = B1/(b0+b2) <==> b1 = B1/(b0+B2/b0)

so we can write :
b0^2 + b1^2 + b2^2 = B0
b0^2 + B1^2/(b0^2+2*B2+B2^2/b0^2) + b2^2/b0^2 = B0

and so :

b0^2 + B1^2/(b0^2+2*B2+B2^2/b0^2) + b2^2/b0^2 - B0 = 0

The paper says it's a root finding problem. I found b0^2, take the square root of it, and get next b1 & b2.

But in fact, the finding of B0, B1, B2, A1, A2 always give good results. But I'm not shure of the unicity of the response.
I mean, we need here a good precision, and sometimes
x + B1^2/(x+2*B2+B2^2/x) + b2^2/x - B0 don't get any result at all.

So i'm a bit confused, because 5 constraints, five Eq's gives five coefficients.
But is it normal "mathematicaly" that the founded coefficients don't give a formula with any root ?

Maybe there is a leak in the method ?

What do you think about

(I will take some time to answer I will travel next week !)

P.S. : I've made an implementation of my work, it's in PJ (XCode project but there is sources,but take care the code is really messy !)

https://dl.dropbox.com/u/4419982/BiquadTools.zip

Post

lereme wrote:sometimes
x + B1^2/(x+2*B2+B2^2/x) + b2^2/x - B0 don't get any result at all.
it's a nonlinear equation for x (= b0^2), parametrized on the 3 variables B0, B1, B2. when i plotted the graph for various choices of these 3 parameters, i observed that it may cross the x-axis once or twice or not at all (depending on the values of the parameters). i'm not sure anymore, but maybe sometimes there are even 3 solutions. when there's more than one root, one is related to another by reflection of poles/zeros in the unit circle and/or sign changes - which makes sense because these operations do not affect the magnitude. sometimes there is no solution at all. in these cases, you have requested a magnitude response which is not possible to obtain with a biquad. consider a requirement to have a notch somewhere - i.e. a zero on the unit circle. this actually requires two complex conjugate zeros - so it already uses up the two available zeros of the biquad which implies that you can't obtain zero amplitude at some additional frequency anymore. so, with this design method, it's vital to select constraints that are actually possible to meet.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Okay thank you Robin.
It's really interesting, but remembering lagrange and taylor series, you have a 5 degree equation she can go by 4 points !

so 6 parameters of liberty (6 coefficients) => 5 points

But in fact I know you can't with a biquad do 2 notch filters...

can you be a bit more clear in your explainations ? I'm really interested in !

Post

lereme wrote:Okay thank you Robin.
It's really interesting, but remembering lagrange and taylor series, you have a 5 degree equation she can go by 4 points !

so 6 parameters of liberty (6 coefficients) => 5 points
i'm not sure why you appeal to taylor series in this context and what about lagrange (i know him mainly from his "multipliers" in constrained optimization). it's just 5 (linear) equations for 5 unknowns (B0, B1, B2, A1, A2) with some additional mumbo jumbo to obtain the filter coefficients from these intermediate variables.
can you be a bit more clear in your explainations ?
i'm not sure what exactly the question is.
Last edited by Music Engineer on Sun Jan 20, 2013 10:16 pm, edited 1 time in total.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

for exemple : I design my constraints

Code: Select all

double w0 = 2.f*M_PI*cutoff/samplingRate;
	double Dw = w0/(Q);
	double G = gain;
	double GB = gain-3db

	double F = fabs(G*G - GB*GB);
	double G00 = fabs(G*G - G0*G0);
	double F00 = fabs(GB*GB - G0*G0);

	double num = G0*G0 * powf(w0*w0 - M_PI*M_PI,2.f) + G*G * F00 * M_PI*M_PI * Dw*Dw / F;
	double den = powf(w0*w0 - M_PI*M_PI,2.f) + F00 * M_PI*M_PI * Dw*Dw / F;

	double G1 = sqrt(num/den);
so I got DC gain, nyquist gain G0, -3db gain, and peak gain...
What's the différence with orfanidis ?

Post

lereme wrote: What's the différence with orfanidis ?
if you refer to his "prescribed nyquist frequency gain" approach to bell filters, the main difference is that his constraints include a zero-derivative (of the magnitude response) at the peak/dip together with 4 constraints on the magnitude itself whereas my approach uses 5 constraints on the magnitude itself.

edit: aah - and furthermore he operates mainly in the analog s-domain with the anticipation to use the bilinear transform whereas my approach is directly in the z-domain.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Sorry you answered before I sent this. No I was wrong, it's pretty clear.

The principal design flaw is to be really precise in fact.

I was thinking in fact the z equation if of the 5th degree (y[n] depends of x[n],x[n-1],x[n-2] and y[n-1], y[n-2], so I considered strangely we were in a 5th degree system.

Bu for the B0,B1,B2, A1,A2, we completely are in a linear system.

so shurely we can find some solutions with are not solvable at the stage of the root founding problem.but with realistic constraints we can found everything.

The question was about the complex conjugates, and the unit circle.

I didn't get any algorithm output with complex roots... when this append ?

Post

And Thank you for this paper, it's really interesting to play with these concepts.

Post

lereme wrote:I didn't get any algorithm output with complex roots... when this append ?
well, if i understand this right - you would not expect to get complex results for any of the A,B or a,b variables for any workable filter. a,b are the real coefficients. still, these real coefficients may of course correspond to complex poles/zeros of the filter as a whole. ..which is why you should also solve the corresponding quadratic equation for these poles/zeros, look at whether they are inside or outside the unit circle and reflect them, if outside.

edit: nah - that's actually not the reason why. real poles outside the unit circle are just as bad. you should do it, however.
Last edited by Music Engineer on Sun Jan 20, 2013 10:36 pm, edited 1 time in total.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

lereme wrote:And Thank you for this paper, it's really interesting to play with these concepts.
:tu: indeed, indeed. it's very satisfying to sit a couple of hours or days at the desk manipulating equations and finally arrive at some piece of code that works.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Okay. I've got difficulties with the english math vocab.

But my point with the [edit]coefficients[/edit] founded at the end (a0,a1,a2,b0,b1,b2)

How can you see the imaginary / real part of 6 real numbers ?

I understand how to take the complex conjugate, or reflect inside the circle a complex number... but an real ?

My head's burning since 4 days actually !!

It's not really intutive to understand how behave zeros and poles intutively...
Last edited by lereme on Sun Jan 20, 2013 10:46 pm, edited 1 time in total.

Post

Robin from www.rs-met.com wrote:which is why you should also solve the corresponding quadratic equation for these poles/zeros, look at whether they are inside or outside the unit circle and reflect them, if outside.
I've never heard about that, neither see any documentation about that.

Theres a lot of things from s to z. but there I'm in deep void !

Post

lereme wrote:
I've never heard about that, neither see any documentation about that.

Theres a lot of things from s to z. but there I'm in deep void !
this has nothing to do with s to z transforms. i mention it only briefly in the final remark on "stable, minimum-phase, normalized" filters. basically, you have 3 coefficients b0, b1, b2 for the numerator (and likewise for the denominator). the quadratic equation b0 + b1*z + b2*z^2 has two roots which are the zeros of the filter. if they are outside the unit circle, the filter is not minimum phase. worse, if the same is true for the denominator, i.e. a0 + a1*z + a2*z^2 has roots outside the unit circle, the filter will be unstable. you can, however, reflect them into the unit circle without changing the magnitude response. this reflection actually just amounts to taking the reciprocal of the root in question (actually, reciprocal of the complex conjugate, but in this case, you'll have two complex conjugate zeros/poles anyway, so it doesn't matter)
Last edited by Music Engineer on Sun Jan 20, 2013 11:24 pm, edited 1 time in total.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post

Yeah I see.

I have to found the values of z in b0 + b1*z + b2*z^2 = 0

I found this with b^2-4ac but if delta is complex I got complex solution.

So i've got roots, complex or real. admit complex roots !
I've got now to change my coefficients for the original root to be mapped to his complex conjugate or reflected inside the unit circle ?

(tell me I'm righ please !)

Post

lereme wrote:Yeah I see.

I have to found the values of z in b0 + b1*z + b2*z^2 = 0

I found this with b^2-4ac but if delta is complex I got complex solution.

So i've got roots, complex or real. admit complex roots !
I've got now to change my coefficients for the original root to be mapped to his complex conjugate or reflected inside the unit circle ?

(tell me I'm righ please !)
yes - well almost. "delta" doesn't have to be complex - negative delta is enough for the root to be complex.
My website: rs-met.com, My presences on: YouTube, GitHub, Facebook

Post Reply

Return to “DSP and Plugin Development”