% The fundamental tool of signal processing is the FFT, or fast Finite
% Fourier Transform.
%
% Here is an example Using FFT to Calculate Sunspot Periodicity. This file
% uses the FFT function to analyze the variations in sunspot activity over
% the last 300 years.
%
% This file is based on the file "sunspots.m" of Matlab.7 but slightly
% modified to point out the different steps of the simulation.
%
% Sunspot activity is cyclical and reaches a maximum about every 11 years.
% The next plot shows a quantity called the Wolfer number, which measures
% both number and size of sunspots. (Astronomers have tabulated this number
% for almost 300 years).
% sunspot.dat is included in Matlab7.0, however we placed this file in the
% same repertory as the m-file that you are currently reading.
load sunspot.dat
year=sunspot(:,1);
wolfer=sunspot(:,2);
plot(year,wolfer);
title('Sunspot Data');
grid;
disp(' ')
disp('Press Any Key to Continue. ')
disp(' ')
pause
% View of the first 50 years.
plot(year(1:50),wolfer(1:50),'b.-');
grid;
disp(' ')
disp('Press Any Key to Continue. ')
disp(' ')
% Compute the FFT of the sunspot data. The first component of Y, Y(1), is
% simply the sum of the data, and can be removed.
pause
Y = fft(wolfer);
Y(1)=[];
% A plot of the distribution of the Fourier coefficients in the complex
% plane is difficult to interpret.
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');
grid;
disp(' ')
disp('Press Any Key to Continue. ')
disp(' ')
pause
% The complex magnitude squared of Y is called the power, and a plot of power
% versus frequency is a "periodogram".
n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')
grid;
disp(' ')
disp('Press Any Key to Continue. ')
disp(' ')
pause
% The scale in cycles/year can be replaced by the scale years/cycle and we
% can estimate the length of one cycle.
plot(freq(1:40),power(1:40))
xlabel('cycles/year')
% Plot the power versus period for convenience (where period=1./freq).
% There is a very prominent cycle with a length of about 11 years.
period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');
grid;
pause
% Finally, we can fix the cycle length a little more precisely by picking
% out the strongest frequency. The red dot locates this point.
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;
disp(' ')
disp('Press Any Key to Continue. ')
disp(' ')
pause
disp('Period=')
[mp,index] = max(power);
period(index)
% Original file, sunspots.m:
% Copyright 1984-2004 The MathWorks, Inc.
% $Revision: 5.14.4.2 $ $Date: 2004/04/10 23:25:48 $