Xem mẫu
Digital Signal Processing: Laboratory Experiments Using C and the TMS320C31 DSK Rulph Chassaing
Copyright © 1999 John Wiley & Sons, Inc. 6 Print ISBN 0-471-29362-8 Electronic ISBN 0-471-20065-4 Fast Fourier Transform
O The fast Fourier transform using radix-2 and radix-4
O Decimation or decomposition in frequency and in time O Programming examples
The fast Fourier transform (FFT) is an efficient algorithm that is used for con-verting a time-domain signal into an equivalent frequency-domain signal, based on the discrete Fourier transform (DFT). A real-time programming example is included with a main C program that calls an FFT assembly function.
6.1 INTRODUCTION
The discrete Fourier transform converts a time-domain sequence into an equiva-lent frequency-domain sequence. The inverse discrete Fourier transform per-forms the reverse operation and converts a frequency-domain sequence into an equivalent time-domain sequence. The fast Fourier transform (FFT) is a very ef-ficient algorithm technique based on the discrete Fourier transform, but with fewer computations required. The FFT is one of the most commonly used oper-ations in digital signal processing to provide a frequency spectrum analysis [1–6]. Two different procedures are introduced to compute an FFT: the decima-tion-in-frequency and the decimation-in-time. Several variants of the FFT have been used, such as the Winograd transform [7,8], the discrete cosine transform (DCT) [9], and the discrete Hartley transform [10–12]. Programs based on the DCT, FHT, and the FFT are available in [9].
6.2 DEVELOPMENT OF THE FFT ALGORITHM WITH RADIX-2
The FFT reduces considerably the computational requirements of the discrete Fourier transform (DFT). The DFT of a discrete-time signal x(nT) is
165
166 Fast Fourier Transform
N – 1
X(k) = ^ x(n) Wnk k = 0, 1, . . . , N – 1 (6.1) n = 0
where the sampling period T is implied in x(n) and N is the frame length. The constants W are referred to as twiddle constants or factors, which represent the phase, or
W = e–j2p/N (6.2)
and is a function of the length N. Equation (6.1) can be written for k = 0, 1, . . . , N – 1, as
X(k) = x(0) + x(1)Wk + x(2)W2k + . . . + x(N – 1)W(N–1)k (6.3)
This represents a matrix of N × N terms, since X(k) needs to be calculated for N values of k. Since (6.3) is an equation in terms of a complex exponential, for each specific k there are approximately N complex additions and N complex multiplications. Hence, the computational requirements of the DFT can be very intensive, especially for large values of N.
The FFT algorithm takes advantage of the periodicity and symmetry of the twiddle constants to reduce the computational requirements of the FFT. From the periodicity of W
Wk+N = Wk (6.4)
and, from the symmetry of W
Wk+N/2 = –Wk (6.5)
Figure 6.1 illustrates the properties of the twiddle constants W for N = 8. For ex-ample, let k = 2, and note that from (6.4), W10 = W2, and from (6.5), W6 = –W2.
FIGURE 6.1 Periodicity and symmetry of twiddle constant W.
6.3 Decimation-in-Frequency FFT Algorithm with Radix-2 167
For a radix-2 (base 2), the FFT decomposes an N-point DFT into two (N/2)-point or smaller DFT’s. Each (N/2)-point DFT is further decomposed into two (N/4)-point DFT’s, and so on. The last decomposition consists of (N/2) two-point DFT’s. The smallest transform is determined by the radix of the FFT. For a radix-2 FFT, N must be a power or base of two, and the smallest transform or the last decomposition is the two-point DFT. For a radix-4, the last decomposi-tion is a four-point DFT.
6.3 DECIMATION-IN-FREQUENCY FFT ALGORITHM WITH RADIX-2
Let a time-domain input sequence x(n) be separated into two halves:
a) x(0), x(1), . . . , x1N – 12 (6.6) and
b) 1}2, x1} + 12, . . . , x(N – 1) (6.7) Taking the DFT of each set of the sequence in (6.6) and (6.7),
(N/2) – 1 N – 1
X(k) = ^ x(n)Wnk + ^ x(n)Wnk (6.8) n = 0 n = N/2
Let n = n + N/2 in the second summation of (6.8), X(k) becomes
(N/2) – 1 (N/2) – 1
X(k) = ^ x(n)Wnk + WkN/2 ^ x n + } Wnk (6.9) n = 0 n = 0
where WkN/2 is taken out of the second summation because it is not a function of n. Using,
WkN/2 = e–jkp = (e–jp)k = (cos p – jsin p)k = (–1)k
in (6.9), X(k) becomes
(N/2) – 1
X(k) = ^ x(n) + (–1)kx n + } Wnk (6.10) n = 0
Because (–1)k = 1 for even k and –1 for odd k, (6.10) can be separated for even and odd k, or
168 Fast Fourier Transform
for even k:
for odd k:
(N/2) – 1
X(k) = ^ x(n) + x n + } Wnk (6.11) n = 0
(N/2) – 1
X(k) = ^ x(n) – x n + } Wnk (6.12) n = 0
Substituting k = 2k for even k, and k = 2k + 1 for odd k, (6.11) and (6.12) can be written as, for k = 0, 1, . . . , (N/2) – 1,
(N/2) – 1
X(2k) = ^ x(n) + x n + } W2nk (6.13) n = 0
(N/2) – 1
x(2K + 1) = ^ x(n) – x n + } WnW2nk (6.14) n = 0
Because the twiddle constant W is a function of the length N, it can be repre-sented as WN. Then, WN can be written as WN/2. Let
a(n) = x(n) + x(n + N/2) (6.15)
b(n) = x(n) – x(n + N/2) (6.16)
Equations (6.13) and (6.14) can be more clearly written as two (N/2)-point DFT’s, or
(N/2) – 1
X(2k) = ^ a(n)WN/2 (6.17) n = 0
(N/2) – 1
X(2k + 1) = ^ b(n)WNWN/2 (6.18) n = 0
Figure 6.2 shows the decomposition of an N-point DFT into two (N/2)-point DFT’s, for N = 8. As a result of the decomposition process, the X’s in Figure 6.2 are even in the upper half and they are odd in the lower half. The decomposition process can now be repeated such that each of the (N/2)-point DFT’s is further decomposed into two (N/4)-point DFT’s, as shown in Figure 6.3, again using N = 8 to illustrate.
The upper section of the output sequence in Figure 6.2 yields the sequence X(0) and X(4) in Figure 6.3, ordered as even. X(2) and X(6) from Figure 6.3 rep-resent the odd values. Similarly, the lower section of the output sequence in Fig-ure 6.2 yields X(1) and X(5), ordered as the even values, and X(3) and X(7) as the odd values. This scrambling is due to the decomposition process. The final
6.3 Decimation-in-Frequency FFT Algorithm with Radix-2 169
FIGURE 6.2 Decomposition of N-point DFT into two (N/2)-point DFT’s, for N = 8.
order of the output sequence X(0), X(4), . . . in Figure 6.3 is shown to be scram-bled. The output needs to be resequenced or reordered. A special instruction us-ing indirect addressing with bit-reversal, introduced in Chapter 2 in conjunction with circular buffering, is available on the TMS320C3x to reorder such a se-quence. The output sequence X(k) represents the DFT of the time sequence x(n). This is the last decomposition, since we have now a set of (N/2) two-point DFT’s, the lowest decomposition for a radix-2. For the two-point DFT, X(k) in
(6.1) can be written as
FIGURE 6.3 Decomposition of two (N/2)-point DFT’s into four (N/4)-point DFT’s, for N = 8.
...
- tailieumienphi.vn
nguon tai.lieu . vn