camsr
KVRAF

6805 posts since 16 Feb, 2005
How can I turn this into a 4th order filter?

yk = a0 xk + a1 xk−1 + a2 xk−2 + b1 yk−1 + b2 yk−2

Here is my code so far, I am drawing a blank

Code: Select all
`    while (j < 2)    {        for (int i = 0; i < blockSize; i++)        {            blp1[2] = blp1[1];            blp1[1] = blp1[0];            blp1[0] = dynbuf[i];            dynbuf1[i] = (blpc1[0] * blp1[0]) + (blpc1[1] * blp1[1]) + (blpc1[2] * blp1[2]) + (blpc1[3] * blp1[3]) + (blpc1[4] * blp1[4]);        }        j++;    }`

thanks
LemonLime
KVRist

231 posts since 15 Apr, 2012, from Toronto, ON
Here's how I handle the cascading in my filters. I implement using direct form II, though, because then you only need 1 delay line. If you stick with direct form I, same principle applies, but just use two delay lines.

Code: Select all
`for (j=0; j < numCascades; j+=2) {    W = inputSample - b1*z[j] - b2*z[j+1];    outputSample = a0*W + a1*z[j] + a2*z[j+1];    z[j+1] = z[j];    z[j] = W;}`

z is the delay buffer, and has a size = max # of cascades * 2, because 2 is the base order of the biquad difference equation. Also, I have my block loop on the outside of that, not on the inside like you do, so it cascades each sample as it goes through the block.
camsr
KVRAF

6805 posts since 16 Feb, 2005
The problem I realized quickly was that the state was corrupted at the second cascade. Appearantly I will have to cascade per sample or code in some mumbo jumbo to do it my way. Doing it my way worked fine otherwise. I already have too many variables for my liking.
camsr
KVRAF

6805 posts since 16 Feb, 2005
LemonLime wrote:Here's how I handle the cascading in my filters. I implement using direct form II, though, because then you only need 1 delay line. If you stick with direct form I, same principle applies, but just use two delay lines.

Code: Select all
`for (j=0; j < numCascades; j+=2) {    W = inputSample - b1*z[j] - b2*z[j+1];    outputSample = a0*W + a1*z[j] + a2*z[j+1];    z[j+1] = z[j];    z[j] = W;}`

z is the delay buffer, and has a size = max # of cascades * 2, because 2 is the base order of the biquad difference equation. Also, I have my block loop on the outside of that, not on the inside like you do, so it cascades each sample as it goes through the block.

Why are you only using half of the z[] array?
LemonLime
KVRist

231 posts since 15 Apr, 2012, from Toronto, ON
camsr wrote:
LemonLime wrote:Here's how I handle the cascading in my filters. I implement using direct form II, though, because then you only need 1 delay line. If you stick with direct form I, same principle applies, but just use two delay lines.

Code: Select all
`for (j=0; j < numCascades; j+=2) {    W = inputSample - b1*z[j] - b2*z[j+1];    outputSample = a0*W + a1*z[j] + a2*z[j+1];    z[j+1] = z[j];    z[j] = W;}`

z is the delay buffer, and has a size = max # of cascades * 2, because 2 is the base order of the biquad difference equation. Also, I have my block loop on the outside of that, not on the inside like you do, so it cascades each sample as it goes through the block.

Why are you only using half of the z[] array?

Oops, typo. I should have looked at my code closer.
The for loop should be:
Code: Select all
`for (j = 0; j < numCascades * 2; j+=2) // order 2 biquad`

So then, with a 2-stage cascade, z[] has a size of 4, and elements 0,1 are the 1- and 2-sample delays respectively and 2,3 are the 1- and 2-sample delays at the second cascade stage.
camsr
KVRAF

6805 posts since 16 Feb, 2005
Are filter coefs the same for each form?
LemonLime
KVRist

231 posts since 15 Apr, 2012, from Toronto, ON
Yes, exactly the same calculations.
camsr
KVRAF

6805 posts since 16 Feb, 2005
Just to get something moving I tried your approach, but IDK how it is unstable.

Code: Select all
`for (j = 1; j <= 4; j += 2)        {            blp1[0] = dynbuf1[i] - (blpc1[3] * blp1[j]) - (blpc1[4] * blp1[j+1]);            dynbuf1[i] = (blpc1[0] * blp1[0]) + (blpc1[1] * blp1[j]) + (blpc1[2] * blp1[j+1]);            blp1[j+1] = blp1[j];            blp1[j] = blp1[0];        }`

dynbuf1[i] is a block sample and blp1[0] would be W in your code.

Here is the stuff I am using to get the coefs: http://unicorn.us.com/alex/2polefilters.html

EDIT: Coefs were not calculated correctly.
camsr
KVRAF

6805 posts since 16 Feb, 2005
I have the DF1 working, but not cascaded yet. The website I linked previously didn't exactly make the implementation simple, and I had to work out a few things myself. I am still confused by the DF2, even though the flowcharts explain it. Is there a simple way to take DF1 to DF2?
LemonLime
KVRist

231 posts since 15 Apr, 2012, from Toronto, ON
One thing that I run into many times with filters is that some books/resources give the difference equations as:
y(n) = a0x(n) + a1x(n-1) + a2x(n-2) - b1y(n-1) - b2y(n-2)
while others reverse the sign of the recursive section that might require you to switch the polarity of those coefficients (which is probably what you did to make it stable).

The way I sort of think about DF2 is that instead of calculating the feedforward section followed by the feedback/recursive section, it reverses them. This puts the delays in the middle, and if you look at the intermediate diagram here ( http://www.bores.com/courses/intro/iir/5_form2.htm ), it shows that the 4 delay elements are redundant and can be reduced to 2.
mystran
KVRAF

4891 posts since 11 Feb, 2006, from Helsinki, Finland
You need two separate filters, with two separate sets of state variables (also known as "delays" by some lesser people). You just put them in series.

The sensible way to implement this, is to write a class for a generic biquad (structure doesn't matter). Then create two objects of your biquad class, then filter the signal with the first one first, and then with the second one.

edit: if you insist on having a class for a 4-pole, the sensible thing to do looks something like this:

Code: Select all
`class FourPole{private:    TwoPole a, b;public:    // ....    double filter(double in)    {        return a.filter(b.filter(in));    }};`
<- plugins | forum
LemonLime
KVRist

231 posts since 15 Apr, 2012, from Toronto, ON
mystran wrote:The sensible way to implement this, is to write a class for a generic biquad (structure doesn't matter). Then create two objects of your biquad class, then filter the signal with the first one first, and then with the second one.

edit: if you insist on having a class for a 4-pole, the sensible thing to do looks something like this:

Code: Select all
`class FourPole{private:    TwoPole a, b;public:    // ....    double filter(double in)    {        return a.filter(b.filter(in));    }};`

This did occur to me as well. Is there anything inherently wrong or "un-optimal" the way I did it, by handling the cascading within the class?

I haven't put it though rigorous tests yet, but with preliminary performance tests it operates quite fast. I take about a 30% hit to execution time with each cascade stage. And for IIR filters, I see little reason to cascade more than 4 times (essentially an 8-order filter). At least for my purposes.

Though I'm always interested in ways to optimize.
mystran
KVRAF

4891 posts since 11 Feb, 2006, from Helsinki, Finland
You should always optimize developer productivity first, and code performance second. I would argue that the method I suggested is significantly easier to maintain in the long run (say, you realize how awful the numerical performance of direct forms is, and want to replace them with something slightly less catastrophic; with my suggested method all you need is to drop in another biquad implementation).

Anyway, if you observe significant run-time performance difference (one way or another), there's probably something wrong with your compiler settings. Normally the main bottleneck with modern computers tends to be pipeline latency.
<- plugins | forum
camsr
KVRAF

6805 posts since 16 Feb, 2005
LemonLime wrote:One thing that I run into many times with filters is that some books/resources give the difference equations as:
y(n) = a0x(n) + a1x(n-1) + a2x(n-2) - b1y(n-1) - b2y(n-2)
while others reverse the sign of the recursive section that might require you to switch the polarity of those coefficients (which is probably what you did to make it stable).

The way I sort of think about DF2 is that instead of calculating the feedforward section followed by the feedback/recursive section, it reverses them. This puts the delays in the middle, and if you look at the intermediate diagram here ( http://www.bores.com/courses/intro/iir/5_form2.htm ), it shows that the 4 delay elements are redundant and can be reduced to 2.

Actually, I just used the difference equation on that website in full. I think any polarity mangling is in the coef calculation, which might explain why when I used your code, it was unstable.
camsr
KVRAF

6805 posts since 16 Feb, 2005
mystran wrote:You need two separate filters, with two separate sets of state variables (also known as "delays" by some lesser people). You just put them in series.

The sensible way to implement this, is to write a class for a generic biquad (structure doesn't matter). Then create two objects of your biquad class, then filter the signal with the first one first, and then with the second one.

edit: if you insist on having a class for a 4-pole, the sensible thing to do looks something like this:

Code: Select all
`class FourPole{private:    TwoPole a, b;public:    // ....    double filter(double in)    {        return a.filter(b.filter(in));    }};`

Yes this is good for testing purposes. I started out writing a class but knew the implementation was going to be difficult for me, so I skipped it and put it in a function. Now I have to make a component.
Next

Moderator: Moderators (Main)