Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

The Fast Fourier Transform Historians of mathematics have found that Gauss was aware of a method for computing DFT sums in O(N log N ) steps. It was rediscovered and published as “An algorithm for machine calculation of complex Fourier series” by J.W. Cooley and J.W. Tuckey Mathematics of Computation, 19, 297 (1965), and it rapidly became very popular. It was selected in 1999 as one of the Top 10 Algorithms of the 20th Century. The discrete Fourier transform and its inverse are defined by N −1 1 X kn gk = √ W fn , N n=0 N −1 1 X −nk W gk , fn = √ N k=0

where W is the N -th root of unity W ≡ ei2π/N . The problem with these sums is that they involve O(N 2 ) operations, because there are O(N ) sums, and each sum involves O(N ) operations. A discrete Fourier transform of a time series with N points can be viewed as multiplication of the N -component column vector of values fn by an N × N matrix with components W kn to yield the N values of gk . The key to improving this quadratic dependence on N is to note that there are only N different values W ≡ ei2πn/N ,

n = 0, 1, . . . , N − 1 ,

i.e., the N -th roots of unity, among the N 2 matrix elements of the DFT considered as a matrix multiplication. In addition, there are only N different time series values. The problem is to find a way to factor these sums to reduce the number multiplications and additions. PHY 410-505 Computational Physics I

1

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

The Danielson-Lanczos Lemma The FFT can be understood using the Danielson-Lanczos lemma, which states that a discrete Fourier transform of N points can be decomposed into two transforms of N/2 points each as follows: N/2−1

gk =

X

N/2−1

W

k(2n)

f2n +

n=0

X n=0

N/2−1

= ≡

X n=0 gkeven

W k(2n+1) f2n+1

N/2−1 2 kn

(W ) f2n + W

k

X

(W 2 )kn f2n+1

n=0

+

W k gkodd

Note that W 2 = exp(i2π/(N/2)) is the appropriate phase factor for a series with N/2 points. If N = 2m with m an integer, this decomposition into even and odd Fourier transforms can be carried out recursively until one is left with Fourier transforms of series with 1 point, which are trivial to perform. There are m = log2 N steps in such a recursion. At each step in the recursion, the split into even and odd transforms requires one complex multiplication and addition to compute each of the N quantities gk . Thus the total number of complex operations required by the algorithm is on the order of N log2 N . Bit-Reversal Magic It is somewhat tricky to keep track of the decomposition into even and odd series in each step of the recursion. An ingenious algorithm for doing this was discovered by Cooley and Tuckey. This algorithm involves a “bit-reversal” reordering of the set of data points {fi }, followed by succesive applications of the Danielson-Lanczos formula with transforms of lengths 2, 4, . . . , N . The significance of the bit-reversal reordering is easy to understand by considering some small values of N . When N = 1, the Fourier transform is trivial: g0 = f0 . PHY 410-505 Computational Physics I

2

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

When N = 2, W = exp(i2π/2) = −1. The Fourier transform is given by 

g0 g1



 =

f0 + f1 f0 − f1



From this formula we see that g even ≡ g 0 =

 =



f0 f0

f0 f0



 +



(−1)0 0

0 (−1)1

g odd ≡ g 1 =

,



f1 f1



f1 f1

 .

 .

In this formula, we have abbreviated “even” by the “bit” 0, and “odd” by the “bit” 1; we see that the “bit patterns” 0 and 1 are identical with the corresponding subscripts of f , namely 0 and 1 in the even and odd transforms. Note also that the even and odd transforms are periodic with period N/2 = 1. even odd In the general case, it is easy to show that gk+(N/2) = gkeven and gk+(N/2) = gkodd . This periodicity can be exploited, using W k+(N/2) = −W k , to reduce the number of complex multiplications from N to N/2 in the even-odd decomposition:

gk = gkeven + W k gkodd ,

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

gk+(N/2) = gkeven − W k gkodd ,

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

The need for bit-reversal becomes apparent when N = 4. In this case, W = exp(i2π/4) = i. The Fourier transform is thus     g0 f0 + f1 + f2 + f3  g1   f0 + if1 − f2 − if3   =  , g2 f0 − f1 + f2 − f3 g3 f0 − if1 − f2 + if3 whence,  f0 + f2  f − f2  ≡ g0 =  0  , f0 + f2 f0 − f2 

g even

PHY 410-505 Computational Physics I

 f1 + f3  f − f3  ≡ g1 =  1  . f1 + f3 f1 − f3 

g odd

3

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

These even and odd transforms are periodic with period 2 in the index of g. Each of these periodic parts has precisely the same structure as an N = 2 transform. Let us therefore decompose gk0 = gk00 + (−1)k gk01 and gk1 = gk10 + (−1)k gk11 , where         f0 f2 f1 f3 f  f  f  f  g 00 =  0  , g 01 =  2  , g 10 =  1  , g 11 =  3  . f0 f2 f1 f3 f0 f2 f1 f3 If we express the subscripts of f in this formula in binary notation, we see that there is the following correspondence between the bit patterns labeling the g’s and the corresponding subscripts: Even − odd pattern n in fn (decimal) n (binaryvalue) 00 0 00 01 2 10 10 1 01 11 3 11 Each bit pattern in the even-odd decomposition is identical to the bit-reversed binary form of the corresponding subscript n! The strategy behind the Cooley-Tuckey algorithm is to organize the recursive Danielson-Lanczos decomposition steps so that the final 1-dimensional transforms (in the first column above) are in increasing binary order. If the N data points fn have been reordered by bit-reversal, then the arrays representing the data points and the 1-dimensional transforms will be identical! The N -dimensional transform is then built up recursively with transforms of length 2, 4, . . . , N . This procedure can be illustrated for the case N = 4: {f0 , f1 , f2 , f3 } ↓ PHY 410-505 Computational Physics I

4

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

{f0 , f2 , f1 , f3 } | {z } | {z } N =2

N =2

↓ {f0 + (−1)0 f2 , f0 − (−1)0 f2 , f1 + (−1)0 f3 , f1 − (−1)0 f3 } | {z } N =4

↓ {(f0 + f2 ) + i0 (f1 + f3 ), (f0 − f2 ) + i1 (f1 − f3 ), (f0 + f2 ) − i0 (f1 + f3 ), (f0 − f2 ) − i1 (f1 − f3 ), } k {f0 + f1 + f2 + f3 , f0 + if1 − f2 − if3 , f0 − f1 + f2 − f3 , f0 − if1 − f2 + if3 } The first bit-reversal step effectively performs 4 N = 1 transforms. Two N = 2 transforms, with W = −1 are performed in the second step using the Danielson-Lanczos lemma, and an N = 4 transform with W = i in the final step completes the process. Note that the transformation has been done “in place”: the input data array fn is transformed into the output gk array, and no extra storage arrays are needed. Power Spectrum or Periodogram The spectrum of frequencies in the DFT is ωk =

2k ωc , N

k=−

N N N N , − + 1, . . . , −1, 0, 1, . . . , − 1, , 2 2 2 2

where ωc = π/δt is called the Nyquist critical frequency. The one-side discrete power spectrum, or periodogram, is defined as follows:  2 |g  h 0| i k=0  1 2 2 |gk | + |gN −k | k = 1, 2, . . . , N2 −1 P (ωk ) = N  gN/2 2 k = N2 PHY 410-505 Computational Physics I

5

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

Power Spectrum of Nonlinear Pendulum The following program computes the power spectrum of the nonlinear driven and damped pendulum from Chapter 3 using the FFT. pendulumPower.cpp #include #include #include #include



#include "Vector.hpp"

1 2 3 4 // vectors and FFT algorithm

const double pi = 4 * std::atan(1.0); double q = 0.5; double omega_0 = 2.0/3.0; double b = 0.9;

8

// damping coefficient // frequency of driving force // amplitude of driving force

cpl::Vector f(cpl::Vector y) { double t = y[0]; double theta = y[1]; double omega = y[2]; cpl::Vector f(3); // Vector with 3 components f[0] = 1; f[1] = omega; f[2] = - q * omega - std::sin(theta) + b * std::cos(omega_0 * t); return f; }

PHY 410-505 Computational Physics I

6

6

10 11 12 14 15 16 17 18 19 20 21 22 23

Friday, October 29, 2004

Chapter 5

Spectral Analysis and Gaussian Quadrature

Lecture 2

void RK4Step(cpl::Vector& y, double h) { cpl::Vector k1 = h * f(y); cpl::Vector k2 = h * f(y + 0.5 * k1); cpl::Vector k3 = h * f(y + 0.5 * k2); cpl::Vector k4 = h * f(y + k3); y += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0; }

25 26 27 28 29 30 31

int main() { std::cout