| 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'
|
