D Details on FFT codes provided by the Numerical recipes book[2]

The following is about a specific FFT subroutine provided by the Numerical recipes book[2]. This is not a general case.

The input and output of the DFT are usually complex numbers. In the implementation of FFT algorithm provided by Numerical recipes book[2], to avoid using complex numbers, the algorithm adopts the real number representation of the complex numbers. In this scheme, two elements of a real number array are used to store one complex number. To store a complex array “cdata” of length N, we will need a real number array “rdata” of length 2N. The first elements of array “rdata” will contain the real part of “cdata(1)”, the second elements of “rdata” will contain the imaginary part of  “cdata(1)”, and so on.

To test the correctness of the above statement, I generated a real number array with length N = 2 × 28 by using a random generating routine and calculate the DFT of the array with two methods. The real array generated by the random generator are considered to be a real number representation of a complex array with length N∕2. Using the real array as the input of the FFT routine (the code in ˜/project_new/fft). To check the correctness of my understanding of the input and output of the FFT, I manually convert the real number of length N to a complex array with length N∕2, and use directly the summation in Eq. (??) to calculate the DFT. The output I got is obviously a complex array with length N∕2. Then I manually convert the complex array to a real number array of length N and plot the output in Fig. 6 with dashed line. The results in Fig. 6 indicates the results given by the FFT and the naive method used by me agree with each other well. This proves that my understanding of the input and output of FFT (especially the storage arrangement) is correct.


pict

Figure 6: The output of a FFT routine (solid line) agrees well with those (dashed line) calculated by directly evaluate the summation in Eq. (??) (the latter is coded by me). The agreement indicates my understanding of the output of FFT (especially the storage arrangement) is correct. Here the time domain data is generated by a random generating routine.

To clearly show the output of FFT, we recover the real and imaginary part of DFT from the output of FFT and plots the data as a function of their corresponding frequency. The results are given in Fig. 7.


pict pict

Figure 7: The real (left) and imaginary (right) part of the discrete Fourier transformation as a function of the frequency. The point n = 128 corresponds to frequency fs2 while the point n = 128 corresponds to frequency +fs2, where fs = 1Δ with Δ is the sample interval in time domain. The time domain data is the same as used in Fig. 6. Why is there a spike at zero frequency?

 

 

 D.1 Computing Fourier integrals using FFT (not finished, to be deleted)