IIR to FIR, is a best fit polynomial usually necessary?

Question

I'm messing around with IIR/FIR filters and want to convert the former to the latter.

I set up a classic impulse response calculation.

        X[4] = 1.0

        Y[0] = 0.0
        Y[1] = 0.0

        for n in range( 2, L ):
          Y[n] = 0.5 * X[n] + 0.3 * X[n-1] + 0.2 * Y[n-1] + 0.1 * Y[n-2]

and (shout out to Dan B and Matt L) using the scipy "lfilter" and "dimpulse" functions. When using initial values of zeros, they match.

        Y2     = sig.lfilter(    [ 0.5 , 0.3 ], [ 1, -0.2, -0.1], X )
        T3, Y3 = sig.dimpulse( ( [ 0.5 , 0.3 ], [ 1, -0.2, -0.1], 1 ) )
        
        for n in range( 20 ):
          print( "%4d  %10.5f   %10.5f %10.5f" % \
                 ( n, Y3[0][n].real, Y2[n].real, Y[n].real ) )

Here are the values.

   0     0.00000      0.00000    0.00000
   1     0.50000      0.00000    0.00000
   2     0.40000      0.00000    0.00000
   3     0.13000      0.00000    0.00000
   4     0.06600      0.50000    0.50000
   5     0.02620      0.40000    0.40000
   6     0.01184      0.13000    0.13000
   7     0.00499      0.06600    0.06600
   8     0.00218      0.02620    0.02620
   9     0.00094      0.01184    0.01184
  10     0.00041      0.00499    0.00499
  11     0.00017      0.00218    0.00218
  12     0.00008      0.00094    0.00094
  13     0.00003      0.00041    0.00041
  14     0.00001      0.00017    0.00017
  15     0.00001      0.00008    0.00008
  16     0.00000      0.00003    0.00003
  17     0.00000      0.00001    0.00001
  18     0.00000      0.00001    0.00001
  19     0.00000      0.00000    0.00000

The obvious way to get the FIR coefficients directly is to do the polynomial division.

$$ \begin{align} H(z) &= \frac{B(z)}{A(z)} \\ &= \frac{b_0 + b_1 z + b_2 z^2 ...}{ 1 + a_1 z + a_2 z^2 .... }\\ &= h[0] + h[1] z + h[2] z^2 .... \end{align} $$

So I did some searching and found numpy.polydiv( B, A ), but was disappointed it doesn't work the way I wanted. It stops at "whole values" instead of "calculating the fractional part".

I wrote a routine to do this (included here for anybody else's benefit).

import numpy as np

#=============================================================================
def main():

        B = np.array( [ 0.5 , 0.3 ] )
        
        A = np.array( [ 1, -0.2, -0.1] )
        
        print( B )
        print( A )
        
        Q, R = DividePolynomials( B, A, 15 )
        
        print( Q )
        print( R )

#=============================================================================
def DividePolynomials( ArgNum, ArgDen, ArgLength ):

        Q = np.zeros( ArgLength * 2, dtype=complex )  
        R = np.zeros( ArgLength * 2, dtype=complex )  
        S = np.zeros( ArgLength * 2, dtype=complex )  
        
        R[0:len(ArgNum)] = ArgNum
        
        for d in range( ArgLength ):
          rd = R[d] / ArgDen[0]
          
          Q[d] = rd
          
          S.fill( 0.0 )
          
          S[d:d+len(ArgDen)] = rd * ArgDen
          
          R -= S

        return Q[0:ArgLength], R[ArgLength:]

#=============================================================================
main()

Here is the output:

[ 0.5  0.3]
[ 1.  -0.2 -0.1]
[  5.00000000e-01+0.j   4.00000000e-01+0.j   1.30000000e-01+0.j
   6.60000000e-02+0.j   2.62000000e-02+0.j   1.18400000e-02+0.j
   4.98800000e-03+0.j   2.18160000e-03+0.j   9.35120000e-04+0.j
   4.05184000e-04+0.j   1.74548800e-04+0.j   7.54281600e-05+0.j
   3.25405120e-05+0.j   1.40509184e-05+0.j   6.06423488e-06+0.j]
[  2.61793882e-06+0.j   6.06423488e-07+0.j   0.00000000e+00+0.j
   0.00000000e+00+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j
   0.00000000e+00+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j
   0.00000000e+00+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j
   0.00000000e+00+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j]

The coefficients match nicely to the expected values that came from the impulse analysis and the remainder gives me an idea of how converged it is.

Of course, I did some searching and found this:

Is there a way to derive an FIR filter using an IIR filter?

In the linked question, the selected answer involved curve fitting, the other answers were consistent with what I was expecting. However, adding the criteria that you want to keep the filter order low, of course makes a better fit polynomial possible than a truncated $H(z)$. I didn't follow the paper references. IEEE papers are usually behind some paywall. But I see this as the identical math problem that we have had around here of "What is the best polynomial to fit $\sin(x)$ from $a$ to $b$" with the quotient of $B(z)/A(z)$ playing the role of the Taylor series.

• Question 1: Is there a polynomial division function I missed in numpy/scipy that does what I want. [Solved: See Olli's answer]

• Question 2: In "real life", what are typical FIR lengths for typical IIR to FIR conversions, and is this extra polynomial fitting step generally needed/beneficial?

I realize that I am dealing with a small rather well behaved IIR in my example.

Answer 1

You can still use NumPy's polydiv, if you first zero-pad B. In Python, after your numpy import and A and B initialization:

print(np.polydiv(np.pad(B, (0, 10)), A)[0])

In older NumPy versions, numpy.pad needs an additional parameter mode='constant', which was made the default value since NumPy 1.17. Running the above prints a sequence of numbers that is identical to what you got by the other means:

[5.00000e-01 4.00000e-01 1.30000e-01 6.60000e-02 2.62000e-02
 1.18400e-02 4.98800e-03 2.18160e-03 9.35120e-04 4.05184e-04]

Answer 2

Question 2: In "real life", what are typical FIR lengths for typical IIR to FIR conversions, and is this extra polynomial fitting step generally needed/beneficial?

That depends highly on what your IIR filter does. In my neck of the woods (audio) the answer is typically "a few thousand".

It really depends on the "frequency" resolution. At what frequency intervals does something "interesting" happen? Here is is a simple example: Let's say you have a 3rd order butter-worth highpass at 40 Hz sampled at 44.1kHz. Dividing the two gives 1000 which is in the ballpark. Turns out 1024 is pretty bad, 2048 is "ok" and 4096 is "good".

Speaking more formally: it really depends on the location of your poles. The lower the frequency and the higher the Q, the more FIR coefficients you need.

I don't think polynomial division helps much here. You either need to truncate the IIR impulse response (with potentially some windowing/tampering at the tail end) or do a straight FIR fit of the transfer function, where you can play around with the error criteria. Trying to match a specific IIR response may not be useful: cut out the "middle man" and design your FIR filter directly from the requirements

Answer 3

I am not sure why you are trying to approximate an IIR with an FIR, but an efficient way to do that is with Truncated IIR Filters.

Might be worth exploring.

If you know your 1st or 2nd-order IIR filters and want to truncate their impulse response, this shows you how to do it. Higher order IIRs can be factored into 2nd-order sections, but the finite lengths of the TIIR will add. If you want a higher order IIR with a given truncation length, the thing to do is to parallelize the 2nd-order sections (using partial fraction expansion) and then truncate each parallel 2nd-order section to the same common FIR length.