16: Analyze data using FFT

Use the FFT (Fast Fourier Transform) to analyze and filter data

1. Read sunspot data s from file sunspot.dat
2. Calculate the FFT of s, and and graph abs(FFT) -- the FFT itself is complex
3. graph the FFT ``power'' of s, (abs(FFT))2, vs. frequency, and vs. period

Mathematica:

(1)
data = ReadList["data/sunspot.dat", ];
year = data[[All,1]]; s = data[[All,2]];

(2)
ffts = Fourier[s];
ListPlot[Abs[ffts], PlotLabel->"Abs[ffts]"];

(3)
len = Length[s]/2
freq = 0.5 Range[len] / len;
power = (Abs[ffts]^2)[[Range[2,len+1]]];
ListPlot[ Transpose[],
    PlotJoined->True, PlotRange->All,
    AxesLabel->];
period=1/freq;
ListPlot[ Transpose[],
    PlotJoined->True, PlotRange->,All},
    AxesLabel->,
    PlotLabel->"Periodogram, years/cycle"];

Matlab:

(1)
load sunspot.dat
year = sunspot(:,1);
s = sunspot(:,2);

(2)
ffts = fft(s);
plot(abs(ffts))

(3)
n=length(ffts);
power = abs(ffts(2:n/2)).^2;
freq = (2:n/2)/n;
plot(freq,power), grid on
title('Periodogram, power vs cycles/year')

(4)
period = 1./freq;
plot(period,power), axis([0,40,0,2e7]), grid on
title('Periodogram, power vs years/cycle')

Maple:

(1)
xxx

(2)
xxx

(3)
xxx

(4)
xxx

IDL:

(1)
len=288 & data=fltarr(2,len)
OpenR, 27, 'sunspot.dat'
ReadF, 27, data
Close, 27
year=data[0,*] & s=data[1,*]

(2)
ffts = fft(s,1)
plot, abs(ffts), title='abs(ffts)'

(3)
len = len/2
freq = (findgen(len)+1)/len/2
power = (abs(ffts[1:len]))^2
plot,freq,power,title='Power vs cycles/year'

(4)
period = 1./freq
plot,period,power,xrange=[0,40], $
    xtitle='years/cycle',ytitle='Power'