How to split high-order IIR to multiple biquads?

[KG_is_back]

the tittle says it all... let's say I have coefficients of a N-th order IIR filter and I what to split it into cascade of biquads. How do I compute the coeffiecients? Is there any particular algorithm that can be implemented automatically? (like: I enter IIR coefficients and it outputs an array of multiple biquad coefs)

[martinvicanek]

KG, you need to locate the poles and zeros of the original nth order filter, then collect the complex conjugate pairs to form biquads. This amounts to a factorization of polynomials. There are ready-to-use routines for that, some are even accessible online. Dont't know about DLLs which you could integrate with FS, though.

[KG_is_back]

To combine two filters one just have to multiply their Transfer functions. That means (b_0 + b_1*z^1+b_2*z^2)*(b_0 + b_1*z^1+b_2*z^2) and the same with A coefficients. Maybe the process can be reversed by dividing the transfer function with biquad transfer function so it eliminates two last coefficients... It would be hard to find such function though (and even harder to create program that finds it). Is there any algorithm that might help? maybe iterative one...

[martinvicanek]

Maybe the process can be reversed by dividing the transfer function with biquad transfer function so it eliminates two last coefficients... It would be hard to find such function though (and even harder to create program that finds it). Is there any algorithm that might help? maybe iterative one...

Yes, that is basically what is known as deflation:

Code: Select all

1. Start with a nth order polynomial P(x) = a0 + a1*x + a2*x^2 + ... + an*x^n
2. Find one root x0 [Note: x0 is a root of P if P(x0) = 0. There are n roots in total.]
3. Divide the original polynomial by (x - x0) to get a new polynomial of degree n-1
4. Repeat steps 2.-3. until all roots are found

Find more information and programs here:
http://en.wikipedia.org/wiki/Root-findi ... olynomials
http://www.akiti.ca/PolyRootRe.html
http://www.hvks.com/index.html
http://apps.nrbook.com/c/index.html go to page 369

[KG_is_back]

This is cool! I can actually do this!
To find the root is simple - Newton method. It would be best if it could be done in ruby as it uses DP-math. I'll probably have to start learning ruby after all.
The division of the polynomials can probably also be done very simply as you are alway dividing by (1*x-x0) where x0 is a constant.
The formula (1*x-x0) can be also written as (b0*z^0 + b1*z^-1) where b0=1 and b1=-x0 so it is nominator of the transfer function of first order filter. The same algorithm would be used for A coefficients.

This way N-th order filter will be split into N first-order filters and they can easily be summed to N/2 biquads.

[KG_is_back]

Finally after such a long time I have it figured out. I've created module, that receives coefficients of higher order IIR and creates 5 arrays of biquad coefficients (how much is needed). After casual testing it shows pretty good - the difference is under 60dB across the whole spectrum so it's basically identical.

Here is an example schematic: it has 3 parallel chains:
First one has two biquads of which response you control.
second one is a 4th order filter combined from those original biquads.
third one has two biquads that are reconstructed from the 4th order filter coefficients. Note that these 2 biquads may have different set of coefficients then the original one, because there are multiple ways you can combine 4poles and zeroes into two 2pole-2zeroes filter (biquad). However the filter responses are (almost) identical, so it works fine.

[MyCo]

Interesting, but the schematic is probably the worst I've seen in a while

To calculate the difference, I think you have to divide not subtract. You are operating in linear scale (not log)

[KG_is_back]

Interesting, but the schematic is probably the worst I've seen in a while

To calculate the difference, I think you have to divide not subtract. You are operating in linear scale (not log)

I've set it up in like 5 minutes, just to find out if it's actually working (and it does, so hureyyyyy ).

I'm pretty sure that difference=subtraction not division. Simply because some of the values might be zero (and probably are, near the end of the IR) and dividing by zero gives NaN. I'm working out the difference in time domain, not frequency domain. Still when you subtract the original and reconstructed you get something that peaks around 10^-6 which is basically a rounding error.

[martinvicanek]

Respect, that's quite some serious stuff there, KG! Although I haven't tested it yet (too much spaghetti for me ) Do you think it could be done with a 14th order allpole IIR at hop(1024)? I am thinking of an LPC application...

[KG_is_back]

Respect, that's quite some serious stuff there, KG! Although I haven't tested it yet (too much spaghetti for me ) Do you think it could be done with a 14th order allpole IIR at hop(1024)? I am thinking of an LPC application...

Yes, it should work for 14th order allpole IIR. However, I do not know what you mean with that hop(1024). The thing is done in ruby via event method. But it could possibly be moded to work in sync with frames.
By the way the A and B coefficients are calculated completely separately, so you may look into the module and delete a bunch of stuff so you can input A and B separately. Just remember for the individual A and B sections you have to input ALL coefficients - biquads are usually normalized to a0=1, but here you have to input that too - it is not done automatically. In case your schematic is not normalized to a0 do so or replace the A calculator with B calculator and then normalize the output biquads.