E. D. Cokelet, NOAA/PMEL

12 Jan 2000

FFTPACK is "a package of Fortran subprograms for the fast Fourier transform of periodic and other symmetric sequences" written by Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, Colorado. FFTPACK can be downloaded for no charge from a national repository of numerical algorithms archived at www.netlib.org. The purpose of these notes is to explain and expand the documentation provided by Swarztrauber in his 4 Apr 1985 version. The 1985 documentation shows how to use the Fortran code, but leaves somewhat unclear what is being calculated and how to interpret the results. (Note: subroutine RFFTF is used by the Ferret external functions ffta and fftp.)

To fix ideas, let us suppose that we have a real time series, r(t), and we wish to compute its Fourier coefficients, i.e. its forward Fourier transform. Swarztrauber’s routine RFFTF takes as input the real sequence of N realizations in time of r(t), i.e.

r(t) = r(iÆt) = r_i+1 i=0,1,..., N-1 (1)

where Æt is the constant time increment. T=NÆt is the fundamental period of the observations such that r₁=r_N+1, etc. (Note: one gives the subroutine N values, r₁ to r_N, and the assumed periodicity implies r₁=r_N+1.) The returned Fourier coefficients will have angular frequencies, w_k = k 2¹/T, k=0, 1, ..., [N/2], which are harmonics of the fundamental angular frequency 2¹/T. (The maximum value of k must be an integer. In the case of integer division only, we adopt the notation [N/2] = "the integer part of N/2." Therefore [32/2]=[33/2]=16.)

The subroutines compute the answer "in place" which means that the output sequence of Fourier coefficients is overwritten into the input seqence of time values. Swarztrauber labels both of these as r_i. To avoid confusion we shall refer to the input sequence (the time series) as r_i and the output sequence (the Fourier coefficients) as R_i.

How do the returned Fourier coefficients relate to the cosine and sine amplitudes? It is instructive to discuss the backward Fourier transform routine, RFFTB. From Swarztrauber’s documentation, but in my notation, RFFTB returns for N even:

which we can write as

where we have used the relation

and have defined

(4)

The factor of N on the left-hand side of (2) enters because "...this transform is unnormalized since a call to RFFTF followed by a call to RFFTB will multiply the input sequence by N." At first glance, (3) seems to imply there are N+1 unique Fourier coefficients for N points in the time series. This is not true because b_N/2 is arbitrary since it is multiplied by

sin[(i-1)¹]=0. It is wise to set b_N/2 to zero to avoid spurious terms when Fourier series are summed later. For N odd RFFTB returns

for which we can write

with

(7)

For compactness we can combine the even and odd results to give

What time series are we approximating by the discrete Fourier series in (8)? The answer can be seen as follows. If we let t_i=iÆt and since T=NÆt then (8) can be written as

or if we replace i-1 by i, we have

Equation (10) is the discrete approximation to

which itself is the truncated approximation to the infinite Fourier series.

If the input time series starts at t=t_s rather than t=0, (11) can be written as

which means that the user inputs

to subroutine RFFTF.

How do the time series look in amplitude and phase notation? Equation (10) can be written as

where

What do the forward Fourier transforms computed by RFFTF look like? From Swarztrauber’s documentation, but again in my notation, RFFTF returns:

(15)

(16)

and if N is even

(17)

Therefore the even-numbered R_j’s are the cosine transforms, and the odd-numbered R_j’s are the sine tranforms.