Plug-ins, Hosts, Apps,
Hardware, Soundware
Developers
(Brands)
Videos Groups
Whats's in?
Banks & Patches
Music Search
KVR

KVR Forum » DSP and Plug-in Development
Goto page 1, 2  Next
camsr
KVRAF
 Posted: Tue Nov 20, 2012 5:06 pm reply with quote
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
----
^ Joined: 16 Feb 2005  Member: #58183
LemonLime
KVRist
 Posted: Tue Nov 20, 2012 5:23 pm reply with quote
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.
^ Joined: 15 Apr 2012  Member: #278696  Location: Toronto, ON
camsr
KVRAF
 Posted: Tue Nov 20, 2012 7:00 pm reply with quote
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.
----
^ Joined: 16 Feb 2005  Member: #58183
camsr
KVRAF
 Posted: Tue Nov 20, 2012 7:05 pm reply with quote
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.

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?
----
^ Joined: 16 Feb 2005  Member: #58183
LemonLime
KVRist
 Posted: Tue Nov 20, 2012 7:37 pm reply with quote
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.

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:

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.
^ Joined: 15 Apr 2012  Member: #278696  Location: Toronto, ON
camsr
KVRAF
 Posted: Tue Nov 20, 2012 8:12 pm reply with quote
Are filter coefs the same for each form?
----
^ Joined: 16 Feb 2005  Member: #58183
LemonLime
KVRist
 Posted: Tue Nov 20, 2012 8:23 pm reply with quote
Yes, exactly the same calculations.
^ Joined: 15 Apr 2012  Member: #278696  Location: Toronto, ON
camsr
KVRAF
 Posted: Tue Nov 20, 2012 10:03 pm reply with quote
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.
----
^ Joined: 16 Feb 2005  Member: #58183
camsr
KVRAF
 Posted: Wed Nov 21, 2012 12:44 am reply with quote
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?
----
^ Joined: 16 Feb 2005  Member: #58183
LemonLime
KVRist
 Posted: Wed Nov 21, 2012 6:15 am reply with quote
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.
^ Joined: 15 Apr 2012  Member: #278696  Location: Toronto, ON
mystran
KVRAF
 Posted: Wed Nov 21, 2012 9:07 am reply with quote
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));
}
};
----
<- my plugins | my music -> @Soundcloud
^ Joined: 11 Feb 2006  Member: #97939  Location: Helsinki, Finland
LemonLime
KVRist
 Posted: Wed Nov 21, 2012 11:50 am reply with quote
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:

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.
^ Joined: 15 Apr 2012  Member: #278696  Location: Toronto, ON
mystran
KVRAF
 Posted: Wed Nov 21, 2012 1:19 pm reply with quote
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.
----
<- my plugins | my music -> @Soundcloud
^ Joined: 11 Feb 2006  Member: #97939  Location: Helsinki, Finland
camsr
KVRAF
 Posted: Wed Nov 21, 2012 3:11 pm reply with quote
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.
----
^ Joined: 16 Feb 2005  Member: #58183
camsr
KVRAF
 Posted: Wed Nov 21, 2012 3:13 pm reply with quote
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:

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.
----
^ Joined: 16 Feb 2005  Member: #58183
 KVR Forum Index » DSP and Plug-in Development All times are GMT - 8 Hours Printable version Page 1 of 2 Goto page 1, 2  Next