IDL and FORTRAN FFT Comparison

Index

Original Spectrum Data | FFT Input Data | FFT Output Data | Data Files and Source Code

Seeing as we here at CIMSS/SSEC do lots and lots of FFT's for spectrum -> interferogram transformation and vice versa, I put out the following data so that we can be sure that the particular FFT algorithm (in whatever language) does not introduce artifacts into our work.

The software that gets shipped with the AERI systems uses the Numerical Recipes FFT code REALFT and FOUR1. A description of how these algorithms work can be found on pp498-508 in Press et al, "Numerical Recipes in FORTRAN. The Art of Scientific Computing", 2nd Ed, Cambridge University Press 1992. The IDL FFT employs a prime factor algorithm that returns the complex discrete Fourier transform of the input.

While the Fortran and IDL FFT's produce essentially the same result, the structure of the input and output data arrays, while consistent within algorithms, is not always clear. The purpose of this web page is to explain how to input data to these routines and provide an example.

Original Spectrum Data

The spectral data used was an average of the AERI-TWP1 hot blackbody (HBB) uncalibrated spectrum, shown in figure 1, derived from the forward scan data. The data was collected from 200-2000cm^-1 and placed in an array that extended from 0cm^-1 to the Nyquist frequency, in this case 7899.5cm^-1 for a laser wavenumber of 15799cm^-1.

The number of points in the real and imaginary spectra is 16385. The first point corresponds to 0cm^-1 and the last point, 16385, corresponds to the Nyquist frequency. The total number of complex valued points is thus 2¹⁴ + 1.

Figure 1.

To referring text


Figure 1.

FFT Input Data

The input to the IDL and Fortran FOUR1 FFT algorithms is required to include both the positive and negative frequencies of the real and imaginary components. figure 2 shows how the data was reflected about the Nyquist frequency to give a power of 2 number of points, 2¹⁵.

Figure 2.

To referring text


Figure 2.

To provide a real, asymmetric interferogram upon Fourier transformation, the input spectrum must be Hermitian. Thus the imaginary component is multiplied by -1 after being reflected. This gives the FFT input spectrum as shown in figure 3.

Figure 3.

To referring text


Figure 3.

The data as shown in figure 3 is what is used directly in the IDL FFT function. For the Fortran FOUR1 the real and imaginary parts of the reflected data are provided in a REAL*4 array with the real and imaginary components alternating as shown in figure 4.

Figure 4.

To referring text


Figure 4.

This format of spectrum input is also used by REALFT except that only the positive spectral frequencies are input and the Nyquist frequency is not included, i.e. only 2¹⁴ points are used with the real and imaginary components of each point alternating in the REALFT input array. It is not altogether clear as to what happens to the Nyquist point. In this example the Nyquist point is zero (due to the spectral extension) but this may not always be the case.

FFT Output Data

The output of the FFT's are shown in figure 5. As is evident they do not agree. However, a number of things can be done. Firstly, both the IDL and FOUR1 FFT results should be divided by 2 as they use twice the number of input points and thus the 1/N scaling that occurs in the inverse FFT will be different from REALFT by a factor of 2. The result of scaling the interferograms is shown in figure 6. Scaling the FOUR1 result makes it agree with the REALFT result. In addition, the scaled IDL FFT result magnitudes is now comparable to the other two interferograms. Closer inspection of figure 6 shows that the FOUR1/REALFT and IDL interferograms appear to be reversed with respect to each other. By flipping the FOUR1/REALFT interferograms about ZPD, we achieve agreement between all three results, as shown in figure 7! Yay. There thus appears to be some ambiguity between the FOUR1/REALFT and IDL FFT as to what are exactly the negative and positive parts of the inverse transform. This is not a problem as the magnitudes are exactly the same - you just have to be careful about how you interpret the results, which is true any time you use Fourier transforms.

Figure 5.

To referring text


Figure 5.

Figure 6.

To referring text


Figure 6.

Figure 7.

To referring text


Figure 7.

Data Files and Source Code

The following links provide the original input data file and the Fortran and IDL code I used to perform the above analysis.

Input data (binary), fft_input.dat
IDL FFT test program, test_idlfft.pro
Fortran FFT test code,
Ancillary IDL code,
- Function to read input data, read_fft_input.pro
- Function to read REALFT output, read_realft.pro
- Function to read FOUR1 output, read_four1.pro
- Procedure to create the data plots, plot_fft_test.pro

[ TOP OF PAGE | UP ]

This page maintained by Paul van Delst
Last updated January 22, 1997