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

    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

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.


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.

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.

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:

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.

Are filter coefs the same for each form?

Yes, exactly the same calculations.

Just to get something moving I tried your approach, but IDK how it is unstable.

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.

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?

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.

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:


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.

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.

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.

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.