functions:FFT
syntax: FFT( y,`COS&SIN' )
The Fourier coefficients are returned. Suppose that the cosine
coefficients are H, and the sine coefficients are G, then H[1]/2
is the mean value of the input data. These coefficients can be
used for smooth interpolation. Suppose xi[j] is the interpolation
location, and 2N is the number of original data points.
yi[j]=H[1]/2+SUM(H[k]*COS((k-1)*xi[j])+G[k]*SIN((k-1)*xi[j]),k,2:N+1)
See the example. The amplitudes, A, and the phases, P, are related to
H and G: A[1]=H[1], A[N+1]=ABS(H[N+1]), A[j]=sqrt(H[j]^2+G[j]^2) for
j=2,3,...,N, and P[1]=0, P[N+1]={0 if H[N+1]>0, 180 if H[N+1]<0},
P[j] = arctan(G[j]/H[j]) and if P[j] < 0 then P[j]=P[j]+360 for
j=2,3,...,N. Conversly, H[j] = A[j]*cos(P[j]) and G[j] = A[j]*sin(P[j])
for j=1,2,...,N+1
The input vector y is assumed to be periodic with an even number of
points. That is, if y is a vector with length N, then it is assumed
that y[1] = y[N+1] = y[2*N+1] = ...
Suppose that N is the length of y. The largest prime factor of N
must not be greater than 23. There must not be more than 11 distinct
prime factors of N. The product of the square-free prime factors of
N must be <= 210. The speed is enhanced by using a value of N with
small prime factors, particularly powers of 2. The PFACTORS function
returns the prime factors of a constant or scalar.
N = 16 ! generate some `data'
X = [0:N-1]*2*PI/N !
Y = SIN(X)+5*RAN(X) !
M=FFT(Y,`COS&SIN') ! calculate Fourier coefficients
XI = [0:2*PI:.05] ! indep. variable for smoothed function
SCALAR\DUMMY K ! define K to be dummy variable for SUM
SET PCHAR 0 ! no plotting symbols
YI=M[1,1]/2+SUM(M[K,1]*COS((K-1)*XI)+M[K,2]*SIN((K-1)*XI),K,2:N/2+1)
GRAPH\NOAXES XI YI
GRAPH\NOAXES [0;2*PI] [M[1,1]/2;M[1,1]/2] ! overlay the mean value