Real-Time Digital Filters: Infinite Impulse ... - American Chemical Society

Digital. Filters: Infinite. Impulse. Response. Filters. Stephen E. Bialkowski. Department of Chemistry and Biochemistry. Utah State University. Logan,...
0 downloads 0 Views 10MB Size
A/C INTERFACE

Real-Time Digital Filters: Infinite Impulse Response Filters Stephen E. Bialkowski Department of Chemistry and Biochemistry Utah State University Logan, Utah 84322

Real-time digital filters fall into two main classes: finite impulse response (FIR) and infinite impulse response (IIR) (1-6). In part one of this series, digital FIR filters were discussed (2). In this A/C INTERFACE the design, operation, and interpretation of digital IIR filters will be discussed. Two main types of digital IIR filters will be examined: the frequency-domain filter, which uses the reiterative estimation process to result in digital filter equivalents of analog frequency-domain filters, and the time-domain IIR filter, which is used for the estimation of signals with known time-domain dependence. This article is by no means comprehensive and will probably only whet the appetite of those who are interested in designing specific IIR filters. IIR filters are indeed difficult to understand. But a little understanding can go a long way in allowing one to design and use this type of filter. 0003-2700/88/0360-403A/$01.50/0 © 1988 American Chemical Society

Basic to understanding both FIR and IIR filters is the impulse response function. The impulse response function is the time-dependent reaction of the filter to a short-duration impulse input. Effectively, the digital filter performs a convolution between the input and the filter impulse response, resulting in the filter output. Convolution of the filter impulse response with the measurement input should not result in some arbitrary output signal. Rather, it should result in a signal estimate with an enhanced signal-to-noise ratio. Thus the signal estimate is highly dependent on the impulse response of that filter. There are many types of filters and many functions that a filter can perform. Playing the odds by an arbitrary choice of filter will most likely result in poor signal estimation. If filter design is based on expected signal and noise components, then the optimum signal estimate can be obtained. In a typical measurement situation, a certain form of the signal is expected. Knowledge of the expected signal and noise components will allow the design of a filter impulse response that is most suited to the measurement application at hand. The impulse response of the FIR fil-

ters discussed in part one of this series was shown to consist of weighting functions with a finite number of terms in the vector representation. Because of the finite time domain of these filter functions, filtering could be accomplished by performing a convolution of the impulse response with the input data. However, it was pointed out that the digital FIR filters were poor estimators of low-frequency components of the signal. The reason for this is simple: The finite duration of the FIR filters results in a low-frequency limit at least equal to the inverse time duration of the impulse response function. This limitation is not present—at least in theory—for the digital IIR filter. As the abbreviation implies, digital IIR filters have infinite impulse response durations. Clearly, convolution of the measurement input with an infinitely long impulse response function, as was done for the FIR filter, cannot be accomplished with an IIR filter. How then does one perform convolution with an infinitely long filter impulse response in real time? The answer lies in how the output is calculated. Digital IIR filters are recursive in nature; both current measurement inputs and past outputs are used to esti-

ANALYTICAL CHEMISTRY, VOL. 60, NO. 6, MARCH 15, 1988 · 403 A

mate the current input. Thus the digi­ tal IIR filters are of infinite reaction or response duration because the current output is a function of all previous in­ puts. Two applications of the digital IIR filter for real-time estimation will be discussed below. The first is for per­ forming the digital equivalent of lowfrequency or wide-band analog signal filtering. This application is not strictly dependent on prior knowledge of the time-dependent character of the sig­ nal, but rather on its frequency compo­ nents. Knowledge of the frequency components of the signal and noise constituents of the measurement is typical. For example, in a gas chromatograph or in a scanning spectrophotometric instrument, the signals are known to vary smoothly over time and thus are of low frequency. Highfrequency components of these mea­ surements can be discarded with little, if any, loss of information. Why use the digital frequency filter if analog filters can be constructed? First, the digital frequency filter gener­ ally is more flexible and easier to con­ struct and control than is the analog circuit equivalent. Second, there are no sources of noise in the digital filter oth­ er than those caused by quantized arithmetic error. Third, the use of modern instrumentation often results in data that are in the digital format. The analog signal is just not available. The second application of digital IIR filters will address optimal filtering of signals of which prior knowledge of the experimental system is in hand. This particular filter, called the Kalman fil­ ter, is based on what is known as the system state description of the particu­ lar experimental system under study. Use of the Kalman filter requires an estimate of the initial state of the sys­ tem, the noise components, and the equations that describe how the system will evolve in time. Background When a digital computer is used to pro­ cess time-dependent analog data, an analog-to-digital converter samples the analog data at equal time intervals, T, thereby obtaining measurement data sequences of the form x[fcT]. The defi­ nition of the operation of a digital filter is ultimately dependent on the data representation. To summarize the data representation, the input to the digital filter is a discrete-time, periodically sampled measurement, x[feT], consist­ ing of the signal, s[feT], combined with noise, v[kT] (1): x[kT] = s[kT] + v[kT]

(1)

where the time is an integer factor, k, of the sample time interval, T. The out­ put of a digital filter is represented as the convolution of the impulse re­

sponse, h[t], with the measurement,

y[kT] = Y h[nT]x[(k - n)T] (2) η

However, because h [t] is infinite for the IIR filters, it is not practical to convo­ lute h[t] with the measurement (as in the FIR filter case). Rather, the current output of the filter, y[feT], is formulat­ ed based on the past filter output, y[(k — n)T], and the current measure­ ments. This filter is called a recursive estimator because the current output is recursively based on the previous out­ puts. Again, it is the recursive nature that allows the IIR filter to be realizable without the input convolution with an infinite duration filter function. Digital IIR filters may be represent­ ed by difference equations (4). Com­ prehensive understanding of these fil­ ters requires knowledge of both timeand frequency-domain processes. A general transform method for these digital filters is required. The Fourier transform may be used to relate the time- and frequency-domain filter functions, but it is best suited to differ­ ential equations and continuous func­ tions. To use the Fourier transform on sampled sequential data, a sampling function must be included along with each function definition. The discrete Fourier transform (DFT) is a compro­ mise. It does operate on discrete data sequences but is difficult to use with causal, one-sided data for which the Laplace transform is most utilitarian. Moreover, DFT transform pairs are both discrete in nature. The ideal transform method would operate natu­ rally on causal, sampled sequential data, would be easily implemented for the solution of difference equations, and would result in functional repre­ sentations of the frequency domain. The Z-transform is just such a method (4,5).

The Z-transform method will first be illustrated by application to a typical IIR filter equation. Consider the case of a recursive filter defined by the differ­ ence equation (4). y\m

= y[(k - 1)T] + K k (x[fcT]-y[(fe-l)T])

(3)

where Kk is the gain for the fcth esti­ mate. The stable range of Kk is between 0 and 1. For Kk = 1, the current output,

,y[&T], is the current input, x[feT]. For Kk = 0, the output is a constant equal to the previous output value. Because this filter is recursive, it is difficult to determine an impulse response by in­ spection of the difference equation. However, using the Z-transform meth­ od, determination of the impulse re­ sponse of this or any other discrete representation is relatively easy. First, one finds the transform of the filter equation. Second, the terms are rear­ ranged so that input and output trans­ formed variables are separated. All terms that are independent of the mea­ surement input and the filter output will make up the filter transfer func­ tion. The third step is to find an inverse transform that results in the convolu­ tion between the input and the impulse response. In particular, the inverse transform of the transfer function is the impulse response of the filter. The Z-transform provides a simple technique for converting difference equations of a sequence of numbers into an algebraic equation of complex functions. It is related to the DFT of causal data. The DFT of y[kT], which extends over a finite range of time, is defined by F\y[kT]] = Y(nil) Ν

= ^y[AT]e-' n S i * T

where F\ } denotes the discrete Fourier transform of the Ν + 1 point data, Ω = 2π/ΝΤ is the discrete angular frequen­ cy, i is the 0 6(k) &(k - m) af /ca* frV (*+1Xfc+2)... (k + m)a*/ml β-"* ko'""

F{z) =

ΣΗ*)*'*

1 z'm zl(z - a) az/{z — a)2 az(z + a)/{z - a)3 zm+V(z-a)m+1 z/(z - e_«) ze-"/(z - β~α)2

the filter. It is at this frequency that the filter is most responsive. The inverse is most easily found by examining the Z-transform pairs given above. In this case, h[kT] = Kkexp(-feT/r)

(10)

where the equation 1 - Kk = exp(-T/r)

(11)

defines the exponential time constant, τ. In fact, using a series expansion, one may find that τ « T/Kk. Equation 11

Figure 1. Illustration of how the filter system transfer function can be used to deter­ mine the algorithm for the real-time digital filter. The numerator terms of Hz) will act on the measurement Input, X(z), and the denominator will act on the output, Y[z). The Inverse ζ Is a delay operation. Using the transfer function In the filter equation allows the system diagram to be constructed, which can be used as a flow chart for computer algorithm develop­ ment. Alternatively, the Z-transform, which yields the explicit equation for the time-domain filter, can be obtained.

also implies that Kk must be positive and less than 1 for a realizable filter. It is important to note that the computa­ tionally simple act of feeding back a constant fraction of the past output to the current input results in a filter with an exponential impulse response. The larger the feedback gain, given by 1 — Kk, the greater the time constant. This impulse response is essentially in­ finite in duration and would be a timeand memory-consuming process using FIR filter algorithms. In the context of digital IIR filters, the Z-transform technique for the de­ termination of filter characteristics is more than just a mathematical conve­ nience. Because of the simplicity of the shift property, the transfer function can be used to determine the algorithm for the digital filter. In fact, the scheme for the entire algorithm is contained in the transfer function. To design an al­ gorithm, one simply reduces the form of the transfer function to one that is dependent only on inverse ζ powers. All terms in the numerator of this function will result in operations on the input measurements, and all terms in the de­ nominator will result in operations on the output. Consequently, transfer functions without inverse ζ-dependent denominators are strictly FIR filters.

406 A · ANALYTICAL CHEMISTRY, VOL. 60, NO. 6, MARCH 15, 1988

The determination of a digital filter al­ gorithm based on the transfer function of Equation 9 is illustrated in Figure 1. Determination of the algorithm is not difficult in this particular instance. But often one has a set of frequency-do­ main filter design characteristics that are to be implemented in a real-time digital filter. The Z-transform method is of particular value in these cases. Design of digital frequency fitters

Recursive filters have poles that can be placed anywhere within the Nyquist limited frequency range (4). Consider the filter given by Ν

^any[k-n]=x[k]

(12)

n=0

Here the sample interval, T, is dropped for simplicity because it is implied in the sequential data. There can be as many poles as there are different past estimates in the differential equation. Because there are several past esti­ mates in this differential equation, the resulting filter generally will be equiva­ lent to an analog multiple-pole filter. The Z-transform of this filter is Ν

^ αηζ~ηΥ(ζ) = X(z) n=0

(13)

and the transfer function is /

H(z) = l

Z

U = Χ

α ζ

Σ " ~"

(14)

I rc=0

Equation 14 can be expressed in terms of Ν products,

H(z) = CzN

Y[(z-zn)

n = 1,2,3... Ν

y[k] = x[k] - aiy[k -

(18)

Three cases can arise: The filter is un­ derdamped when αϊ 2 < 4a 2 , critically damped when αϊ 2 = 4a2, and overdamped when αϊ 2 > 4α2. For the criti­ cally damped case, both roots are equal and

(15)

(aj2)zl(z

(16)

The Z-transform is 2

+ a2z~ ) = X(z)

(17a)

It follows that H(z) is H(z) = 22/(22 + αλζ + α2) = [z/(z - Zl)][z/(z - z2)]

=

2

- a / 2 ) + z/(z - axl2)

h[k] = (k + 1MV2)*

where the two roots are determined from the solution to the quadratic equation,

(20)

Filter stability requires that αχ < 2. For small values of οχ, this impulse re­ sponse will resemble an exponential decay. Large values of oi result in a more complicated impulse response. The impulse response and the filter re­ sponse to a step function are illustrated in Figure 2. In the underdamped case there are two complex conjugate roots. The problem may be simplified by defining z2 + αγζ + a2 = (z - cei0)(z - ce" i0 ) (21) where a! = - 2 c c o s 0

(22a)

a2 = c2

(22b)

and

With these definitions, the system transfer function is H(z) = z2l{z - ceiB)(z - ce"1'6)

(17b)

(19)

The impulse response can be deter­ mined directly from the table of Ztransforms,

1]a&[k - 2]

Y(z)(l + aiz~

ι ± ι/(αι 2 - 4