Consider an experiment. It could be an automotive crash test, an atom smashing in an accelerator, or a new surgical procedure. The common factor is that you get only one shot at collecting data, so the data acquisition (DA) must be flawless. You’ve shielded the signal leads and made sure that proper grounding procedures are used. You’ve placed anti-aliasing filters at the inputs of your DA cards, which are themselves of the highest caliber. You’ve tested everything in a dry run before the big event. At the appropriate time, you press the start button, and—the resulting data are contaminated with noise.

Back in the old days, when data were logged on paper with strip-chart recorders, you would be out of luck. But because the collected data—noise and all—now reside in machine-readable format on a PC, there is an escape clause called a digital filter.

A digital filter is a mathematical construct that takes a sequence of numbers and accentuates or attenuates certain of its characteristics. In some ways, digital filters are similar to their analog counterparts. This makes it possible for you to filter data after they have entered the digital domain. In many cases, digital filtering permits processing that would not be practical with an analog filter. One example of a complex filter would be a high-end audio equalizer that adapts its response to accommodate the acoustic characteristics of your living room. While it is possible to implement a filter to perform this function with analog circuitry, it’s easier and more economical to implement it as a segment of computer code executing on a DSP. Performing filtering functions numerically provides the additional benefit of consistency. A digital filter will always process a given signal in the same way because component tolerances, drift, and aging don’t affect numerical algorithms as they do analog components.

Numerical Simulations

One of the more intuitive methods of developing a digital filter is by numerically simulating its analog counterpart. Figure 1 shows a resistance-capacitance (RC) low-pass filter. Such a circuit can be readily simulated on a computer with a circuit analysis program (e.g., Simulation Program with Integrated Circuit Emphasis, or SPICE). By representing current and voltage as numerical quantities, the simulated circuit will behave analogously to the real thing.

Figure 1. An RC low-pass filter (A) can be numerically simulated by deriving the continuous-time differential equation (B) describing its behavior and then approximating it with an iterative equation (C) that can be solved numerically by a computer.

The basis of numerical simulation lies in converting the continuous differential equations that describe the behavior of an actual circuit into discrete time-difference equations. The general procedure for doing this is similar to differential calculus, but in reverse. In calculus, you begin with a small but finite Dt and move toward an infinitesimal dt, with the aid of much hand waving. To turn differential equations into difference equations, you have to translate the problem from the dt of calculus to the Dt of numerical analysis, also aided by a few gestures. Such circuit simulators as SPICE perform the conversion automatically. For a simple circuit, such as the low-pass filter of Figure 1, it’s also easy to do this manually. The differential equation relating output voltage y(t) to input voltage x(t) is:

In this case, you substitute D’s for d’s, yielding:

Notice that x and y are now just sequences of numbers, with members of the sequence indicated by subscripts, as opposed to continuous functions. The term Dy/Dt represents the amount of change in the variable y between discrete “ticks” spaced t units of time apart. A bit of rearrangement and sleight of hand will give you an equation that is useful for predicting the next value of y from both its present value and the present value of the input x:

Figure 2. Numerical models attempt to simulate reality as a sequence of numbers. The step response of a digital low-pass filter will be a series of discrete points, not a continuous function.

A difference equation in this form is easily expressed for computer solution. In BASIC, the equation might be written:

Y_new=y_old+(1/(R*C))*(x-y_old)*dt

You can also wrap some code around the statement to make it a generally useful function:

STATIC Function lo_pass(V_in, R, C, dt)

DIM Y

’ return present value of ‘Y’

Lo_pass=y

Y=y+(1/(R*C))*(x-y)*dt

End function

The STATIC keyword defines all the variables in the function to maintain their values between function calls, allowing the function to “remember” the state of the filter.

Figure 3. Another way of visualizing linear difference equations is by flow graphs built from delay units, adders/subtractors, and fixed gain stages. This graph represents a numerical model of a low-pass resistance-capacitance filter.

If this routine is called iteratively in a loop, you get behavior much like that of an analog filter. The following program simulates the low-pass filter with R = 1 ž and C = 1 F.

Dt=0.1

For t=0 to 2 step dt

If t>=0.5 then x=1 else x=0

Y=lo_pass(x, 1, 1, dt)

Print t, x, y

Next t

The input is a step from 0 to 1, with the step occurring when t = 0.5 s. The response of this filter can be seen in Figure 2. The input and output are represented by a series of points, a reminder that you are no longer working with continuous functions.

Difference equations can also be represented graphically (see Figure 3). Linear difference equations, the type most commonly used in digital filters, can be represented using only three types of function blocks: adders, gain stages (multiplication by a constant), and delay units. Because the output value of the delay updates only at discrete times (every Dt tick), it can be easier to think of these blocks as S/H de vices rather than as delays.

Because our numerical simulation of an RC filter is just an approximatio

Figure 4. Under certain situations, a digital filter can exhibit qualitatively different behavior than that of the analog filter on which it is based. By making the simulation time step too large, you can make a simulation of a resistance-capacitance low-pass filter become unstable and numerically “explode.”

n of the real thing, there are situations where that ap proxi mation will break down. One such case is when the value of Dt >2RC. Under this condition, the digital filter will initiate ever-increasing oscillations in response to any input (see Figure 4). This behavior—referred to as numerical instability—is qualitatively different than anything you would see in a real RC filter. Despite numerous methods of reducing its presence, numerical instability is both a common and an annoying feature of numerical simulations, often appearing in more subtle ways than runaway oscillation. Numerical instability can also serve as a reminder of the difference between a virtual world and the real world.

One variation of the notion of digitally simulating a filter is to forget the notion of explicitly accounting for Dt and just assume that the interval is a unit delay. This lets you write a generalized equation relating the values of the output sequence y_n to the values of an input sequence x_n. This equation is known as the generalized infinite impulse response (IIR) filter equation:

This equation gives you the freedom to define many types of filters, including low-pass, high-pass, bandpass, and band-stop. Higher order filters with specific characteristics can also be implemented. The a_i and b_k coefficients are determined by the characteristics of the filter desired and the value of Dt. How you derive these coefficients is a subject considerably beyond the scope of this article.

Because you limited yourself to building filters with linear functions (addition, delays, and multiplication by a constant), the resultant systems are linear. An important consequence of linearity is that it doesn’t matter if you first filter a number of input signals and add them, or first add them and then filter the sum. This property can be represented:

Another consequence of linearity in an IIR filter is that the filtering can be described by what is called the impulse response (see Figure 5), or the sequence of numbers that comes out of the filter in response to a single 1 fed into the input (the impulse). Because of linearity, all input signals can be viewed as the sum of a number of successive im pulses of varying magnitude. The impulse response therefore completely describes the behavior of a given filter.

The reason the IIR filter is called an infinite impulse response filter is that an im pulse at the input can result in a non-zero output forever, as numbers reverberate around numerical feedback loops. For practical purposes, however, a numerically stable filter will have an impulse response that dies down below a given threshold in a finite number of time steps.

A Filter an Accountant Could Love

Figure 6. The finite impulse response filter lets you build a filter function from a unit impulse response. As the unit impulse works its way through the delay line, it forces an output equal to the gain at that delay stage.

If you can obtain the impulse response of a desired filter, there is a conceptually simple computational structure to implement it. If you create a tapped line of delay elements (see Figure 6), you can create an arbitrary impulse response by simply setting the gains on each tap to the value of the impulse response at the corresponding time. As an impulse moves through the delay line, it will set the output to the value of the gain at that tap.

Because you can’t implement an infinite number of delays and taps, you can’t create infinite impulse responses. For most useful filters, however, the impulse response dies out after a number of time steps; if the response is truncated after several points, it can be realized through this structure. Because the filter will have an impulse response limited in duration to a finite number of delay stages, it is called a finite impulse response, or FIR filter. The equation for an FIR filter with K taps is:

A BASIC function to implement the filter can be written:

Static Function FIR1(x, K, a())

Dim xp(0 to MAXPOINTS)

For I=K to 1 step -1

Xp(I) = XP(I-1)

Next I

Xp(0) = x

Sum=0

For i=0 to k-1

Sum=sum+a(k)*xp(k)

Next I

FIR1=sum

End function

The number of delay stages is given by K, and the array containing the FIR coefficients is given by the array A().

Other than conceptual simplicity, the FIR filter offers the advantage of guaranteed stability. Because there are no feedback paths, there are no opportunities for the filter to numerically explode.

The FIR filter, however, is not without drawbacks. Although neither IIR nor FIR filters can meaningfully process signals with a frequency greater than half the sampling rate (the Nyquist limit), an IIR filter can exhibit interesting behavior at arbitrarily low frequencies because the filter’s feedback path provides it with a potentially infinite memory of past signal inputs. Because an FIR filter has no feedback, its memory is limited to the total number of delays. For an FIR filter to have response variations set at widely variable frequencies, it needs a lot of taps—and therefore lots of computation.

The simplest case of FIR filtering is when the weight of all n taps is set to 1/n. This results in what is known as a moving average filter (see Figure 7). These filters are often used to smooth data because they’re easy to understand and easy to implement. One application for this filter is in smoothing daily stock market prices. By averaging a day’s closing price with the closing prices of the two preceding days and the two following days, you get a clearer picture of a stock’s price trend.

This application illustrates an important notion of filtering—using future data to develop the present output value. A five day moving average filter that relies solely on past and present stock data would show a trend line that lags the data by two days, as shown in Figure 8B. Allowing the use of filter data results in more useful and representative output.

If you expand the definition of the FIR filter to allow for the use of future input points, it becomes even more powerful. While you obviously can’t obtain values of inputs before they exist, you can often accept a delay in the output data. In the stock market example, the delay is two days—too long to get rich quick. For many signal processing applications, however, especially those operating on already-collected data, a delay of quite a few sample periods can be perfectly acceptable.

For example, when digitally processing an audio signal from a CD player, even a significant delay (e.g., 100 ms), will be undetectable by the listener. In this situation, the person playing the CD has no idea what is physically being read from the disk at any point in time; a small delay from read head to speakers will not be noticed.

Ad Hoc Filtering

All the techniques discussed in this article have described linear filters. These filters are popular not only because they are useful but because their behavior is analog-like and can be analyzed rigorously. When your goal is to remove unwanted signals from a data set, nonlinear techniques (e.g., curve fitting) can be viewed as ad hoc filters. Depending on the characteristics of the data, ad hoc filters can often be constructed to clean up a dirty signal better than a linear filter of comparable complexity.

One data filtering scenario that I’ve often encountered has been that of false readings from A/D converters on DA cards (especially the really inexpensive ones). When reading data from a DA card, the signal ends up consisting of three superimposed features:

• the actual signal

• a few LSBs of noise

• a few wild readings resulting from misconversions

Although simple linear filtering techniques (e.g., a moving average filter) are generally adequate to suppress small amounts of noise, they don’t work as well at eliminating the effects of spurious 0 or full-scale readings. One technique I’ve found to work is to simply discard the readings identified as bad and interpolate from adjacent readings.

Figure 9 shows a filter performing this function. This filter reads five successive samples, sorts from lowest to highest, discards the high and low samples, and averages the remaining three. By ignoring exceptional samples, you have a filter that can deal much more effectively with random spikes than a simple three-point moving average filter (see Figure 10). A BASIC function to implement this filter can be written:

STATIC FUNCTION Best3of5 (x)

DIM a(5) b(5)

‘ Shift Samples into delay line

FOR i=5 TO 2 STEP -1

a(i)=a(i-1)

NEXT i

a(1)=x

‘ transfer to temp array & sort

FOR i=1 TO 5

b(i)=a(i)

NEXT i

‘ Bubble sort ok for only 5 items

FOR i=1 TO 4

FOR j=1 TO 4

IF b(j)>b(j+1) THEN SWAP b(j),

B(j+1)

NEXT j

NEXT i

‘ Take average of middle 3 values

Best3of5=(b(2)+b(3)+b(4))/3

END FUNCTION

There are certainly situations where this approach would be disastrous, such as when spurious spikes are actually part of the signal of interest or when the spikes occur frequently. As with any filtering technique, your mileage may vary depending on the characteristics of the signals you’re trying to process.

Conclusion

Digital filters are mathematical constructs that are useful for processing signals after they have been sampled through a DA system. Although many digital filters are based on approximate simulations of analog filters, the flexibility provided by numerical methods also permits processing techniques unique to the digital domain and often difficult or impossible to implement with traditional analog filters.

Home | Contact Us | Advertise © 2009 Questex Media Group, Inc.. All rights reserved. Reproduction in whole or in part is prohibited. Please send any technical comments or questions to our webmaster.