Digital Signal Processing. The Scientist and Engineers Guide. Steven W. Smith
.pdf
Chapter 12The Fast Fourier Transform |
239 |
domain decomposition is stopped two stages early, when each signal is composed of only four points. Instead of calculating the last two stages, highly optimized code is used to jump directly into the frequency domain, using the simplicity of four point sine and cosine waves.
Another popular algorithm eliminates the wasted calculations associated with the imaginary part of the time domain being zero, and the frequency spectrum being symmetrical. In other words, the FFT is modified to calculate the real DFT, instead of the complex DFT. These algorithms are called the real FFT and the real Inverse FFT (or similar names). Expect them to be about 30% faster than the conventional FFT routines. Tables 12-6 and 12-7 show programs for these algorithms.
There are two small disadvantages in using the real FFT. First, the code is about twice as long. While your computer doesn't care, you must take the time to convert someone else's program to run on your computer. Second, debugging these programs is slightly harder because you cannot use symmetry as a check for proper operation. These algorithms force the imaginary part of the time domain to be zero, and the frequency domain to have left-right symmetry. For debugging, check that these programs produce the same output as the conventional FFT algorithms.
Figures 12-10 and 12-11 illustrate how the real FFT works. In Fig. 12-10,
(a) and (b) show a time domain signal that consists of a pulse in the real part, and all zeros in the imaginary part. Figures (c) and (d) show the corresponding frequency spectrum. As previously described, the frequency domain's real part has an even symmetry around sample 0 and sample N/2 , while the imaginary part has an odd symmetry around these same points.
4000 'INVERSE FFT FOR REAL SIGNALS |
|
|
4010 |
'Upon entry, N% contains the number of points in the IDFT, REX[ ] and |
|
4020 |
'IMX[ ] contain the real and imaginary parts of the frequency domain running from |
|
4030 |
'index 0 to N%/2. The remaining samples in REX[ ] and IMX[ ] are ignored. |
|
4040 |
'Upon return, REX[ ] contains the real time domain, IMX[ ] contains zeros. |
|
4050 ' |
|
|
4060 ' |
|
|
4070 |
FOR K% = (N%/2+1) TO (N%-1) |
'Make frequency domain symmetrical |
4080 |
REX[K%] = REX[N%-K%] |
'(as in Table 12-1) |
4090 |
IMX[K%] = -IMX[N%-K%] |
|
4100 NEXT K% |
|
|
4110 ' |
|
|
4120 |
FOR K% = 0 TO N%-1 |
'Add real and imaginary parts together |
4130 |
REX[K%] = REX[K%]+IMX[K%] |
|
4140 NEXT K% |
|
|
4150 ' |
|
|
4160 |
GOSUB 3000 |
'Calculate forward real DFT (TABLE 12-6) |
4170 ' |
|
|
4180 FOR I% = 0 TO N%-1 |
'Add real and imaginary parts together |
|
4190 |
REX[I%] = (REX[I%]+IMX[I%])/N% |
'and divide the time domain by N% |
4200 |
IMX[I%] = 0 |
|
4210 NEXT I%
4220 '
4230 RETURN
TABLE 12-6
240 The Scientist and Engineer's Guide to Digital Signal Processing
|
|
Time Domain |
|
|
|
Frequency Domain |
||
2 |
|
|
|
8 |
|
|
|
|
|
|
a. Real part |
|
|
|
|
c. Real part (even symmetry) |
|
|
|
4 |
Amplitude |
1 |
|
Amplitude |
0 |
|
|
|
|
|
0 |
|
|
|
-4 |
-1 |
|
|
|
|
|
|
|
|
|
|
|
|
-8 |
|
|
|
|
|
|
|
0 |
16 |
32 |
48 |
643 |
0 |
16 |
32 |
48 |
63 |
|||||||||||
|
|
|
|
|
Sample number |
|
|
|
|
|
|
|
Freqeuncy |
|
|
|
|
|||
2 |
|
|
|
|
|
|
|
|
|
|
|
|
8 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
b. Imaginary part |
|
|
|
|
|
|
|
|
d. Imaginary part |
(odd symmetry) |
|
|
|
|||||
|
|
4 |
Amplitude |
1 |
|
Amplitude |
0 |
|
|
|
|
|
0 |
|
|
|
-4 |
-1 |
|
|
|
|
|
|
|
|
|
|
|
|
-8 |
|
|
|
|
|
|
|
|
|
|
|
|
0 |
16 |
32 |
48 |
643 |
0 |
16 |
32 |
48 |
63 |
||||||||||||||||
|
|
|
|
|
Sample number |
|
|
|
|
|
|
|
|
|
|
Frequency |
|
|
|
|
|
||||
FIGURE 12-10
Real part symmetry of the DFT.
Now consider Fig. 12-11, where the pulse is in the imaginary part of the time domain, and the real part is all zeros. The symmetry in the frequency domain is reversed; the real part is odd, while the imaginary part is even. This situation will be discussed in Chapter 29. For now, take it for granted that this is how the complex DFT behaves.
What if there is a signal in both parts of the time domain? By additivity, the frequency domain will be the sum of the two frequency spectra. Now the key element: a frequency spectrum composed of these two types of symmetry can be perfectly separated into the two component signals. This is achieved by the even/odd decomposition discussed in Chapter 5. In other words, two real DFT's can be calculated for the price of single FFT. One of the signals is placed in the real part of the time domain, and the other signal is placed in the imaginary part. After calculating the complex DFT (via the FFT, of course), the spectra are separated using the even/odd decomposition. When two or more signals need to be passed through the FFT, this technique reduces the execution time by about 40%. The improvement isn't a full factor of two because of the calculation time required for the even/odd decomposition. This is a relatively simple technique with few pitfalls, nothing like writing an FFT routine from scratch.
244 |
The Scientist and Engineer's Guide to Digital Signal Processing |
is illustrated in Fig. 13-1, where various shapes of short input pulses produce exactly the same shape of output pulse. Second, the shape of the output waveform is totally determined by the characteristics of the system, i.e., the value and configuration of the resistors, capacitors and inductors. Third, the amplitude of the output pulse is directly proportional to the area of the input pulse. For example, the output will have the same amplitude for inputs of: 1 volt for 1 microsecond, 10 volts for 0.1 microseconds, 1,000 volts for 1 nanosecond, etc. This relationship also allows for input pulses with negative areas. For instance, imagine the combination of a 2 volt pulse lasting 2 microseconds being quickly followed by a -1 volt pulse lasting 4 microseconds. The total area of the input signal is zero, resulting in the output doing nothing.
Input signals that are brief enough to have these three properties are called impulses. In other words, an impulse is any signal that is entirely zero except for a short blip of arbitrary shape. For example, an impulse to a microwave transmitter may have to be in the picosecond range because the electronics responds in nanoseconds. In comparison, a volcano that erupts for years may be a perfectly good impulse to geological changes that take millennia.
Mathematicians don't like to be limited by any particular system, and commonly use the term impulse to mean a signal that is short enough to be an impulse to any possible system. That is, a signal that is infinitesimally narrow. The continuous delta function is a normalized version of this type of impulse. Specifically, the continuous delta function is mathematically defined by three idealized characteristics: (1) the signal must be infinitesimally brief, (2) the pulse must occur at time zero, and (3) the pulse must have an area of one.
Since the delta function is defined to be infinitesimally narrow and have a fixed area, the amplitude is implied to be infinite. Don't let this bother you; it is completely unimportant. Since the amplitude is part of the shape of the impulse, you will never encounter a problem where the amplitude makes any difference, infinite or not. The delta function is a mathematical construct, not a real world signal. Signals in the real world that act as delta functions will always have a finite duration and amplitude.
Just as in the discrete case, the continuous delta function is given the mathematical symbol: *( ) . Likewise, the output of a continuous system in response to a delta function is called the impulse response, and is often denoted by: h ( ) . Notice that parentheses, ( ), are used to denote continuous signals, as compared to brackets, [ ], for discrete signals. This notation is used in this book and elsewhere in DSP, but isn't universal. Impulses are displayed in graphs as vertical arrows (see Fig. 13-1d), with the length of the arrow indicating the area of the impulse.
To better understand real world impulses, look into the night sky at a planet and a star, for instance, Mars and Sirius. Both appear about the same brightness and size to the unaided eye. The reason for this similarity is not
246 |
The Scientist and Engineer's Guide to Digital Signal Processing |
Convolution
Just as with discrete signals, the convolution of continuous signals can be viewed from the input signal, or the output signal. The input side viewpoint is the best conceptual description of how convolution operates. In comparison, the output side viewpoint describes the mathematics that must be used. These descriptions are virtually identical to those presented in Chapter 6 for discrete signals.
Figure 13-2 shows how convolution is viewed from the input side. An input signal, x (t), is passed through a system characterized by an impulse response, h (t), to produce an output signal, y (t). This can be written in the familiar mathematical equation, y (t) ' x (t) th (t). The input signal is divided into narrow columns, each short enough to act as an impulse to the system. In other words, the input signal is decomposed into an infinite number of scaled and shifted delta functions. Each of these impulses produces a scaled and shifted version of the impulse response in the output signal. The final output signal is then equal to the combined effect, i.e., the sum of all of the individual responses.
For this scheme to work, the width of the columns must be much shorter than the response of the system. Of course, mathematicians take this to the extreme by making the input segments infinitesimally narrow, turning the situation into a calculus problem. In this manner, the input viewpoint describes how a single point (or narrow region) in the input signal affects a larger portion of output signal.
In comparison, the output viewpoint examines how a single point in the output signal is determined by the various values from the input signal. Just as with discrete signals, each instantaneous value in the output signal is affected by a section of the input signal, weighted by the impulse response flipped left-for-right. In the discrete case, the signals are multiplied and summed. In the continuous case, the signals are multiplied and integrated. In equation form:
EQUATION 13-1 |
|
%4 |
|
|
y (t ) ' |
m |
x (J) h (t & J) d J |
||
The convolution integral. This equation |
||||
defines the meaning of: y (t ) ' x (t )th (t ). |
|
|
||
|
|
&4 |
|
This equation is called the convolution integral, and is the twin of the convolution sum (Eq. 6-1) used with discrete signals. Figure 13-3 shows how this equation can be understood. The goal is to find an expression for calculating the value of the output signal at an arbitrary time, t. The first step is to change the independent variable used to move through the input signal and the impulse response. That is, we replace t with J (a lower case
















































































































































































