% 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 \$