Christopher E. Reid is a senior software engineer at Cadence Design Systems, formerly of MITRE Corporation. He is modeling signal integrity problems on printed circuit boards for Cadence. At MITRE he developed speech-recognition systems, and designed and simulated digital receivers for operation during local thunderstorms. He has extensive experience building software for digital signal processing.
Thomas B. Passin is on the technical staff at MITRE Corporation. His work includes simulations of the effects of an electromagnetic pulse on various antenna/receiver combinations. He developed new techniques for this work which have been verified by experiment. Before joining MITRE, Mr. Passin spent ten years designing and building various detectors of nuclear readiation.
The cornerstone of digital filtering is the digital simple harmonic oscillator, expressed by a two-term recursion formula
yt+1 = ayt + byt-1 + cxt+1where xt is the sequence of input values and yt is the sequence of output values, and a, b, and c are constants. This simple formula can produce a long clear bell tone of any desired frequency, or a thudding drum like noise, or even a the sound of a cheap tin whistle depending only on the values of the constants and the input values. This article explains why this simple formula works, shows how to calculate the frequency response of this simple filter, and introduces the methods required to implement a filter with any desired frequency response. You will be able to test these concepts yourself if you have a sound board on your system.
The Oscillator
The general equation for a digital simple harmonic oscillator, shown in Equation 1, is surprisingly simple. Given a sequence of input points, xt, and three constants, a, b, and c, the output points yt are defined by the two-term recursion formula. This equation can oscillate, and, from another point of view, it is also a simple filter. We call it a filter because given an input signal, xt, it produces an output signal that depends on the input in certain predictable ways. (We explore this point of view in the discussion of the ocsillator as a filter, later in the article.)In words, Equation 1 says the position of the oscillating body at time yt+1 depends on the applied force at that instant, xt+1, and on two previous position values. The time value, t, is an integer. This equation is the discrete version of a physical oscillator such as a pendulum, a simple bell, or a tuning fork. This article implicitly assumes the time between sample points is Dt = 1.
The impulse response of any system is its behavior after being hit with an ideal, instantaneous force. A quick blow of a stick on a drum, or a quick, staccato note on a piano approximates an impulse on each of those systems. In the digital world the ideal impulse is truly attainable. It is a driving force, xt, which is zero at all times except the instant of the impulse. As a matter of convenience later, we assume the impulse occurs at time t = 1, so the impulse is defined by
The function in Listing 1 calculates N points of the impulse response of Equation 1 and writes the output to a disk file. To see what is happening, examine the output of this algorithm for the special case
so the parameters are
Calculating a few points by hand reveals that the outputs are
and then the pattern repeats forever. This can be expressed by noting that, for this example
Of course this simple result is not an accident. The values for a and b were carefully chosen for this example. However, the result is typical. In this special case, the waveform is a sinusoid with a period of eight data points. The data rate on a CD player is 44,100 samples per second. That means that if this waveform was stored on a CD, you would hear it as a pure tone with the frequency 44100/8 = 5812.5 Hz (cycles per second), which is approximately four octaves above middle-E on the piano, just off the top of the keyboard.
It is instructive to play around with this algorithm. Try it for various values of a and b. Notice that c is simply a scale factor. For some combinations of a and b, the waveform amplitude decays steadily to zero. For other values, you will see it grow until the computer signals a floating-point overflow. You may find some values for which the waveform does not oscillate at all, but climbs steadily upward, maybe turning back down, maybe not. If you have a sound board you can even listen to the output, although you will first have to scale the waveform appropriately and use only integer values. Figure 1 shows the output of this function for the equation
yt+1 = 1.6019yt - 0.9802yt-1 + xt+1
Calculating the Coefficients
So what is the relationship between the coefficients of Equation 1 and its impulse response? A full investigation of that question takes more space and time than this article allows. For the whole story, see chapters one through three of our book [2], Signal Processing in C. The final result is as follows. It can be shown that the equation
yt+1 = 2e-a cos(w)yt - e-2ayt-1 + cxt+1has impulse response
In other words, when plucked, the system produces a damped sinusoid of frequency w radians per data point with a damping factor of exp(-a) per data point. If the damping is slow (a close to zero) this waveform sounds like a long clear bell tone. If the damping is quick (a large) it sounds more like a drum beat, or a dull thud. You can obtain a waveform of any desired frequency, say f Hz, by using the relationship
where s is the number of samples per second produced by your D/A converter. Equation 3 uses a = 0.01, and w = 0.2p, which produces the damped sinusoid of Figure 1.
It is worthwhile to examine some special cases of Equation 5. First, if a = 0, the waveform does not decay, but remains the same amplitude indefinitely. If a > 0, the waveform decays more or less quickly towards zero. If, however, a < 0 then the waveform grows exponentially until your computer signals a floating-point overflow.
As long as a >= 0 the impulse response is well behaved no matter what the value of w. There is one special case to consider, w = 0. In that case, Equation 5 does not seem to work. However, the recursion equation 4 is well defined for this case. Its impulse response can be obtained from Equation 5 by taking the limit as w approaches 0, which gives the proper impulse response in this case,
This waveform does not oscillate. Figure 2
plots this function for the case a = 0.01, the same damping used for .
The Oscillator as a Filter
The function in Listing 2 takes N input points and produces N output points using Equation 1 as a filter (with c = 1). For example, Figure 3 shows the waveform produced by Equation 3 when driven by white noise using this algorithm. When played through a D/A converter this waveform sounds somewhat like a cheap tin whistle.
Filters are not very useful without some method of designing them to give the frequency response required for a given application. The first step of that process is calculating the frequency response of Equation 1. It is not particularly difficult to derive the required formulas, but it would take quite a few pages. (For the details see chapter 14 of [2].) The formula for the gain of Equation 1 is simple to remember by writing the formula with all three output terms on the same side like this,
yt+1 - ayt - byt-1 = cxt+1The left-hand side of this equation looks like a polynomial, and the resemblance is more than superficial. The gain of this filter as a function of frequency w is given by
This gain is given in terms of power (voltage squared) because that is the quantity of primary concern.
The appearance of the complex exponential in Equation 7,
ejw = cos(w) + j sin(w)implies the gain is the same for any frequency Q = w + 2pk for any integer k. That means that the gain only needs to be calculated for some fixed range of w of length 2p. Conventionally the frequency range considered is -p to p. Because only the squared-magnitude of the complex polynomial is used in Equation 7, we also have
In other words, the gain at -w is the same as the gain at w. The solid line in Figure 4 plots the gain of equation 3, the same equation used for Figure 1 and Figure 2. The dashed line in Figure 4 is for an oscillator with the same frequency, but with damping term a = 0.05 instead of 0.01 in Equation 3. The peak of the gain spreads as the damping a increases.
A More General Filter
The simple harmonic oscillator is only the simplest possible equation of an infinite family of possibilities. The general recursion formula for this type of filter is given by (with a change of notation from Equation 1)
yt+1 = -c1yt - c2yt-1 - . . . - cnyt-n+1 + kxt+1There can be any number of fixed, real coefficients for this type of equation. The gain of this equation as a function of frequency is, in a way, no more difficult to calculate than the special case already considered. The equation can be rearranged as before,
yt+1 + c1yt + c2yt-1 + . . . + cnyt-n+1 + kxt+1which makes the left-hand side resemble a polynomial. Just as in the simpler case, the gain of this equation can be written as
It is possible to approximate any desired frequency-response curve to any specified accuracy using this polynomial. There are, however, practical problems in some cases. These difficulties are typical of numerical calculations in general, and will not be discussed in this article. In practice, it is necessary to test your implementation to make sure the frequency response is as desired. Special design techniques exist for the most commonly desired filter types. Park's book [1] is an excellent reference in this area.
Filters and Circular Buffers
The code of Listing 1 and Listing 2 was sufficient for demonstration purposes, but something more is required to implement the more general filter of Equation 8. The simple harmonic oscillator uses a two-term recursion formula which means it had to keep track of two previous output values. The general formula, Equation 8, requires the handling of n previous output values. The simple procedure of Listing 2 is not efficient when generalized to these more complicated filters. Furthermore, the previous code could filter only one vector of input values. More often it is necessary to keep the set of previous output values for use on the next vector to be filtered. Physically, these values represent energy stored in the oscillator. Throwing them away and starting over again would be like turning off your television every time you changed channels. A general purpose algorithm must save these values.The best way to store the old filter outputs is in a circular buffer. In general, a circular buffer can be any length, N. The array index, i ranges from 0 to N-1, and then wraps around back to 0. Mathematically, the array index is calculated modulo N. Such a circular buffer is exactly what is needed to hold the history of n previous output values of the filter, Equation 8. We can simply use a circular buffer of length N = n, store the filter outputs in order, and keep an index to the most recent value. All the values can then be accessed by incrementing the index from that initial value modulo N. Once a new filter output is calculated, that value overwrites the oldest previous output value in the buffer, and the index to the most recent value is adjusted accordingly.
This technique does work, but there is one problem. Every time an index is calculated, code that checks its value and calculates a new value modulo N when required must be executed. This is not efficient. It requires one, maybe two comparisons, a jump of the program counter if nothing needs to be done, and an addition, subtraction, or modulus operator if the value must be adjusted. This is a significant problem because this code will be executed n times for each output value of the filter calculated.
The process is much more efficient if N is chosen to be a power of 2, N = 2m. Of course, this often requires allocating a larger block of memory for the circular buffer than would otherwise be required, but this is usually not a problem. The key to writing more efficient code is the relationship
i & M = i modulo N, with M = 2m - 1where i is any integer and & is the masking operator (bitwise AND) of C. Thus, for example, the code to increment i by one and calculate its value modulo N reduces to
i = (i + 1) & MThe masking operator is very fast on most machines, so this code is very efficient.
To see why this masking operation is the same as calculating i modulo N = 2m, remember that the binary representation of M = 2m - 1 is m ones, filling bits 0 through m-1 of the integer representation. Thus any part of the integer that is a multiple of 2m is masked to zero (bits m and higher in the integer). This also works for negative numbers, as some reflection along the same lines will reveal. All this assumes the processor uses 2's complement binary arithmetic. There are very few processors that do not meet this requirement.
Listing 3 shows the structure definitions for the circular buffer and the filter routine. The circular buffer contains a member, index, which is used as the index into the buffer member where the next value should be placed. The mask member is the mask value as discussed above. It is assumed that this structure has been initialized with a length equal to a power of 2, the smallest power of 2 needed to contain the required output values. Usually the buffer is initialized to hold zeros.
The Filter structure in Listing 3 contains a pointer to a CircularBuffer and an array of filter coefficients. These coefficients are the {c1, c2, ..., cn} of Equation 8. This structure is used by the filter routine shown in Listing 4. The code is a straightforward implementation of Equation 8 except for the circular buffer which has already been explained. Notice that a local copy of the circular buffer structure is constructed. A glance at the structure definition in Listing 3 shows only two integers and one pointer are copied. If integers and pointers are both four bytes that means only 12 bytes are moved. This little bit of overhead saves execution time later by removing one level of indirection to retrieve these values from inside the loop.
The simplified code presented in this article should be enough to get started using digital filters. More general code requires a better way to handle vectors of data, and more flexible filtering routines that allow other filter techniques not discussed here. Reference [2] gives a full account of these considerations. Further, a source code library of DSP routines is available as well as an interactive DSP environment for testing many DSP concepts. Contact Reid Associates Companion Software, (800) 374-7343.
References
[1] T. W. Parks and C. S. Burrus. Digital Filter Designs. Wiley-Interscience, John Wiley and Sons, Inc., 1987.[2] Christopher E. Reid and Thomas B. Passin. Signal Processing in C. John Wiley and Sons, Inc., New York, 1992.