3

.

THE FAST FOURIER TRANSFORM

3–1

Chapter 3: THE FAST FOURIER TRANSFORM

3–2

§3.1 FFT IN A NUTSHELL

Even though this course is not about any numerical recipe, in order for us to utilize the numerical tools in intelligent ways, it pays to understand the essnetial features of the tools we would frequently use. the Fast Fourier Transform (FFT) is one such numerical tool. According to G. Strang, FFT is perhaps the most widely used numerical tool today in scientific and engineering computations. To c version the Fourier transform: start our discussion of FFT, let us recall MATLAB

X (k) =

N −1 

x(s)e− j2πkn/N

, k = 0, 1, . . . , N − 1

(3.1)

n=0

where the denominator 1/N is omitted for brevity. The one thing we have taken for granted is that for every kthe Fourier coefficient term, X (k), one must perform N-complex term multiplications and N additions, thus requiring N 2 multiplications in order to obtain the complete coefficients. Hence, for a measurement data set of 100 with 1,024 sample points, one would need approximately a total of 108 complex-function multiplications just to obtain the Fourier coefficients, not to mention the inverse Fourier transform. Now, if the fundamental frequency is about 1 Hz and one wishes to capture up to 1 KHz frequency components, then the minimum number of samples would be about 4,096. This then would require about 1.6 × 109 operations. An FFT algorithm would reduce the multiplications by a factor of 6N : N2 224 ≈ 682! = 6N 6 × 212

(3.2)

To appreciate the tremendous computational efficiency improvements possible by FFT, let us consider the case of N = 4 (4 samples per period, a very crude sampling!) and introduce the following notation w = e− j2π/N , or

w = e− jπ/2 for N = 4

(3.3)

Hence, the four equations from (3.1) can be written as n=0: n=1: n=2: n=3:

 X (0) = x((0)w0 + x((1)w0 + x((2)w0 + x((3)w0    0 1 2 3 X (1) = x((0)w + x((1)w + x((2)w + x((3)w  X (2) = x((0)w0 + x((1)w2 + x((2)w4 + x((3)w6     0 3 6 9 X (3) = x((0)w + x((1)w + x((2)w + x((3)w

which can be written in matrix form as 3–2

(3.4)

3–3

§3.1

FFT IN A NUTSHELL

{X4 } = F4 {x4 } X4 = [ X (0) x4 = [ x(0)  0 w  w0 F4 =  0 w w0

X (1)

X (2)

X (3) ]T

x(1) x(2) x(3) ]T  w0 w0 w0 w1 w2 w3   w2 w4 w6 w3 w6 w9

(3.5)

Note from (3.3) we have w6 = w2 w0 = 1



e− jπ/2

since and

6

2   = e− jπ = e− jπ/2

w4 = w0

so that F4 is simplified to: 

1 1  1 w1 F4 =  1 w2 1 w3

1 w2 1 w2

 1 w3   w2 w1

(3.6)

Observe from the above equation that the computation of X (k) will require N 2 complex multiplications and N (N − 1) complex additions. Hence, for a large N , computations can become quickly prohibitive. (Remember FFT is often used in conjunction with real-time identification problems!) The key in reducing the computational overhead drastically is the following product form of F4 : 

1 0 0 1 F4 =  1 0 0 1

1 0 −1 0

 0 1 w0 w1   1 w2  0 0 0 1 0 0 −w

0 0 1 w0

 0 1 0 0  w0 0 2 w 0

0 0 1 0

0 1 0 0

 0 0  0 1

(3.7)

Notice that the third matrix in the above equation is just a permutation operator; the even-numbered samples are stacked first in the vector x4 , which is followed by the odd-numbered vectors. Specifically, we have 

  x(0) 1  x(2)   0  = x(1) 0 x(3) 0

0 0 1 0

0 1 0 0

Thus, the following pattern emerges from this equation 3–3

  0 x(0) 0   x(1)    0 x(2) 1 x(3)

3–4

Chapter 3: THE FAST FOURIER TRANSFORM



I F4 = 2 I2

W2 −W2



1 W2 = 0

0 w1



0 F2

F2 0



even/odd permutation operator



 (3.8)

w0 w2

1 F2 = w0





It is the generality of the preceding product-form decomposition that renders FFT an efficient tool. For example, for N = 512, we have the following product form:



F512

I = 256 I256

W256 −W256



0 F256

F256 0



even/odd permutation operator

 (3.9)

W256 = diag(w0 , w1 , . . . , w255 )

The successive inversions of the first product matrix is trivial since W is diagonal:



I256 I256

W256 −W256

−1

 =

1 2

I256 W∗256

I256 −W∗256



where W∗ is the complex conjugate of W.

The key aspect of the FFT algorithm is in the successive decomposition of F512 into: 3–4

(3.10)

3–5

§3.1 

F256

I = 128 I128

W128 −W128



F128 0

0 F128



FFT IN A NUTSHELL

even/odd permutation operator



W128 = diag(w0 , w1 , . . . , w128 ) 

F128

I = 64 I64

W64 −W64



0 F64

F64 0



even/odd permutation operator



W64 = diag(w0 , w1 , . . . , w64 ) 

F64

I = 32 I32

W32 −W32



F32 0

0 F32



even/odd permutation operator



(3.11)

W64 = diag(w0 , w1 , . . . , w32 ) ... ...     I2 W2 F2 0 even/odd permutation F4 = I2 −W2 operator 0 F2 W = diag(w0 , w1 )

The result of the FFT algorithm outlined above reduces the discrete Fourier transform from O(N 2 ) to O( 12 N log2 N ) algorithm, a marvelous achievement by numerical analysts during the 1960s. It turns out that there is a competitive transform that works exclusively from time-domain signals to time-domain expansion coefficients. This is called wavelet transform, which we will study later √ in the course. Observe that each column in F N , if normalized by N , constitutes an orthonormal vector. However, each basis function cos 2πkt in the Fourier series expansion is non-compact. The wavelet basis functions are, on the other hand, compact. Thus, a race has begun between FFT and FWT (Fast Wavelet Transform). References for FFT 1.

G-AE Subcommittee on Measurement Concepts. “What is the Fast Fourier Transform?” IEEE Trans. Audio and Electroacustics. (June 1967), Vol. AU-15, pp. 45-55. Also Proc. IEEE (Oct. 1967), Vol. 55, pp. 1664-1674.

2.

Bergland, G.D. “A Guided Tour of the Fast Fourier Transform.” IEEE Spectrum (1969), Vol. 6, No.7, pp. 41-52.

3.

Brigham, E.O., The Faster Fourier Transform and Its Applications, Prentice-Hall, 1988.

4.

Brigham, E. O., and R. E. Morrow. “The Fast Fourier Transform.” IEEE Spectrum (December 1967), Vol. 4, pp. 63-70. 3–5

Chapter 3: THE FAST FOURIER TRANSFORM

3–6

5.

Cooley, J. W., and J. W. Tukey. “An Algorithm for Machine Calculation of Complex Fourier Series.” Math Computation (April 1965), Vol. 19, pp. 297-301.

6.

Cooley, J. W., P. A. W. Lewis, and P. D. Welch. “The Finite Fourier Transform.”IEEE Trans. on Audio and Electroacoustics (June 1969), Vol. AU17, No. 2, pp. 77-85.

7.

Cooley, J. W., P. A. W. Lewis, and P. D. Welch. “The Fast Fourier Transform and its Applications.” IEEE Trans on Education (March 1969), Vol. 12, pp. 27-34.

8.

Cooley, J. W., P. A. W. Lewis, and P. D. Welch. “The Fast Fourier Transform Algorithm: Programming Consideration in the Calculation of Sine, Cosine and Laplace Transforms.” J. Sound Vib. (July 1970), Vol. 12, No. 3, pp. 315-337.

9.

Durvasula, S., B. Panda, and R. K. Kapania, “Impulse Excitation Technique For Determining Vibration Characteristics Of Aeroelastic Models of Tall Structures,”Inst Eng India Part ND 2, Vol. 61, May 1981.

10. Eshleman, R. L., “Machinery Diagnostics And Your FFT,” S. V. Sound and Vibration, Vol. 17, No. 4, pp. 12-18, April 1983. 11. Gentleman, W.M., and G. Sande. “Fast Fourier Transform for Fun and Profit,” AFIPS Proc., 1966 Fall Joint Computer Conf., Vol. 29, pp. 563-678, Washington, DC: Spartan, 1966. 12. Halliwell, G., “Using Fourier Transform Analysers To Investigate Component Vibration,” Nuclear Engineering International (GB), Vol. 25, 296, pp. 42-45, March 1980. 13. Harris, F. J. “On the Use of Windows for Harmonics Analysis with the Discrete Fourier Transform.” Proc. IEEE (January 1978), Vol. 66, No. 1, PP. 51-83. 14. Karlsson, L., “Numerical Analysis of Damped Transient Beam Vibration by Use of Fourier Transforms,” Int. J. Numer. Methods Eng. (GB), Vol. 21, No. 4, pp. 683-89, April 1985. 15. Parmenter, W. W., and R. G. Christiansen, “Recovery of Modal information Form a Beam Undergoing Random Vibration,” Trans. ASME Ser. B., Vol. 96, No. 4, pp. 1307-13, November 1974. 16. Radcliffe, C. J., and C. D. Mote, Jr., “Identification and Control of Rotating Disk Vibration,” Trans ASME J. Dyn. Syst. Meas. and Control, Vol. 105, No. 1, pp. 39-45, March 1983. 17. Ramirez, R. W. The FFT Fundamentals and Concepts. Englewood Cliffs, NJ: Prentice-Hall, 1985. 18. Taylor, T. D., R. S. Hirsh, and M. M. Nadworny, “Comparison of FFT, Direct Inversion and Conjugate Gradient Methods For Use in Pseudo-Spectral Methods,” Computers & Fluids, Vol. 12, No. 1, pp. 1-9, 1984. 3–6

3–7

§3.1

FFT IN A NUTSHELL

19. Webb, R., “Frequency Somain Instrumentation Techiniques For The Condition Monitoring of Rotating Machinery,” Noise and Vib. Control Worldwide (GB), Vol. 14, No. 8, pp. 215-19, October 1983.

3–7