|
Loading LOADSUB ALL FROM "FFT.HTS"
or LOADSUB FROM "MATHLIB.HTS"
or LOADSUB Fft FROM "MATHLIB.HTS"
Usage INTEGER Logn
REAL A(*)
COMPLEX F(*)
CALL Fft(Logn,A(*),F(*))
Description
Fft calculates the discrete Fourier transform of the sequence in the array A and stores the result in the array F. Logn is the base-2 log of the number of points in the sequence. The array A must contain at least 2Logn elements and the array F must contain at least 2Logn-1 elements. If they have more than the required number of elements, the extra elements are ignored and unmodified. The number of elements in A denoted by each permitted value of Logn is shown in the table below:
Logn No. Elements (2Logn)
2 4
3 8
4 16
5 32
6 64
7 128
8 256
9 512
10 1024
11 2048
12 4096
13 8192
14 16384
If the values in A are taken to be values of a continuous signal, a(t), sampled at constant intervals of T (time, distance, or whatever units apply), and if the signal sampled contained no terms at or above the frequency 1/2T, then the coefficients in the array C are the coefficients of the Fourier sine series that describes a(t). A(t) can be reconstructed from the elements of F through the following formula: 
where
If the signal a(t) contains components at or above the frequency 1/2T, the situation is complicated by aliasing, which is explained in most signal processing textbooks.
Some of the more common operations done using discrete Fourier transforms, such as convolution, correlation, filtering, and finding power spectral densities are available as separate CSUBs; see the entries for Autocorrelation, Convolve, Correlate, Filter, Rfilter, and Power_Spectrum for details on their use. The inverse of Fft is performed by the Ifft subroutine. A discrete Fourier transform for complex sequences is computed by the Cfft routine.
Errors
Fft causes an HTBasic error if its arguments are not of the types shown in the USAGE section, above, if Logn is not between 2 and 15, inclusive, if the size of A is smaller than 2Logn, or if the size of F is smaller than 2Logn-1.
Examples
Often, the discrete Fourier transform calculated by the Fft and Cft Math Library subroutines are used to obtain frequency information about continuous signals that have been sampled at discrete points. Such is the case with signal brought into the computer by an A/D converter or by an image scanner. The examples in this section use HTBasic Math Library subroutines to explain some of the pitfalls that arise when a continuous signal is represented by discrete samples and how they are commonly avoided. These pitfalls happen because the samples presented to the computer contain no information about the continuous signal's behavior between the sample points or outside the interval covered by the array of samples.
The frequency spectrum of a continuous, periodic signal, a(t), of period P can be represented by a set of coefficients, ck, such that
The coefficients ck are called the complex Fourier series.
The Nyquist criterion. Consider a periodic sawtooth wave. One period of the wave is shown in the figure 1. If the sampling interval is one period of the waveform, the waveform can be represented in a computer by the 16 samples shown as diamonds in the figure. Note that the samples contain no information about what happens to the wave between samples, so the same samples could represent an infinite number of other waveforms.
The complex Fourier sine/cosine transform of the sawtooth can be calculated or found in a math table. The coefficients, ck, of the sine/cosine transform are
where i is the square root of -1.
The following BASIC program calculates the discrete Fourier transform of the sampled waveform.
10 LOADSUB ALL FROM "FFT.CSB"
20 LOADSUB ALL FROM "WAVEFORM.CSB"
30 INTEGER I
40 REAL A(0:15),Anew
50 COMPLEX B(0:7)
60 CALL Waveform(16.0,1.0,0.,0.,3,A(*))
70 CALL Fft(4,A(*),B(*))
80 FOR I=0 TO 7
90 PRINT I,B(I)
100 NEXT I
110 END
The values produced by the discrete Fourier transform are plotted in figure 2 as diamonds. The first four nonzero values of the complex Fourier sine/cosine transform are plotted in the same figure as squares. Note that the plot is on a logarithmic scale.
The coefficients fk of the discrete Fourier transform are related to those of the complex sine/cosine transform by the following formula:
where "*" represents the complex conjugate operation and N represents the number of samples in the sampling interval (16 in the above example). Therefore, the spectrum calculated by the discrete Fourier transform represents the ordinary Fourier spectrum only when the signal in question has nonzero spectral components for one value of m or n. This is called the Nyquist criterion. Application of this criterion resolves the ambiguity mentioned earlier of which of the infinite number of possible waveforms the samples represent.
Usually, this criterion is satisfied by allowing nonzero spectral components only for m = 0. The Nyquist criterion can then be stated as requiring samples to be taken at twice the highest frequency present in the signal being sampled. Note that the sawtooth waveform used in this example does not obey the Nyquist criterion, since is has spectral components for an infinite number of frequencies.
The following BASIC program produces two continuous waveforms that have the same 16 samples as the sawtooth waveform shown at the beginning of this example. The first waveform has nonzero spectral components for m = 0 and the second for n = 0. The waveforms are plotted after the program listing.
10 LOADSUB ALL FROM "FFT.HTS"
20 LOADSUB ALL FROM "WAVEFORM.HTS"
30 INTEGER I,J
40 REAL A(0:15),Anew1,Anew2,Theta1,Theta2
50 COMPLEX B(0:7)
60 RAD
70 CALL Waveform(16.0,1.0,0.,0.,3,A(*))
80 CALL Fft(4,A(*),B(*))
90 FOR I=0 TO 1024
100 Anew1=0.
120 Anew2=0.
120 FOR J=0 TO 7
130 Theta1=2.0*PI*I*J/1024.0
140 Theta2=2.0*PI*I*(16-J)/1024.0
150 Anew1=Anew1+REAL(B(J))*COS(Theta1)+ IMAG(B(J))*SIN(Theta1)
160 Anew2=Anew2+REAL(B(J))*COS(Theta2)- IMAG(B(J))*SIN(Theta2)
170 NEXT J
180 PRINT I/1024.0*16.0,Anew1,Anew2
190 NEXT I
200 END
Windowing. Consider the function
Obviously, this function contains no frequency components above 14/T, in radians per unit time, distance, or whatever. If the waveform described by this function is sampled every 5T/4096 units, the sample frequency is 24096/5T or 8192/5T radians per unit, which is well above twice the highest frequency component in the waveform, so the Nyquist criterion is satisfied with this sample spacing. A sampling interval of 1024 samples of this waveform is shown in figure 5.
When this waveform is sampled for use in the computer, no information is provided on the behavior before or after the sample interval. When a discrete Fourier transform is done on the test waveform, the discrete Fourier transform algorithm assumes that the data set presented to it represents one period of a periodic waveform. Because the set of samples did not cover exactly one period of the waveform, the discrete Fourier transform connects the last sample in the data set with the first; that is, it computes the transform of a function like that shown by the solid line in figure 6. This introduces high-frequency components in the frequency spectrum. The discrete Fourier transform for this waveform is calculated by the BASIC program below and is plotted in figure 8 on the next page. Figure 7 shows the Fourier sine/cosine transform of the original, continuous waveform for comparison purposes.
10 LOADSUB ALL FROM "FFT.CSB"
20 LOADSUB ALL FROM "POLAR.CSB"
30 INTEGER I,J
40 REAL A(0:1023),Amp(0:511),Phase(0:511),Anew,Theta
50 COMPLEX B(0:511)
60 RAD
70 FOR I=0 TO 1023
80 Theta=2.0*PI*I*1.25/1024.0
90 A(I)=(SIN(Theta)-SIN(3.0*Theta)/9+
SIN(5.0*Theta)/25-SIN(7.0*Theta)/49)*8.0/(PI*PI)
100 NEXT I
110 CALL Fft(10,A(*),B(*))
120 CALL Polar(B(*),"D",Amp(*),Phase(*))
130 FOR I=0 TO 511
140 PRINT I,Amp(I),Phase(I)
150 NEXT I
160 END
One way to reduce the high-frequency components in the spectrum of the sampled signal is to multiply the signal by a window function, such as one of the several provided in the math library. The following BASIC program windows the test function using the Kaiser-Bessel window with parameter 4.0 and calculates the discrete Fourier transform of the windowed waveform. The dashed line in figure 6 shows the windowed function and figure 9 shows the discrete transform of the windowed function.
10 LOADSUB ALL FROM "FFT.HTS"
20 LOADSUB ALL FROM "BESMC.HTS"
30 LOADSUB ALL FROM "POLAR.HTS"
40 INTEGER I
50 REAL A(0:1023),Amp(0:511),Phase(0:511)
60 COMPLEX F(0:511)
70 FOR I=0 TO 1023
80 Theta=2.0*PI*I*1.25/1024.0
90 A(I)=(SIN(Theta)-SIN(3.0*Theta)/9+
SIN(5.0*Theta)/25-SIN(7.0*Theta)/49)*8.0/(PI*PI)
100 NEXT I
110 CALL W_kaiser(A(*),4.0,A(*))
120 FOR I=0 TO 1023
130 PRINT I,A(I)
140 NEXT I
150 CALL Fft(10,A(*),F(*))
160 CALL Polar(F(*),"D",Amp(*),Phase(*))
170 FOR I=0 TO 511
180 PRINT I,Amp(I),Phase(I)
190 NEXT I
200 END
Fig 7. The Fourier sine/cosine transform of the test waveform.
Fig 8. The discrete Fourier transform of the test waveform shown in figure 5.
Fig. 9. The discrete Fourier transform of the windowed function.
Note that the magnitudes of the high-frequency components of the spectrum shown in figure 9 are greatly reduced. Note also the wildly-varying phases of the values shown in figure 9. Windowing functions usually reduce the spurious high-frequency components in a sample at the expense of inaccuracies in the phase.
See Also
Cfft, Convolve, Correlate, Filter, Ifft, Power_spectrum, Rfilter
Notes
If the related Cfft subroutine is applied to a sequence of real values, it outputs twice as many coefficients as the Fft routine. For a real input sequence, the real parts of these coefficients are symmetric about the N/2-1st and N/2th coefficients (beginning subscripts with the 0th coefficient) and the imaginary parts are antisymmetric around these same coefficients. The coefficients output by Fft for the same input sequence are double the first N/2 coefficients output by Cfft except for the 0th or d. c. coefficient, which have the same value.
| 
