2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I think scipy.signal.dimpulse is what you're looking for. $\endgroup$ – Matt L. Aug 1 at 8:01
  • $\begingroup$ @MattL. Thanks Matt. I am actually working on the general X case, not just an impulse, so this is a special case, but good to know and I'll include it. I figured if I can find the exact IIR to any X and Y then I can find a corresponding FIR, so full solution, compact IIR or extended FIR. In the process, I wanted to confirm the principle illustrated in my question, and couldn't find a polynomial (we used to call it synthetic division) routine. I thought that odd. So now it turns out if I can find any IIR in general I can also solve the best fit IIR to any FIR. Correct? $\endgroup$ – Cedron Dawg Aug 1 at 10:57
  • $\begingroup$ I'm a bit confused at your problem statement. In the above question you seem to be looking for a way to approximate a given IIR filter by an FIR filter. In your comment under Hilmar's answer, you seem to be looking for an IIR filter that produces a given output sequence from a given input sequence. Those two problems are different, aren't they? $\endgroup$ – Matt L. Aug 1 at 16:36
  • $\begingroup$ Yes, my question concerns a subset of my problem. I should be able to convert back and forth between a best fit IIR and a best fit FIR. My initial inclination was IIR->FIR easy, FIR->IIR hard, then I found the link and the IIR->FIR was made harder by the order constraint. I'm just trying to get an idea of the landscape of the arena I am in. I have heard of "thousands of FIR coefficients before", confirmed again by Hilmar. A truncated FIR->IIR is going to be harder than a longer one I would think. Numerically I can find a best fit IIR, exactly with the right rank set, but more freedom... $\endgroup$ – Cedron Dawg Aug 1 at 16:46
  • $\begingroup$ ... and a's and b's "swap values" for a better local fit. I think I need to ensure wider bandwidth in by test signal. I can also find and exact FIR fit for a FIR filter, but I don't think that is anything special. Neither might my IIR technique be. $\endgroup$ – Cedron Dawg Aug 1 at 16:50
1
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks, I am trying to "cut out the middle man" in the other direction, wondering if it had any practical value and the dimensions of real life systems. Suppose a system is governed by an unknown IIR. My shallow understanding (just starting to explore) is that the usual procedure is to find a FIR (H(z)) then use methods to derive the best fit IIR (A(z) and B(z)). I set about trying to find the best IIR directly from an arbitrary input and it's corresponding output. Some success so far. The polynomial division is part of a double check. $\endgroup$ – Cedron Dawg Aug 1 at 13:35
  • $\begingroup$ In my experience you never use a middle man. If you want FIR you design this directly, if you want IIR, you design that directly. There is an interesting middle ground: "Warped FIR" it's an IIR but can be designed using FIR design methods in a "warped" frequency domain. $\endgroup$ – Hilmar Aug 1 at 16:24
  • $\begingroup$ I'm trying to solve the opposite black box problem. Given a stretch of the input and and corresponding output of a filter, figure out whether it is IIR or FIR and what the coefficients are. These considerations are part of that, real life parameters will give me targets to test against. So, first I have to generate some ideal noiseless inputs and outputs to work against. The above code is part of that. There can be noise in the signals, vs noise in my data collection for real data. $\endgroup$ – Cedron Dawg Aug 1 at 16:39
1
$\begingroup$

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]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Okay, I put zeros in the "B =" line of my code and the resulting numbers do match. Thanks. Q1 solved. $\endgroup$ – Cedron Dawg Aug 1 at 14:31
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Trying to get a feel for ball park numbers. I am supposing some one is measuring something and modeling it with an IIR. Then they want to reproduce it over a fairly small range of frequencies, so a small FIR might be a better fit and safer bet for their operating region. What are typical real life, in your experience, solving real problems, sizes of filters and do you ever convert them back and forth like that? So, I appreciate the link, I'll remember where it is. Matt left a good one too. How many to model voice degradation in air over a distance of inches? $\endgroup$ – Cedron Dawg Aug 3 at 2:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.