What is KVR Audio? | Submit News | Advertise | Developer Account

Options (Affects News & Product results only):

OS:
Format:
Include:
Quick Search KVR

"Quick Search" KVR Audio's Product Database, News Items, Developer Listings, Forum Topics and videos here. For advanced Product Database searching please use the full product search. For the forum you can use the phpBB forum search.

To utilize the power of Google you can use the integrated Google Site Search.

Products 0

Developers 0

News 0

Forum 0

Videos 0

Search  

cascading help

DSP, Plug-in and Host development discussion.

Moderator: Moderators (Main)

KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Tue Nov 20, 2012 5:06 pm

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
Image
KVRist
 
219 posts since 15 Apr, 2012, from Toronto, ON

Postby LemonLime; Tue Nov 20, 2012 5:23 pm

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.
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Tue Nov 20, 2012 7:00 pm

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.
Image
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Tue Nov 20, 2012 7:05 pm

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?
Image
KVRist
 
219 posts since 15 Apr, 2012, from Toronto, ON

Postby LemonLime; Tue Nov 20, 2012 7:37 pm

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.
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Tue Nov 20, 2012 8:12 pm

Are filter coefs the same for each form?
Image
KVRist
 
219 posts since 15 Apr, 2012, from Toronto, ON

Postby LemonLime; Tue Nov 20, 2012 8:23 pm

Yes, exactly the same calculations.
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Tue Nov 20, 2012 10:03 pm

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.
Image
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Wed Nov 21, 2012 12:44 am

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?
Image
KVRist
 
219 posts since 15 Apr, 2012, from Toronto, ON

Postby LemonLime; Wed Nov 21, 2012 6:15 am

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.
KVRAF
 
3927 posts since 11 Feb, 2006, from Helsinki, Finland
 

Postby mystran; Wed Nov 21, 2012 9:07 am

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));
    }
};
Image <- plugins | forum
KVRist
 
219 posts since 15 Apr, 2012, from Toronto, ON

Postby LemonLime; Wed Nov 21, 2012 11:50 am

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.
KVRAF
 
3927 posts since 11 Feb, 2006, from Helsinki, Finland
 

Postby mystran; Wed Nov 21, 2012 1:19 pm

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.
Image <- plugins | forum
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Wed Nov 21, 2012 3:11 pm

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.
Image
KVRAF
 
3920 posts since 16 Feb, 2005

Postby camsr; Wed Nov 21, 2012 3:13 pm

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.
Image
Next

Moderator: Moderators (Main)

Return to DSP and Plug-in Development