Infinite Impulse Response Filters

Digital filters which must be implemented recursively are called Infinite Impulse Response (IIR) filters because, theoretically, the response of these filters to an impulse never settles to zero. In practice, the impulse response of many IIR filters approaches zero asymptotically, and may actually reach zero in a finite number of samples due to the finite word length of digital computers.

One method of designing digital filters starts with the Laplace transform representation of an analog filter with the required frequency response. For example, the Laplace transform representation (or continuous transfer function) of a second order notch filter with the notch at f0 cycles per second is:

```
y(s)/u(s) = (f(0)/(2*pi) + s^2) / (1 + 2*s*f(0)/(2*pi) + s^2)

```

where s is the Laplace transform variable. Then the continuous transfer function is converted to the equivalent discrete transfer function using one of several techniques. One of these is the bilinear (Tustin) transform, where

```(2/d)*(z-1)/(z+1)
```

is substituted for the Laplace transform variable s. In this expression, z is the unit delay operator.

For the notch filter above, the bilinear transformation yields the following discrete transfer function:

```y(z)/u(z) = ((1+c^2)/2 - 2*c*z + ((1+c^2)/2)*z^2) / (c^2 - 2*c*z +
z^2)
```

where c = (1 - p*f0*d) / (1 + p*f0*d).

Enter the following IDL statements to compute the coefficients of the discrete transfer function:

```delt = 0.02
; Notch frequency in cycles per second:
f0 = 6.5
c = (1.0-!PI*F0*delt) / (1.0+!PI*F0*delt)
b = [(1+c^2)/2, -2*c, (1+c^2)/2]
a = [ c^2, -2*c, 1]
```

Alternately, you can run the following batch file to compute the coefficients:

```@sigprc13
```

See Running the Example Code if IDL does not find the batch file.

IIR Filter Implementation

Since an Infinite Impulse Response filter contains feedback loops, its output at every time step depends on previous outputs, and the filter must be implemented recursively with difference equations. The discrete transfer function

y(z) = (( b(0) + b(1)*z + ... + b(nb)z^nb ) / ( a(0) + a(1)*z + ... + a(na)*z^nb )) * u(z)

is implemented with the difference equation

y(k) = ( b(0)*u(k-nb) + b(1)*u(k-nb+1) + ... + b(nb)*u(k)

- a(0)*y(k-na) - a(1)*y(k-na+1) - ... - a(na-1)*y(k-1) ) / a(na)

An IIR filter is stable if the absolute values of the roots of the denominator of the discrete transfer function a(z) are all less than one. The impulse response of a stable IIR filter approaches zero as the time index k approaches infinity. The frequency response function of a stable IIR filter is the Discrete Fourier Transform of the filter's impulse response.

The figure below plots the impulse and frequency response functions of the notch filter defined above using recursive difference equations.

Enter the following commands at the IDL prompt to create the plot:

```; Load the coefficients for the notch filter:
@sigprc13
; Degree of denominator polynomial:
na = N_ELEMENTS(A)-1
; Degree of numerator polynomial:
nb = N_ELEMENTS(B)-1
N = 1024L
; Create an impulse U:
U = FLTARR(N)
U[0] = FLOAT(N)
Y = FLTARR(N)
Y[0] = B[2]*U[0] / A[na]
; Recursively compute the filtered signal:
FOR K = 1, N-1 DO \$
Y(K) = ( TOTAL( B[nb-K>0:nb]*U[K-nb>0:K]) \$
- TOTAL( A[na-K>0:na-1]*Y[K-na>0:K-1] ) ) / A[na]
; Compute spectrum V:
V = FFT(Y)
; F = [0.0, 1.0/(N*delt), ... , 1.0/(2.0*delt)]
F = FINDGEN(N/2+1) / (N*delt)
; Magnitude of first half of V:
mag = ABS(V(0:N/2))
; Phase of first half of V:
phi = ATAN(V(0:N/2), /PHASE)
; Create log plots of magnitude in decibels and phase in degrees.
; Set up for two plots in window:

!P.MULTI = [0, 1, 2]
PLOT, F, 20*ALOG10(mag), YTITLE='Magnitude in dB', \$
XTITLE='Frequency in cycles / second', /XLOG, \$
XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1, \$
TITLE='Frequency Response Function of b(z)/a(z)'
PLOT, F, phi/!DTOR, YTITLE='Phase in degrees', \$
YRANGE=[-180,180], YSTYLE=1, YTICKS=4, YMINOR=3, \$
XTITLE='Frequency in cycles / second', /XLOG, \$
XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1

!P.MULTI = 0
```

 Note
Because the impulse response approaches zero, IDL may warn of floating-point underflow errors. This is an expected consequence of the digital implementation of an Infinite Impulse Response filter.

Alternately, you can run the following batch file to create the plot:

```@sigprc14
```

See Running the Example Code if IDL does not find the batch file.

The same code could be used to filter any input sequence u(k).