FAST FOURIER TRANSFORM (FFT) WITH MAPLE

The command for FFT in Maple is FFT(m,x,y). The FFT command is part of the standard package in Maple 6, the instruction has to be read in with the readlib command in Maple V. In FFT(m,x,y,), m is an integer, x is an array containing the real part of the series that we wish to transform, and y is an array containing the imaginary part of the series. The arrays should contain N= elements. The output of the FFT command is the integer m. After the command is executed x and y will contain, respectively, the real and imaginary part of the discrete Fourier transform, i.e. the initial series is overwritten. The indexing of the arrays is a bit awkward: Maple stores arrays as x=x[1],x[2],x[3]..., that is the first index is 1, while in the formulas for discrete Fourier transform the first entry is . For this reason x[n+1]+I*y[n+1]= while in the output x[n+1]+I*y[n+1]= . Also note the minus sign in the exponent in the equation for the transform.

To check that we got it right

> restart;

If you use Maple V you need to add the command

>readlib(FFT);

> x:=array([1,2,0,0]);y:=array([0,1,0,0]);

Please note that you neead to declare x, and y as arrays! Simply writing

>x:=[1,2,0,0];y:=[0,1,0,0];

will not work. The FFT command will now produce the two arrays x and y which are the real and imaginar parts of the discrete Fourier transform. To print an array you need to use the print command, unlike the case of a list where if you write the name of the list followed by a semi-colon it will print it out.

> FFT(2,x,y);print(x);print(y);

Note that the first element of the Fourier transform has index 1 in the array

> x[1];

The inverse FFT is iFFT.

> iFFT(2,x,y);print(x);print(y);

You should check for yourself that you have understood how the indexing and sign convention works!

As examples of the application of FFT let us calculate discrete fourier series and the power spectrum for some time series

> restart;with(plots):

Warning, the name changecoords has been redefined

Again if you are using Maple V you need to put

>readlib(FFT);

here

Let us pick m=10, so that N= =1024

> m:=10;N:=2^m;

Let us first choose the time series of the function sampled at times t[k]= , k=0..N-1. The sampling frequencies in Hertz

are nu[n]=n/T, or 2*pi*n/T in radians. The fundamental frequency of the signal is then the 50th non-zero sampling frequency

> x:=array([seq(evalf(sin(n*100*Pi/N)),n=0..N-1)]):

> y:=array([seq(0,n=0..N-1)]):

Next, FFT the series

> FFT(m,x,y);

The power spectrum is then fully concentrated in the 50th nonzero frequency as it should:

> listplot([seq([n-1,2*(abs(x[n]+I*y[n]))^2/N],n=1..N/2)],view=[0..100,0..600],
labels=["n","spectral power "]);

THE RANDOM WALK

As our next example let us calculate the power spectrum of a "random walk" in which the value of the n'th term in the time series is incremented from the previous one by a random amount

To be specific we take to be a random number with normal (Gaussian) distribution and variance 1. Such numbers can be generated

using the stats package. If we want a new sequence of random numbers every time we generate the series we can initiate the random number generator by the command randomize(), which output a "seed" from the computer clock. If we want to be able to reproduce the series we initiate the generator by randomize(seed), where the seed is a integer chosen by us.

> restart;with(plots):with(stats):randomize();

Warning, the name changecoords has been redefined

Again if you are using Maple V put

>readlib(FFT);

here

> m:=10;N:=2^m;

Let us first generate N random numbers

> x:=array([stats[random, normald](N)]):

> x[1]:=0;

> for k from 2 to N do
x[k]:=x[k-1]+x[k]:od:

To see what the random walk looks like:

> listplot(x);

> y:=array([seq(0,k=1..N)]):

The discrete Fourier transform assumes the periodic extension of the functions to be transformed. On the other hand the random walk doesn't typically end up where it started. This introduces a spurious discontinuity in the periodic extension of the function, something which will introduce

high frequency components in the Fourier coefficients. Sometimes one avoided this issue by Fourier transforming the "de-trended" series

where the walk starts and ends at the same spot.

> FFT(m,x,y);

The power spectrum is best visualized on a log-log plot omitting the first point

> listplot([seq([log(n-1),log(2*((x[n])^2+(y[n])^2))/N],n=2..N/2)],labels=["log n",
"log(power spectrum) "]);

We note that the on the average the log-log plot of the power spectrum lies along a straight line. By reading off the slope

of the line one finds that the power spectrum is approximately proportional to where f is the frequency. This behavior is typical of the random walk. Indeed, we can generate time series which looks much like the random walk by assuming a fourier transform where the cooefficents

between n=1 and N/2 follow the power law but where the phases are random.

> pi:=evalf(Pi);

> phi:=array([stats[random,uniform[-pi,pi]](N)]):

> x[1]:=0;y[1]:=0;x[N/2+1]:=2/N;y[N/2+1]:=0;

For a real time series the coefficients between N/2 and N are not independent of the coefficients between 1 and N/2. Because of the idiosyncracy

of the Maple indexing mentioned at the beginning, it is a bit tricky to get it right. You may wish to check the formulas yourself!

> for n from 2 to N/2 do
x[n]:=cos(phi[n])/(n-1):
y[n]:=sin(phi[n])/(n-1):
x[N-n+2]:=x[n]:
y[N-n+2]:=-y[n]:od:

> iFFT(m,x,y);

> listplot(x);

You can check by the command listplot(y) that except for very small rounding-off errors this walk is real! Not that because of the periodicity assumed by the discrete Fourier transform the walk will end up where it started. If you don't like this you should add a trend.

The original (non de-trended) walk is special in that the walk after the walker has reached a certain point, the rest of walk is completely

independent of what went on before, there is no memory and the walk cannot be predicted from past data. Let us see

what happens if we let the power spectrum fall off as rather than

> restart;with(plots):with(stats):randomize();m:=10;N:=2^m;

Warning, the name changecoords has been redefined

If you are using Maple V put

>readlib(FFT);

here

> pi:=evalf(Pi);

> phi:=array([stats[random,uniform[-pi,pi]](N)]):

> x[1]:=0;y[1]:=0;x[N/2+1]:=evalf((2/N)^3/2);y[N/2+1]:=0;

> for n from 2 to N/2 do
x[n]:=evalf(cos(phi[n])/(n-1)^(3/2)):
y[n]:=evalf(sin(phi[n])/(n-1)^(3/2)):
x[N-n+2]:=x[n]:
y[N-n+2]:=-y[n]:od:

> x:=array([seq(x[n],n=1..N)]):

> y:=array([seq(y[n],n=1..N)]):

> FFT(m,x,y);

> listplot(x);

This walker is persistent . If the walk is in a certain direction it will tend to continue that way! Let us next consider the

where the power spectrum is proportional to 1/f

> restart;with(plots):with(stats):randomize();m:=10;N:=2^m;

Warning, the name changecoords has been redefined

If you are using Maple V put

> readlib(FFT)

here

> pi:=evalf(Pi);

> phi:=array([stats[random,uniform[-pi,pi]](N)]):

> x[1]:=0;y[1]:=0;x[N/2+1]:=evalf((2/N)^(1/2));y[N/2+1]:=0;

> for n from 2 to N/2 do
x[n]:=evalf(cos(phi[n])/(n-1)^(1/2)):
y[n]:=evalf(sin(phi[n])/(n-1)^(1/2)):
x[N-n+2]:=x[n]:
y[N-n+2]:=-y[n]:od:

> x:=array([seq(x[n],n=1..N)]):

> y:=array([seq(y[n],n=1..N)]):

> FFT(m,x,y);

> listplot(x);

This random walk looks qualitatively quite different! There is "anti-persistance": if the walker has made some steps in a certain

direction, there is a tendency for the next steps to by be in the opposite direction. There is also large scale structures for long times.

1/f noise is seen fairly often in physics, for reasons that are often poorly understood, and there is sometimes a certain amount of mysticism attached to the phenomenon.

>