biquad basics... looking for a lesson

stefan · 9 Feb 2006

They are "equivalent". They have just swapped the a and the b coefficients. They are also swapped in the coefficient calculation, so you end up with the same coeffs at the same places really. But it is confusing nonetheless.

In the first picture, they have left out "a0", which is the "gain" of the entire filter. This is probably because they have normalized the coefficients by a0, which you do by dividing all the other coefficients with a0. or. Actuall, not "all the other", but "all coefficients". a0/a0 ends up being 1, hence, that one can be stripped out from the schematic.

a0 would normally reside right above "-a1", as a gain stage. Sometimes people don't use "-a1" and "-a2", but just "a1" and "a2" instead, and make sure that the actual coefficient have the right sign instead.

Now, to mess things up for you, in the ParametricEQ and the biquad analysers, we use the 2nd picture version, but unnormalized, so b0 is kept and used. The normalization is made within the filter / analysers.

Sorry if biquads are a bit messy. Not my fault there's no "standard" way to do it 🙂

oddson · 9 Feb 2006

Thanks, I almost get it now...

Shouldn't the SM primative use the same conventions as the SM Freq/Phase analyzers?

stefan · 9 Feb 2006

oddson wrote: Thanks, I almost get it now...

Shouldn't the SM primative use the same conventions as the SM Freq/Phase analyzers?

I guess... I didn't even think about the primitive when I implemented the EQ or when I told jorge how to make the analysers. Is it not better to use the implementations in the EQ? (since there are two flavors I mean, stream coeffs and green float data coeffs).

I don't know really 🙂

oddson · 9 Feb 2006

It's better for all the parts of SM to be internally consistent...

I'm using your code for the EQ for learning purposes, I just noticed the lack of a 'b0' coefficient in the 2P2Z today...

I'm still thinking your first response through... expect more questions..

Thanks! 😀

oddson · 9 Feb 2006

As promised, more questions...

From the biquad code in the EQ module...

a0b0 = a0*b0inv; a1b0 = a1*b0inv; a2b0 = a2*b0inv; b1b0 = b1*b0inv; b2b0 = b2*b0inv; ... y0 = a0b0*x0 + a1b0*x1 + a2b0*x2 - b1b0*y1 - b2b0*y2; ... out = y0;

First, since b0 is inverted to become b0inv, how is this the parameter from a gain stage?

Second, why would you caluclate all those intermediaries instead of using the distributive property of multiplication?

My algebra is 20 years rusty but, couldn't you use:

y = b0inv * ( a0*x0 + a1*x1 + a2*x2 - b1*y1 - b2*y2 );

rl · 10 Feb 2006

A little suggesstion to the SM team: Expand the tooltip text of the 2P2Z module, so that it shows the biquad's transfer function. I'm always confused when using it and have to cross-check with an example jorge posted here some months ago.

oddson wrote: My algebra is 20 years rusty but, couldn't you use:

y = b0inv * ( a0*x0 + a1*x1 + a2*x2 - b1*y1 - b2*y2 );

I think so too. Your calculation also uses less multiplies -> less clock cyles. But care has to be taken on lengthy expressions since there is only a limited set of SSE registers. The code module gives no indication when an expression cannot be mapped to the SSE register set.

stefan · 10 Feb 2006

oddson wrote: As promised, more questions...

From the biquad code in the EQ module...

a0b0 = a0*b0inv; a1b0 = a1*b0inv; a2b0 = a2*b0inv; b1b0 = b1*b0inv; b2b0 = b2*b0inv; ... y0 = a0b0*x0 + a1b0*x1 + a2b0*x2 - b1b0*y1 - b2b0*y2; ... out = y0;

First, since b0 is inverted to become b0inv, how is this the parameter from a gain stage?

Actually, I'm not quite sure about this. The different schematics I find _all_ have normalized a0 or b0 to 1.

From RBJ's "audio cookbook" (which is a good read): "The most straight forward implementation would be the "Direct Form 1" (Eq 2):

y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]"

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

I could get the discrete version out of the "before normalization" transfer function, if you really really really need it 🙂 In the end, we will turn out calculating the same thing though.

oddson wrote: Second, why would you caluclate all those intermediaries instead of using the distributive property of multiplication?

My algebra is 20 years rusty but, couldn't you use:

y = b0inv * ( a0*x0 + a1*x1 + a2*x2 - b1*y1 - b2*y2 );

Yes, and it'd be way less stupid than what I did 🙂

Here's a "cleaner" implementation:

streamin in; streamout out;

streamin freq; streamin linGain; streamin Q; streamin samplerate; // mem float y0,y1,y2,x0,x1,x2; // coeffs float a0,a1,a2,b0,b1,b2; // 1/b0 optimisation float w0,alpha,A,inv2Pi,alphaDivA; in = in + 0.00000000001; //hop(1) //{ inv2Pi = 0.5/3.1415926535897932384626433832795; A = linGain; w0 = 2*3.1415926535897932384626433832795*freq/samplerate; alpha = sin1(w0*inv2Pi)/(2*Q); a0 = 1 + alpha*A; alphaDivA = alpha/A; a1 = -2 * cos1(w0*inv2Pi); a2 = 1 - alpha*A; b0 = 1 + alphaDivA; b1 = a1; b2 = 1 - alphaDivA; //}

x0 = in; y0 = (a0*x0 + a1*x1 + a2*x2 - b1*y1 - b2*y2) / b0; // push vars to mem x2 = x1; x1 = x0; y2 = y1; y1 = y0; out = y0; in = in - 0.00000000001;

If you want to use an "optimized" version, enable the hop section and use hop(16) or something 🙂

rl · 10 Feb 2006

BTW: the code snippet

float in; in = in + 0.00000000001;

translates to:

float in=0.0000000000; float F0P0000000000=0.0000000000; movaps xmm0,in; addps xmm0,F0P0000000000; //Assignment> sLeft=xmm0 movaps in,xmm0;

(in=in+0)

so, was this a mistake, or are you using some secret wisdom? 😉

stefan · 10 Feb 2006

mistake 🙁

infuzion · 10 Feb 2006

rl wrote: But care has to be taken on lengthy (Code) expressions since there is only a limited set of SSE registers. The code module gives no indication when an expression cannot be mapped to the SSE register set.

I never thought of that, but it makes sence... how many registers are there? How can we tell if we hit the limit? TIA

rl · 10 Feb 2006

Jorge has given an example here

oddson · 10 Feb 2006

rl wrote: ...so, was this [in = in - 0.00000000001] a mistake, or are you using some secret wisdom?

That's that thing I found when trying to make a A/D envelope for my Karplus... the input variables are read-only, correct?

stefan wrote: ...[the two diagrams] are "equivalent". They have just swapped the a and the b coefficients. They are also swapped in the coefficient calculation, so you end up with the same coeffs at the same places really. But it is confusing nonetheless.

OK, first I know they are functionally equivalent, but are the coefficient groups actually interchangable? I'm dubious about ending up with the "same coeffs at the same places"

Especially given the first direct form y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]"

In this, only a0 is in a denominator position... so why are we dividing by b0?

But if you've used the other version of the output transform, how can you use the coefficient recipes from the cookbook? [e.i. shouldn't a0 be b0 etc...]

But clearly it's working for your peaking filters, but I can't make anything else work by following the cookbook...

😕 ^2

infuzion · 11 Feb 2006

rl wrote: Jorge has given an example here

Dang, that was almost a year ago! Before my time. thanks for the link 🙂

stefan · 11 Feb 2006

oddson wrote: rl wrote: ...so, was this [in = in - 0.00000000001] a mistake, or are you using some secret wisdom?

That's that thing I found when trying to make a A/D envelope for my Karplus... the input variables are read-only, correct?

stefan wrote: ...[the two diagrams] are "equivalent". They have just swapped the a and the b coefficients. They are also swapped in the coefficient calculation, so you end up with the same coeffs at the same places really. But it is confusing nonetheless.

OK, first I know they are functionally equivalent, but are the coefficient groups actually interchangable? I'm dubious about ending up with the "same coeffs at the same places"

Especially given the first direct form y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]"

In this, only a0 is in a denominator position... so why are we dividing by b0?

But if you've used the other version of the output transform, how can you use the coefficient recipes from the cookbook? [e.i. shouldn't a0 be b0 etc...]

But clearly it's working for your peaking filters, but I can't make anything else work by following the cookbook...

😕 ^2

The thing is that we also swap the _calculations_ of the coefficients.

That is, in the second schematic (the picture) b1 is calculated like a1 is calculated for the first schematic. And so forth...

oddson · 11 Feb 2006

Sorry if I'm seeming to be dense, but sometimes I am...

Is this the correct block diagram for the SM biquad freq/phase analyzers? That is, would this diagram lead to the correct interpretation of the coefficients used in the analyzers?

Are the signs correct on the b1 and b2? So positive values of the coefficients result in inverted phase signals?

For this diagram would the (direct, non-normalized) output function be?:

y(n) = a0*b0*x(n) + a1*b0*x(n-1) + a2*b0*x(n-2) - a0*b0*b1*y(n-1) - a0*b0*b2*y(n-2)

Am I correct so far? edit - wrong see below

rl · 12 Feb 2006

I'm usually using this post as a reference: Equalizer

Borogove · 14 Feb 2006

Then there's Direct Form II...

oddson · 14 Feb 2006

Borogove wrote: Then there's Direct Form II...

Right, but that's with the other assignation of the coefficients (ur, um, the 'a's and 'b's are switched). It also lacks b0 (a0 with it's conventions). Direct Form II is the form in Jorge's block flow diagram that rl's post above links to.

sambean · 16 Feb 2006

For this diagram would the (direct, non-normalized) output function be?:

y(n) = a0*b0*x(n) + a1*b0*x(n-1) + a2*b0*x(n-2) - a0*b0*b1*y(n-1) - a0*b0*b2*y(n-2)

Am I correct so far?

Should be:

y(n) = a0*b0*x(n) + a1*b0*x(n-1) + a2*b0*x(n-2) - b0*b1*y(n-1) - b0*b2*y(n-2)

(to get this, follow the paths from the delays to the output..)

Which becomes:

y(n) = a0'*x(n) + a1'*x(n-1) + a2'*x(n-2) - b1'*y(n-1) - b2'*y(n-2)

if a0' = a0*b0 etc...

Try not to think of the coefficients as the 'a coefficients' and the 'b coefficients' as I don't think there's a set convention. To distinguish them think of them as the 'feedforward coefficients' and the 'feedback coefficients'. In most cases the biquads will be normalised so you'll have 3 feedfoward (either a's or b's) and 2 feedback coefficients...

Cheers

Sam

oddson · 16 Feb 2006

y(n) = a0*b0*x(n) + a1*b0*x(n-1) + a2*b0*x(n-2) - b0*b1*y(n-1) - b0*b2*y(n-2)

(to get this, follow the paths from the delays to the output..)

I thought I did... 😕 how do you go through the feedback delays without also going through the a0 amp?

😮 Oh wait! (light goes on) since y(n) has gone through the whole circut, so has y(n-1) etc...

I still haven't worked out how the normalized output function expressed in the non-normalized coefficients is:

y[n] = (a0/b0)*x[n] + (a1/b0)*x[n-1] + (a2/b0)*x[n-2] - (b1/b0)*y[n-1] - (b2/b0)*y[n-2]

...but at least I'm back on track... my misread of the schematic had me stopped dead in getting any further in understanding this.

Thanks Sam 🙂

oddson · 16 Feb 2006

Stefan...

Just to make sure I'm gettting it... to implement an RBJ filter from the cookbook I have to swap the a/b coefficients to use your output function...

It only looks like you've copied it directly because the peaking filter happens to be symetrical in coeffs... right?

sambean · 17 Feb 2006

I still haven't worked out how the normalized output function expressed in the non-normalized coefficients is:

y[n] = (a0/b0)*x[n] + (a1/b0)*x[n-1] + (a2/b0)*x[n-2] - (b1/b0)*y[n-1] - (b2/b0)*y[n-2]

OK, found it...

With a digital system constructed from delays etc, for example a filter, there is a mathematical technique you can use called the z-transform.

The z-transform is an equation which fully describes the filter, with it you can extract the frequency response, impulse response, etc..

Without going into too much detail, the z-transfrom of the non-normalised biquad is:

H(z) =

{a0 + a1 * z^(-1) + a2 * z^(-2)} --------------------------------------- {b0 + b1 * z^(-1) + b2 * z^(-2)}

if we now divide top and bottom by b0 we get:

H(z) =

{a0/b0 + a1/b0 * z^(-1) + a2/b0 * z^(-2)} --------------------------------------- {1 + b1/b0 * z^(-1) + b2/b0 * z^(-2)}

if we now use the inverse z-transform we get the function:

y[n] = (a0/b0)*x[n] + (a1/b0)*x[n-1] + (a2/b0)*x[n-2] - (b1/b0)*y[n-1] - (b2/b0)*y[n-2]

Hey presto!

As this has the same z-transfer function as the non-normalised version it is essentially the same filter..

😀

stefan · 17 Feb 2006

oddson wrote: Stefan...

Just to make sure I'm gettting it... to implement an RBJ filter from the cookbook I have to swap the a/b coefficients to use your output function...

It only looks like you've copied it directly because the peaking filter happens to be symetrical in coeffs... right?

Seems about right, yeah.