SETTING UP THE FAST FOURIER TRANSFORM

Typically we Fourier transform seismograms about a thousand points long. Under these conditions another Fourier summation method works about a hundred times faster than those already given. Unfortunately, the faster Fourier transform program is not so transparently clear as the programs given earlier. Also, it is slightly less flexible. The speedup is so overwhelming, however, that the fast program is always used in routine work.

Flexibility may be lost because the basic fast program works with complex-valued signals, so we ordinarily convert our real signals to complex ones (by adding a zero imaginary part). More flexibility is lost because typical fast FT programs require the data length to be an integral power of 2. Thus geophysical datasets often have zeros appended (a process called ``zero padding") until the data length is a power of 2. From time to time I notice clumsy computer code written to deduce a number that is a power of 2 and is larger than the length of a dataset. An answer is found by rounding up the logarithm to base 2. The more obvious and the quicker way to get the desired value, however, is with the simple Fortran function pad2(). 

integer function pad2( n )
integer n
pad2 = 1
while( pad2 < n )
        pad2 = pad2 * 2 
return; end

How fast is the fast Fourier transform method? The answer depends on the size of the data. The matrix times vector operation in (3) requires N² multiplications and additions. That determines the speed of the slow transform. For the fast method the number of adds and multiplies is proportional to $N \,\log_2 N$.Since 2¹⁰=1024, the speed ratio is typically 1024/10 or about 100. In reality, the fast method is not quite that fast, depending on certain details of overhead and implementation. In 1987 I tested the three programs on a 1024-point real signal and found times

 		slowft 		 153s
		polyft 		 36s
		ftu 		 .7s

Below is ftu(), a version of the fast Fourier transform program.

There are many versions of the program--I have chosen this one for its simplicity. Considering the complexity of the task, it is remarkable that no auxiliary memory vectors are required; indeed, the output vector lies on top of the input vector. To run this program, your first step might be to copy your real-valued signal into a complex-valued array. Then append enough zeros to fill in the remaining space. 

subroutine ftu( signi, nx, cx )
#   complex fourier transform with unitary scaling
#
#               1         nx          signi*2*pi*i*(j-1)*(k-1)/nx
#   cx(k)  =  ---- * sum cx(j) * e
#             sqrt(nx)   j=1             for k=1,2,...,nx=2**integer
#
integer nx, i, j, k, m, istep, pad2
real    signi, scale, arg
complex cx(nx), cmplx, cw, cdel, ct
if( nx != pad2(nx) )    call erexit('ftu: nx not a power of 2')
scale = 1. / sqrt( 1.*nx)
do i= 1, nx
        cx(i) = cx(i) * scale
j = 1;  k = 1
do i= 1, nx {
        if (i<=j) { ct = cx(j); cx(j) = cx(i); cx(i) = ct }
        m = nx/2
        while (j>m && m>1) { j = j-m; m = m/2 }         # "&&" means .AND.
        j = j+m
        }
repeat {
        istep = 2*k;   cw = 1.;   arg = signi*3.14159265/k
        cdel = cmplx( cos(arg), sin(arg))
        do m= 1, k {
                do i= m, nx, istep
                        { ct=cw*cx(i+k);  cx(i+k)=cx(i)-ct;  cx(i)=cx(i)+ct }
                cw = cw * cdel
                }
        k = istep
        if(k>=nx) break
	}
return; end

The following two lines serve to Fourier transform a vector of 1024 complex-valued points, and then to inverse Fourier transform them back to the original data:

call ftu(  1., 1024, cx)
call ftu( -1., 1024, cx)

An engineering reference given at the end of this chapter contains many other versions of the FFT program. One version transforms real-valued signals to complex-valued frequency functions in the interval $0 \le \omega < \pi$.Others that do not transform data on top of itself may be faster with specialized computer architectures.

EXERCISES:

1. Consider an even time function that is constant for all frequencies less than $\omega_0$ and zero for all frequencies above $\omega_0$.What is the rate of decay of amplitude with time for this function?
2. Waves spreading from a point source decay in energy as the area on a sphere. The amplitude decays as the square root of energy. This implies a certain decay in time. The time-decay rate is the same if the waves reflect from planar interfaces. To what power of time t do the signal amplitudes decay? For waves backscattered to the source from point reflectors, energy decays as distance to the minus fourth power. What is the associated decay with time?