Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Digital Signal Processing. The Scientist and Engineers Guide. Steven W. Smith

.pdf
Скачиваний:
167
Добавлен:
02.05.2014
Размер:
7 Мб
Скачать
☆

Chapter 12The Fast Fourier Transform

237

Likewise, when this correct symmetry is present in the frequency domain, the Inverse DFT will produce a time domain that has an imaginary part composes of all zeros (plus round-off noise). These debugging techniques are essential for using the FFT; become familiar with them.

Speed and Precision Comparisons

When the DFT is calculated by correlation (as in Table 12-2), the program uses two nested loops, each running through N points. This means that the total number of operations is proportional to N times N. The time to complete the program is thus given by:

EQUATION 12-1

 

DFT execution time. The time required

ExecutionTime ' k N 2

to calculate a DFT by correlation is

proportional to the length of the DFT

DFT

 

squared.

 

where N is the number of points in the DFT and kD F T is a constant of proportionality. If the sine and cosine values are calculated within the nested loops, kDFT is equal to about 25 microseconds on a Pentium at 100 MHz. If you precalculate the sine and cosine values and store them in a look-up-table, kDFT drops to about 7 microseconds. For example, a 1024 point DFT will require about 25 seconds, or nearly 25 milliseconds per point. That's slow!

Using this same strategy we can derive the execution time for the FFT. The time required for the bit reversal is negligible. In each of the Log2 N stages there are N/2 butterfly computations. This means the execution time for the program is approximated by:

EQUATION 12-2

FFT execution time. The time required to calculate a DFT using the FFT is proportional to N multiplied by the logarithm of N.

ExecutionTime ' kFFT N log2N

The value of kFFT is about 10 microseconds on a 100 MHz Pentium system. A 1024 point FFT requires about 70 milliseconds to execute, or 70 microseconds per point. This is more than 300 times faster than the DFT calculated by correlation!

Not only is N Log2 N less than N 2 , it increases much more slowly as N becomes larger. For example, a 32 point FFT is about ten times faster than the correlation method. However, a 4096 point FFT is one-thousand times faster. For small values of N (say, 32 to 128), the FFT is important. For large values of N (1024 and above), the FFT is absolutely critical. Figure 12-8 compares the execution times of the two algorithms in a graphical form.

238

The Scientist and Engineer's Guide to Digital Signal Processing

FIGURE 12-8

Execution times for calculating the DFT. The correlation method refers to the algorithm described in Table 12-2. This method can be made faster by precalculating the sine and cosine values and storing them in a look-up table (LUT). The FFT (Table 12-3) is the fastest algorithm when the DFT is greater than 16 points long. The times shown are for a Pentium processor at 100 MHz.

 

1000

 

 

 

 

 

 

 

 

 

(seconds)

100

 

 

 

 

correlation

 

 

 

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

time

 

 

 

 

 

 

 

correlation

 

1

 

 

 

 

 

 

w/LUT

 

Execution

0.1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

FFT

 

0.01

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.001

 

 

 

 

 

 

 

 

 

 

8

16

32

64

128

256

512

1024

2048

4096

Number points in DFT

The FFT has another advantage besides raw speed. The FFT is calculated more precisely because the fewer number of calculations results in less round-off error. This can be demonstrated by taking the FFT of an arbitrary signal, and then running the frequency spectrum through an Inverse FFT. This reconstructs the original time domain signal, except for the addition of roundoff noise from the calculations. A single number characterizing this noise can be obtained by calculating the standard deviation of the difference between the two signals. For comparison, this same procedure can be repeated using a DFT calculated by correlation, and a corresponding Inverse DFT. How does the round-off noise of the FFT compare to the DFT by correlation? See for yourself in Fig. 12-9.

Further Speed Increases

There are several techniques for making the FFT even faster; however, the improvements are only about 20-40%. In one of these methods, the time

FIGURE 12-9

DFT precision. Since the FFT calculates the DFT faster than the correlation method, it also calculates it with less round-off error.

 

70

 

 

 

 

 

 

 

60

 

 

 

 

 

 

million)

50

 

 

 

 

 

 

40

 

 

 

 

 

 

per

 

 

 

 

 

 

 

 

 

 

 

 

 

Error (parts

30

 

 

 

 

 

 

20

 

correlation

 

 

 

 

 

 

 

 

 

 

10

 

 

 

 

 

 

 

 

 

 

 

FFT

 

 

 

0

 

 

 

 

 

 

 

16

32

64

128

256

512

1024

Number of points in DFT

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.

Chapter 12The Fast Fourier Transform

241

Time Domain

2

a. Real part

1

Amplitude

0

-1

0

16

32

48

643

Sample number

2

 

 

 

 

 

b. Imaginary part

 

 

 

1

 

 

 

 

Amplitude

 

 

 

 

0

 

 

 

 

-1

 

 

 

 

0

16

32

48

643

Sample number

 

Frequency Domain

 

8

 

c. Real part (odd symmetry)

Amplitude

4

0

 

-4

-8

0

16

32

48

63

Frequency

8

d. Imaginary part (even symmetry)

Amplitude

4

0

 

-4

-8

0

16

32

48

63

Frequency

FIGURE 12-11

Imaginary part symmetry of the DFT.

The next step is to modify the algorithm to calculate a single DFT faster. It's ugly, but here is how it is done. The input signal is broken in half by using an interlaced decomposition. The N/2 even points are placed into the real part of the time domain signal, while the N/2 odd points go into the imaginary part. An N/2 point FFT is then calculated, requiring about one-half the time as an N point FFT. The resulting frequency domain is then separated by the even/odd decomposition, resulting in the frequency spectra of the two interlaced time domain signals. These two frequency spectra are then combined into a single spectrum, just as in the last synthesis stage of the FFT.

To close this chapter, consider that the FFT is to Digital Signal Processing what the transistor is to electronics. It is a foundation of the technology; everyone in the field knows its characteristics and how to use it. However, only a small number of specialists really understand the details of the internal workings.

242

The Scientist and Engineer's Guide to Digital Signal Processing

3000

'FFT FOR REAL SIGNALS

 

3010

'Upon entry, N% contains the number of points in the DFT, REX[ ] contains

3020

'the real input signal, while values in IMX[ ] are ignored. Upon return,

3030

'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.

3040

'

 

3050

NH% = N%/2-1

'Separate even and odd points

3060 FOR I% = 0 TO NH%

 

3070

REX(I%) = REX(2*I%)

 

3080

IMX(I%) = REX(2*I%+1)

 

3090 NEXT I%

 

3100

'

 

3110

N% = N%/2

'Calculate N%/2 point FFT

3120

GOSUB 1000

'(GOSUB 1000 is the FFT in Table 12-3)

3130 N% = N%*2

 

3140

'

 

3150

NM1% = N%-1

'Even/odd frequency domain decomposition

3160 ND2% = N%/2

 

3170 N4% = N%/4-1

 

3180 FOR I% = 1 TO N4%

 

3190

IM% = ND2%-I%

 

3200

IP2% = I%+ND2%

 

3210

IPM% = IM%+ND2%

 

3220

REX(IP2%) = (IMX(I%) + IMX(IM%))/2

3230

REX(IPM%) = REX(IP2%)

 

3240

IMX(IP2%) = -(REX(I%) - REX(IM%))/2

3250

IMX(IPM%) = -IMX(IP2%)

 

3260

REX(I%) = (REX(I%) + REX(IM%))/2

3270

REX(IM%) = REX(I%)

 

3280

IMX(I%) = (IMX(I%) - IMX(IM%))/2

3290

IMX(IM%) = -IMX(I%)

 

3300 NEXT I%

 

3310 REX(N%*3/4) = IMX(N%/4)

 

3320 REX(ND2%) = IMX(0)

 

3330

IMX(N%*3/4) = 0

 

3340 IMX(ND2%) = 0

 

3350 IMX(N%/4) = 0

 

3360

IMX(0) = 0

 

3370

'

 

3380

PI = 3.14159265

'Complete the last FFT stage

3390 L% = CINT(LOG(N%)/LOG(2))

 

3400 LE% = CINT(2^L%)

 

3410 LE2% = LE%/2

 

3420 UR = 1

 

3430

UI = 0

 

3440 SR = COS(PI/LE2%)

 

3450

SI = -SIN(PI/LE2%)

 

3460 FOR J% = 1 TO LE2%

 

3470

JM1% = J%-1

 

3480

FOR I% = JM1% TO NM1% STEP LE%

3490

IP% = I%+LE2%

 

3500

TR = REX[IP%]*UR - IMX[IP%]*UI

 

3510

TI = REX[IP%]*UI + IMX[IP%]*UR

 

3520

REX[IP%] = REX[I%]-TR

 

3530

IMX[IP%] = IMX[I%]-TI

 

3540

REX[I%] = REX[I%]+TR

 

3550

IMX[I%] = IMX[I%]+TI

 

3560

NEXT I%

 

3570

TR = UR

 

3580

UR = TR*SR - UI*SI

 

3590

UI = TR*SI + UI*SR

 

3600 NEXT J%

TABLE 12-7

3610 RETURN

CHAPTER

Continuous Signal Processing

13

Continuous signal processing is a parallel field to DSP, and most of the techniques are nearly identical. For example, both DSP and continuous signal processing are based on linearity, decomposition, convolution and Fourier analysis. Since continuous signals cannot be directly represented in digital computers, don't expect to find computer programs in this chapter. Continuous signal processing is based on mathematics; signals are represented as equations, and systems change one equation into another. Just as the digital computer is the primary tool used in DSP, calculus is the primary tool used in continuous signal processing. These techniques have been used for centuries, long before computers were developed.

The Delta Function

Continuous signals can be decomposed into scaled and shifted delta functions, just as done with discrete signals. The difference is that the continuous delta function is much more complicated and mathematically abstract than its discrete counterpart. Instead of defining the continuous delta function by what it is, we will define it by the characteristics it has.

A thought experiment will show how this works. Imagine an electronic circuit composed of linear components, such as resistors, capacitors and inductors. Connected to the input is a signal generator that produces various shapes of short pulses. The output of the circuit is connected to an oscilloscope, displaying the waveform produced by the circuit in response to each input pulse. The question we want to answer is: how is the shape of the output pulse related to the characteristics of the input pulse? To simplify the investigation, we will only use input pulses that are much shorter than the output. For instance, if the system responds in milliseconds, we might use input pulses only a few microseconds in length.

After taking many measurement, we come to three conclusions: First, the shape of the input pulse does not affect the shape of the output signal. This

243

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

a.

b.

c

Chapter 13Continuous Signal Processing

245

Linear

System

Linear

System

Linear

System

 

*(t)

d.

Linear

System

 

FIGURE 13-1

The continuous delta function. If the input to a linear system is brief compared to the resulting output, the shape of the output depends only on the characteristics of the system, and not the shape of the input. Such short input signals are called impulses. Figures a,b & c illustrate example input signals that are impulses for this particular system. The term delta function is used to describe a normalized impulse, i.e., one that occurs at t ' 0 and has an area of one. The mathematical symbols for the delta function are shown in (d), a vertical arrow and *(t).

obvious, since the viewing geometry is drastically different. Mars is about 6000 kilometers in diameter and 60 million kilometers from earth. In comparison, Sirius is about 300 times larger and over one-million times farther away. These dimensions should make Mars appear more than three-thousand times larger than Sirius. How is it possible that they look alike?

These objects look the same because they are small enough to be impulses to the human visual system. The perceived shape is the impulse response of the eye, not the actual image of the star or planet. This becomes obvious when the two objects are viewed through a small telescope; Mars appears as a dim disk, while Sirius still appears as a bright impulse. This is also the reason that stars twinkle while planets do not. The image of a star is small enough that it can be briefly blocked by particles or turbulence in the atmosphere, whereas the larger image of the planet is much less affected.

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

Соседние файлы в предмете Электроника