FUN WITH SPECTRA Computing Fourier transforms using the Fast Fourier Transform (FFT) algorithm is a standard tool in geophysical analysis. There are many books and articles on Fourier theory and how the FFT works. I am distributing the part of Numerical Recipes that discusses this, which seems as good a treatment as any. You will certainly learn much more about this if you take our Data Analysis class. For brevity in this class I am simply going to focus on what you need to know to use an FFT program and leave the theory for you to pick up elsewhere. We start with a "time series," which is a series of measurements that are made at some regular interval. These measurments are normally stored on the computer in an array. We also need to know the sample interval (the time between sample points, the reciprocal of the digitization rate) and the number of points in the time series. In most cases the time series consists of real (i.e., not complex) numbers and we will assume this is true for the examples shown. The purpose of the FFT is to examine the different frequency waves that may be present in the data. The FFT is a linear transformation that allows us to look at the data in the "frequency domain" rather than the "time domain" of the original data. In fact, the original time series can be expressed exactly as a sum of harmonic waves of different frequency content and phase. The FFT allows us to readily compute what the frequencies and phases of these waves are. The simplest FFT algorithms require that the number of points in the original time series be a power of two. We will assume that this is true here, but you should be aware that alternative FFT methods exist for non-powers of two. Thus, the input to a typical FFT algorithm will consist of an array of data values, the number of time points (n), and the sample interval (dt). For a real time series, the FFT algorithm will then return a complex spectra that will give the amplitude and phase at a series of equally spaced frequency points. There will be n/2 + 1 frequency points and the frequency spacing will be df = 1/(n*dt). The first frequency point is at f = 0.0 The last frequency point is at f = 1/(n*dt) and is called the Nyquist frequency. The Nyquist frequency is the highest frequency that one can resolve in the time series. At the Nyquist frequency, there are two data points in the time series per wavelength. If higher frequencies are present in the data, they will be "aliased" so that they are seen at lower frequencies. From the time series alone, there is no way to discriminate between correctly resolved frequencies and aliased data from higher frequencies. Thus, sometimes filters (anti-aliasing filters) are applied to the data prior to digitization to remove the frequencies above the Nyquist so that this will not be a problem. The FFT algorithm provides the spectra at n/2 + 1 frequency points which may be represented as n/2 + 1 complex numbers. (complex numbers are necessary to represent both the amplitude and phase of each harmonic component). "Wait a minute," the alert reader might cry. We started with n real points in the time domain and now we have n+2 real points in the frequency domain, if we count each complex number as two reals. Have we gained information compared to the original time series? The resolution to this apparent problem comes from the fact that the first and last frequency points always have zero imaginary parts. Thus the total number of independent data points is the same in both the time and frequency domains. In fact, most FFT algorithms return only n/2 frequency points by packing the real part of the Nyquist point into the imaginary part of the zero frequency point. Once we have the frequency domain points we can make an amplitude spectra plot by plotting the absolute values (Cabs function) of the points versus frequency. If we want a power spectrum, then we square the amplitudes. Getting the correct units on the y-axis can be a tricky business (let's save that discussion until your time series analysis class); in many cases, however, we are only interested in the position of the peaks or the shape of the spectrum so we don't need to worry about this. We can also transform back to the time domain by computing an inverse FFT. If your code is working correctly, you should get back exactly the same time series as you started with. Here is a set of routines, adopted from those in Numerical Recipes to do both forward and inverse FFTs. You may obtain these in the /net/fishpub/239 directory under the name fftsubs.c if you don't want to cut and paste from the class web site (use anonymous ftp to mahi.ucsd.edu if you don't have a fishnet account). ------------------------------------------------------------------------ /* getspec computes the spectrum of a real time series using FFT routines that are compatitible with the codes in Numerical Recipes. Requires: taper, realfft, fft MAX_TS_PTS defined as maximum number of time points Inputs: ts = input time series with values from ts[1] to ts[ntime] ntime = number of time points (int) dt = sample interval (s) (float) itap = 1 to apply Hann taper before computing sprectrum = 0 otherwise Returns: spec = complex spectrum stored as floats (see below) nfreq = number of complex frequency points (int, = 1+ntime/2) df = frequency spacing spec is packed as follows: spec[0] = not used spec[1], spec[2] = real & imag parts at frequency of zero spec[3], spec[4] = real & imag parts at frequency of df spec[5], spec[6] = real & imag parts at frequency of 2*df . spec[2*nfreq-1], spec[2*nfreq] = real & imag parts at freq. of (nfreq-1)*df this is the Nyquist frequency (November, 2001, Peter Shearer) */ void getspec(float ts[], long ntime, float dt, int itap, float spec[], long *nfreq, float *df) { void realfft(float data[], long n, int isign); void taper(float ts1[], long n, float ts2[]); float a[MAX_TS_PTS+1]; long i; if (itap != 1) { for (i=1; i<=ntime; ++i) a[i] = ts[i]; } else { taper(ts, ntime, a); } realfft(a, ntime, 1); *nfreq = 1 + ntime/2; *df = 1./(dt*(float)ntime); spec[1] = a[1]; spec[2] = 0.; spec[*nfreq*2-1] = a[2]; spec[*nfreq*2 ] = 0.; for (i=3; i < *nfreq*2-1; ++i) spec[i] = a[i]; } /* getts computes a real time series from a complex spectrum (stored as in the getspec function) using FFT routines that are compatitible with the codes in Numerical Recipes. Requires: realfft, fft MAX_TS_PTS defined as maximum number of time points Inputs: spec = complex spectrum stored as floats nfreq = number of complex frequency points (int, = 1+ntime/2) df = frequency spacing Returns: ts = input time series with values from ts[1] to ts[ntime] ntime = number of time points (int) dt = sample interval (s) (float) spec is assumed to be packed as follows: spec[0] = not used spec[1], spec[2] = real & imag parts at frequency of zero spec[3], spec[4] = real & imag parts at frequency of df spec[5], spec[6] = real & imag parts at frequency of 2*df . spec[2*nfreq-1], spec[2*nfreq] = real & imag parts at freq. of (nfreq-1)*df this is the Nyquist frequency (November, 2001, Peter Shearer) */ void getts(float spec[], long nfreq, float df, float ts[], long *ntime, float *dt){ void realfft(float data[], long n, int isign); void taper(float ts1[], long n, float ts2[]); float a[MAX_TS_PTS+3], scale; long i; *ntime = (nfreq-1)*2; *dt = 1./(df*(float)*ntime); scale = (float)2/(float)*ntime; a[1] = spec[1]*scale; a[2] = spec[nfreq*2-1]*scale; for (i=3; i i) { SWAP(data[j],data[i]); SWAP(data[j+1],data[i+1]); } m = n/2; while (m >= 2 && j > m) { j = j - m; m = m/2; } j += m; } mmax=2; while (n > mmax) { istep = mmax << 1; theta = idir*(6.28318530717959/mmax); wtemp = sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi = sin(theta); wr = 1.0; wi = 0.0; for (m=1; m