Digital filtering without the painA few years ago, while teaching an undergraduate course in electrical engineering, I needed to introduce the concept of digital filtering. The trick was to find a way to present the material without (at least initially) relying on advanced mathematical concepts. This is the approach I settled on.
Remember back in grade school when you learned about averages? They are a powerful tool. Most of the time, an average is defined as:
The idea of adapting averaging to digital filtering is that if you have n "noisy" measurements, the average of those measurements will give you the best "estimate" of the true nature of x. This technique seems like a low-pass filter, but in grade school that term wouldn't have made much sense.
Let's consider the simplest case of averaging readings. In this case you measure the signal at a fixed sampling rate. You could obtain the measurement data using an analog-to-digital converter (ADC) or a voltmeter. The only real requirement is that a fixed time interval occur between each reading. The simplest average you could perform would be:
The implementation is pretty simple. Just add the current measurement to the previous measurement and divide by the number of measurements, in this case simply two.
Moving average = FIR filter
This type of averaging is also a simple example of a symmetric finite impulse response (FIR) digital filter. Why is the response finite? Since only the latest two samples are used, any older samples are not available for use in the calculation. Therefore, the duration of any sample's influence on the output is finite. In DSP land, the equation might look something like this:
The subscripts on x and y simply indicate the relative age of the data and result.
A few comments are in order. First, notice that the coefficients are symmetric (both are 0.5); that is, the list of coefficients is identical when examined frontward or backward. This numerical palindrome doesn't seem like a big deal here with just two coefficients, but some symmetric FIR filters have many, many coefficients.
Also, observe that this filter has unity gain for a direct current (constant) input. After all, the average of two identical readings yields the same value as the readings themselves.
What might not be as obvious is that an input at half of the sample rate cancels out completely, as you can see in Figure 1.
Figure 1: Taken in pairs, the sampled input always sums to 0
Another important characteristic of this filter is that it's unconditionally stable. With no feedback of any kind, this stability must be the case. As you will see later on, not all digital filters can make this claim.
When you sampled the data, you knew what time the samples were taken. For what time is the output of the filter defined; in other words, what is the "delay" of the filter? Intuitively, since the magnitude of the filter output is the average of two sequential samples, it seems like the filter delay should be halfway between the averaged samples. For a signal at half the sample rate, this would be a delay of 1/4th the period of the input, or 90 degrees.
Verifying filter performance using Matlab
We used a little math and a lot of intuition to determine some of the filter performance. A tool like MATLAB can give us a little better picture. The MATLAB program shown in Listing 1 was used.
Listing 1: MATLAB program to determine filter performance
% Program SimpleFIR % Specify Transfer Function Coefficients b = [0.5 0.5]; % Numerator a = ;% Denominator (always 1 for FIR filters)
% generate the filter response points w = 0:pi/255:pi; [h, omega] = freqz (b, a, 512);
% generate the magnitude plot plot (omega/pi, abs(h)); grid; xlabel ('\omega/\pi') ylabel ('Linear Gain') title ('Simple FIR Magnitude Response'); pause
% generate the phase plot plot (omega/pi, 180*unwrap(angle(h))/pi); grid; xlabel('\omega/\pi') ylabel('Phase Angle (deg)') title('Simple FIR Phase Response'); pause
The magnitude plot shown in Figure 2 verifies that the gain is 1.0 at zero frequency and falls to 0.0 when the frequency is half the sample rate.
Figure 2: Simple FIR magnitude response
The phase plot in Figure 3 shows that the filter has linear delay that is a maximum of 90 degrees at a frequency of half the sample rate.
Figure 3: Simple FIR phase response
Manually computing symmetric FIR filter performance
If you're like me, you're suspicious of software tools until you do the math manually. For a simple FIR filter it's not too tough. In Equation 3:
the subscripts on x and y indicated the relative "age" of the incoming data and the filter output. Let's rewrite the equation for this filter using the unit time delay operator z. The z operator is very handy. It allows for mathematical analysis of operations that rely on sampling at constant time intervals. In a mathematical expression when data is multiplied by z, the exponent of z indicates the relative age of the data. While a general discussion of z transforms is beyond the scope of this article, simple applications of z transforms can be very useful. Using the z operator, Equation 4 can be put into the following form:
X(z) indicates a time series of data inputs into the filter; Y(z) indicates a time series of the filter outputs. In this case, the transfer function (the term that converts the input into the output) is defined as:
Cast into the usual DSP terminology, our objective now is to find the characteristics of our transfer function. Rather than spend a semester in a linear systems course, for now let's just pick a few facts about z transforms that we will need. First z can be defined as:
A note about j. No, you can't take the square root of -1, but think about the characteristics that j would have. When j is multiplied by j, the result is -1. Two terms containing j can be added together. So j is really just an operator that would behave like the square root of -1 would.
A slight digression here to talk about the maximum bandwidth of digital filters. The maximum output bandwidth of a digital filter is defined to be less than half of the sample rate and is referred to as baseband. If the sample rate is considered to be 2π, then defining the frequency ω to be between 0 and π should make perfect sense. Remember our objective here: find the filter response as a function of frequency.
The plan, mathematically speaking, is to get rid of the z terms in the equations as quickly as possible. You can do this using a little trig to rewrite Euler's Formula into the form:
Can we get our equation into this form? Yes, if we remember that both coefficients are identical and apply a little math step by step:
About this time you may look at Equation 9 and be tempted to say, "just great. And now what does that mean?" There are three terms to consider.
The constant term 2b0 is simply 2" 0.5 or just 1.
The next term:
has a constant magnitude of 1, but its phase varies linearly with frequency. This variance is why symmetric FIR filters are said to have linear phase performance. By the way, linear phase performance is a very good thing. It minimizes distortion of filtered data.
The last term:
has an amplitude that varies with frequency. When ω is zero, the amplitude is 1. When ω is equal to π, the amplitude is zero. This duplicates our understanding of the filter from intuition. The actual curve is shown in Figure 4.
Figure 4: Two point average magnitude response
As you can see it's pretty easy to predict the performance of a filter with just two "taps." As you might imagine, when the number of coefficients increases, the performance of the filter increases as well. Of course, from the point of view of filter design, we solved the problem backward in the previous example. Now let's try going in the other direction.
Manual design of an FIR filter
You can design FIR filters by hand. Let's add the constraint that when the input frequency is one-quarter of the sample rate, the gain of the filter is 0.5. We still want the filter to be a low-pass filter with unity gain at zero frequency and a gain of zero at half the sample rate.
Because we just added another constraint, we'll need to add another coefficient. Since we want to "pair up" terms, if you initially have an even number of taps, you can add just one coefficient. If you have an odd number of taps, you must add two. Since the original filter had two taps, we will need three; this makes the transfer function:
Since we're working with symmetric FIR filters, the first and last taps get paired up, so the coefficients must be identical. Working the math yields:
Again, the remaining exponential gives only the filter delay; in this case it is one time step. The amplitude of the filter is given by the second term. Some algebra will be required to calculate the two coefficients. The Cos term vanishes when ω = π/2. At that frequency we want the response to be 0.5, so b1 = 0.5. When the frequency is 0, we want the response to be 1, so b0 = 0.25. When plotted, this response yields the curve shown in Figure 5.
Figure 5: Three-tap FIR filter magnitude response
Two characteristics of any length symmetric FIR filter can now be understood. First, the delay through the filter increases as the number of terms (or taps) increases. The delay is half of the total delay inside the filter. Second, as the number of terms in the filter is increased, the filter can be specified at more and more points. Clearly, using this approach it's possible to design low-pass, high-pass, and band-pass filters using algebra. As you add coefficients, the amount of computation increases; eventually it's pretty unworkable to do the math by hand. This process has been automated in software; algorithms often applied are the Parks-McClelland and Remez exchange algorithms. One other comment: true "brick-wall" filters require infinitely many taps, but a symmetric FIR filter designed using several thousand taps can look pretty good indeed!
Designing an FIR filter using MATLAB
Let's design a filter that has an application in amateur radio. The Dovetron MPC-1000R radio teletype demodulator uses separate superheterodyne receivers for the "mark" and "space" tones with an intermediate frequency (IF) of 750Hz. The MATLAB program in Listing 2 produces a 750Hz IF filter using DSP instead of op-amp active filters.
Listing 2: MATLAB program to determine filter performance
% Program DoRemez % Linear Phase FIR BandPass Filter Design using remez
%%%%%%Problem parameters%%%%%% fsamp = 8000.0; % Sample frequency
% Band Edge frequencies f1 = 650.0; % lower stop band limit f2 = 700.0; % lower pass band limit f3 = 800.0; % upper pass band limit f4 = 850.0; % upper stop band limit fedge = [f1 f2 f3 f4];
% magnitude values (you really want 0 gain in stop bands and unity gain in passband) mval = [0 1 0];
% deviations (specify the ripples in each band) rs = 35.0; % stopband ripple in dB rsl = 10^(-rs/20); % linear stopband ripple rpl = 0.1; %linear passband ripple dev = [rsl rpl rsl];
%%%%%%Remez Order Estimation (how many coefficients?)%%%%%% [n, fpts, mag, wt] = remezord (fedge, mval, dev, fsamp);
disp ('estimated filter order'); disp (n);
%%%%%%Remez Algorithm%%%%%% % generate the numerator coefficients< b = remez (n, fpts, mag, wt)
% generate a plot w = (400.0/4000.0*pi):pi/255:(1100.0/4000.0*pi); [h, omega] = freqz (b, 1, w); plot ((4000.0*omega/pi), 20*log10(abs(h))); grid; xlabel ('Hz'); ylabel('Gain, dB') title('DoRemez Magnitude Response'); pause
The program generated the filter using 179 coefficients shown in Figure 6. Of course, anyone who uses analog filters would really hate this filter. The pass-band ripple seems high and the out-of-band rejection is lousy, at 30dB. You can change the ripple parameters and build a better filter demonstrated in Figure 7.
Figure 6: DoRemez magnitude response
Figure 7: DoRemez magnitude response
Of course this improvement in performance requires more coefficients. In this case 369 were needed. That means at least 369 additions and 185 multiplications must be performed 8,000 times per second to keep the filter up to date. Given that a radio teletype demodulator would require two of these filters and the rest of the code, perfect filters may run you out of CPU cycles. You should expect to compromise on filter length (and therefore performance) in real-world applications.
Perhaps there is a more computationally efficient way to implement a digital filter. This requires us to go back to the concept of averages.
"Optimized" average computation
Let's say that you want to compute the average of the latest 100 data samples. This could be done using the following expression:
Using this approach you have to hang on to the most recent 100 samples. To compute the new result, the oldest sample is thrown away and the newest sample inserted. This seems pretty inefficient, computationally speaking. Maybe you can get a similar effect without all the complexity.
If you think about this for a minute, only one percent of the result depends on the latest data sample. The remaining 99% depends on the previous samples. Why not try this instead:
This computation has a DC gain of 1, and the newest reading is weighted at 1% as desired. It's also much faster to compute and only the previous estimate and a new sample are required to update the filter. Because of the feedback (using the previous estimate), this is an infinite impulse response (IIR) filter!
Why is that the case? Well, think about what happens if you get one really bad reading. After each new sample the effect of the bad sample is only reduced by a factor of 0.99. In a real-world DSP implementation, the effect of a bad sample is not truly infinite, but it can be a very long time.
This example of a simple IIR filter is worth remembering. It works quite well when carefully implemented and can have a very long time constant. It makes a great variable-bandwidth lowpass filter; just two coefficients have to be adjusted to change the bandwidth.
The IIR filter transfer function
The equation for what we just accomplished could be written as:
With a little algebra we can arrive at the transfer function:
In a standard format, numerator coefficient vector of the transfer function is called "b"; the denominator coefficient vector is called "a". The following MATLAB program was used to plot the filter response:
% Program SimpleIIR
% Specify Transfer Function Coefficients b = [0.01]; % Numerator a = [1 -0.99]; % Denominator
% generate the filter response w = 0:pi/255:pi; [h, omega] = freqz (b, a, 512);
% generate the magnitude plot plot (omega/pi, 20*log10(abs(h))); grid; xlabel('\omega/\pi') ylabel('Gain (dB)') title('Simple IIR Magnitude Response'); pause
% generate the phase plot plot (omega/pi, 180*unwrap(angle(h))/pi); grid; xlabel('\omega/\pi') ylabel('Phase Angle (deg)') title('Simple IIR Phase Response'); pause
The magnitude response plot in Figure 8 shows a steep low-pass response. As you can see, this filter has impressive performance given the simple math. Can you combine the simple FIR and IIR filters that have been examined? Sure.
Figure 8: Simple IIR magnitude response
Combining (cascading) filters
The FIR transfer function we found basically was:
Equation 16 can be multiplied with the IIR transfer function:
resulting in a new transfer function:
This new first order IIR filter, shown in Figure 9, has the characteristics of both the previous filters.
Figure 9: Combined IIR magnitude response
The characteristics include the steep low-pass curve of the IIR filter with the 0 gain at half the sample frequency of the FIR filter. Since the new filter uses a previous filter output value, the filter must still be an IIR filter. Most practical IIR filters use a time series of both the input samples and the filter outputs. IIR filters don't have linear phase response as you can see in Figure 10.
Figure 10: Combined IIR phase response
Given that it's possible to realize Bessel filters using an IIR design, this limitation can be mitigated in most cases.
IIR filters and stability
Consider the IIR filter example I just described with a constant input amplitude of 1.0. What happens, if due to a numerical error, the 0.99 coefficient is actually 1.0?
As you can see in Table 1, the output is increasing with each sample. Given enough time, the estimated value will become quite large and eventually crash the process. It doesn't take much for an IIR filter to go unstable.
Since an IIR filter uses feedback, it can obviously be unstable if it isn't designed correctly. As you can see from the example, the filter must be implemented correctly as well. It is possible--and pretty simple--to determine if a given IIR filter will be stable in design:
1. Repeatedly multiply all of the terms in the denominator of the transfer function by z until all powers of z are non-zero. This yields a polynomial in z in the denominator.
2. Solve for z in the denominator polynomial.
3. All of the roots of this denominator polynomial (often called "poles") must be less than one in magnitude. This is usually stated as, "The poles due in the denominator of a transfer function must reside inside of the unit disk."
Finally, due to numerical errors, any given DSP system generally has a practical limit to the number of terms (taps) that are "safe" with an IIR filter. It can be as few as 15 or 20 denominator coefficients.
Given a practical limit on the number of terms in an IIR filter, there is a way around the problem if a high-performance IIR filter is needed. A long IIR filter can be broken into smaller filters using partial fraction expansion and cascaded. Or simply cascade shorter identical IIR filter sections until you get the desired performance. Both solutions are commonly used in real-world designs.
Software design of IIR filters
IIR filters can be designed to implement traditional filter characteristics, such as Butterworth, Bessel, Cauer, and Chebyshev filters. Because of the long data persistence in IIR filters, an IIR filter of modest length can have the same basic amplitude performance as an FIR filter with many more taps. Here is a MATLAB program that designs a Bessel filter similar to the FIR-based IF filter we designed before.
% Program DoBessel % Bessel Bandpass Filter Design
%%%%%%Program Input Parameters%%%%%% % frequencies in Hz fsamp = 8000.0; fcenter = 750.0; baud = 45.45; harm = 2.0; fband = baud/2.0 * (0.5 + harm);
W1 = fcenter - fband; W2 = fcenter + fband;
N = input ('Enter filter order (7)?');
%%%%%%Normalize and Prewarp%%%%%% % Determine the normalized angular bandedge frequencies wp1 = 2.0 * pi * W1 / fsamp; wp2 = 2.0 * pi * W2 / fsamp;
% Prewarp the digital edge frequencies wpa1 = tan (wp1 / 2.0); wpa2 = tan (wp2 / 2.0);
Wn = [wpa1 wpa2];
%%%%%%Calculate the analog bandpass filter%%%%%% [nump, denp] = besself(N,Wn);
%%%%%%Perform Bilinear Transform%%%%%% [numd, dend] = bilinear (nump, denp, 0.5);
%%%%%%Plot the Results%%%%%% wd = 0.3927:0.005:0.7853; hd = freqz (numd, dend, wd); md = abs(hd); pd= angle(hd)*180/pi;
plot((8000.0*wd/(2.0*pi)), pd); grid; xlabel('Hertz') ylabel('Degrees') title('Bessel Discrete Bandpass Phase'); pause
plot((8000.0*wd/(2.0*pi)), (20.0 * log10 (md))); grid; xlabel('Hertz') ylabel('Magnitude') title('Bessel discrete bandpass magnitude');
This program generates a transfer function with only 15 numerator coefficients and 15 denominator coefficients. That means the filter update requires only 30 additions and 30 multiplies. The filter performance, shown in Figure 11, is good considering the small number of coefficients.
Figure 11: Bessel discrete bandpass magnitude
Bessel filters are used because of the linear phase characteristics in the pass band. Note how linearly the phase changes from about 720 to 790Hz in Figure 12.
Figure 12: Bessel discrete bandpass magnitude
You need to be aware that the discontinuities are from "phase wrapping" in the plot. The nonlinear portion is best seen from 600 to 700Hz and 800 to 900Hz. With any luck, there won't be much energy passed by the filter in this region.
Back to class
By the way, the final project in the class was to design a complete DSP equivalent of the Dovetron MPC-1000R demodulator. By the time the students finished coding their projects and debugged them, they had FSK modulation and demodulation as well as radio teletype transmission and reception pretty well figured out. More than one expressed an interest in amateur radio!
Gary Geissinger is chief electrical engineer at DigitalGlobe Inc and teaches senior- and graduate-level hardware and software interface classes at the University of Colorado. You can reach him at firstname.lastname@example.org.
Oppenheim, A. V. and R. W. Schafer. Digital Signal Processing. Prentice-Hall, Inc. 1975.
Rabiner, L. R. and B. Gold. Theory and Application of Digital Signal Processing. Prentice-Hall, Inc. 1975.
Stearns, S. D.and R.A. David. Signal Processing Algorithms. Englewood Cliffs: Prentice-Hall, Inc., 1988.
Mitra, S. K. Digital Signal Processing: a computer-based approach. McGraw-Hill, 2001.
Rorabaugh, C. B. Digital Filter Designers Handbook. McGraw-Hill, 1993. Lathi, B. P. Linear Systems and Signals. Berkeley-Cambridge, 1992.