function slide=sunspots
% This is a slideshow file for use with playshow.m and makeshow.m
% To see it run, type 'playshow sunspots', 

%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 5.11 $

if nargout<1,
  playshow sunspots
else
  %========== Slide 1 ==========

  slide(1).code={
   'load sunspot.dat',
   'year = sunspot(:,1);',
   'wolfer = sunspot(:,2);',
   'plot(year,wolfer)',
   ' title(''Sunspot Data'')' };
  slide(1).text={
   ' Press the "Start" button to see a demonstration of the',
   ' FFT function in MATLAB used to analyze the variations',
   ' in sunspot activity over the last 300 years.'};

  %========== Slide 2 ==========

  slide(2).code={
   'plot(year,wolfer)',
   ' title(''Sunspot Data'')' };
  slide(2).text={
   'You are probably aware that sunspot activity is cyclical,',
   ' reaching a maximum about every 11 years. Let''s confirm that.',
   ' Here is a plot of a quantity called the Wolfer number, which',
   ' measures both number and size of sunspots. Astronomers',
   ' have tabulated this number for almost 300 years.',
   ' >> load sunspot.dat',
   ' >> year=sunspot(:,1); wolfer=sunspot(:,2);',
   ' >> plot(year,wolfer)',
   ' >> title(''Sunspot Data'')'};

  %========== Slide 3 ==========

  slide(3).code={
   'plot(year(1:50),wolfer(1:50),year(1:50),wolfer(1:50),''ro'');' };
  slide(3).text={
   ' Here is a closer look at the first 50 years.',
   '',
   ' >> plot(year(1:50),wolfer(1:50),year(1:50),wolfer(1:50),''ro'');'};

  %========== Slide 4 ==========

  slide(4).code={
   'plot(year(1:50),wolfer(1:50),year(1:50),wolfer(1:50),''ro'');' };
  slide(4).text={
   ' The fundamental tool of signal processing is the FFT, or fast',
   ' Finite Fourier Transform. To take the FFT of the sunspot data',
   ' we would type',
   '  >> Y = fft(wolfer);',
   ' The first component of Y, Y(1), is simply the sum of the data,',
   ' and can be removed.',
   ' >> Y(1)=[ ];'};

  %========== Slide 5 ==========

  slide(5).code={
   'Y = fft(wolfer);',
   'n = length(Y);',
   '% The first component, Y(1), is simply the sum of the data',
   '% and can be removed.',
   'Y(1) = [];',
   'plot(Y,''ro'')',
   'title(''Fourier Coefficients in the Complex Plane'');',
   'xlabel(''Real Axis'');',
   'ylabel(''Imaginary Axis'');',
   'set(gca,''YTick'',[],''XTick'',[]);' };
  slide(5).text={
   ' A graph of the distribution of the Fourier coefficients (given',
   ' by Y) in the complex plane is pretty, but difficult to interpret.',
   ' >> plot(Y,''ro'')',
   ' >> title(''Fourier Coefficients in the Complex Plane'');',
   ' >> xlabel(''Real Axis'');',
   ' >> ylabel(''Imaginary Axis'');',
   ' >> set(gca,''YTick'',[],''XTick'',[]);', 
   ' ',
   ' We need a more useful way of examining the data in Y.'};

  %========== Slide 6 ==========

  slide(6).code={
   'power = abs(Y(1:n/2)).^2;',
   'nyquist = 1/2;',
   'freq = (1:n/2)/(n/2)*nyquist;',
   'plot(freq,power)',
   'xlabel(''cycles/year'')',
   'title(''Periodogram'')' };
  slide(6).text={
   ' 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:n/2)).^2;',
   ' >> nyquist = 1/2;',
   ' >> freq = (1:n/2)/(n/2)*nyquist;',
   ' >> plot(freq,power)',
   ' >> xlabel(''cycles/year'')',
   ' >> title(''Periodogram'')',
   ''};

  %========== Slide 7 ==========

  slide(7).code={
   'plot(freq(1:40),power(1:40))',
   'xlabel(''cycles/year'')' };
  slide(7).text={
   'Zoom in on the first 40 components for a closer look.',
   '',
   ' >> plot(freq(1:40),power(1:40))',
   ' >> xlabel(''cycles/year'')',
   ' ',
   'The scale in cycles/year is somewhat inconvenient. We can',
   ' plot in years/cycle and esimate what one cycle is.'};

  %========== Slide 8 ==========

  slide(8).code={
   'period=1./freq;',
   'plot(period,power);',
   'axis([0 40 0 2e+7]);',
   'ylabel(''Power'');',
   'xlabel(''Period (Years/Cycle)'');' };
  slide(8).text={
   ' Here we plot power versus period for convenience (where',
   ' period=1./freq). As expected, 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)'');'};

  %========== Slide 9 ==========

  slide(9).code={
   'hold on;',
   'index=find(power==max(power));',
   'mainPeriodStr=num2str(period(index));',
   'plot(period(index),power(index),''r.'', ''MarkerSize'',25,''EraseMode'',''none'');',
   'text(period(index)+2,power(index),[''Period = '',mainPeriodStr], ''EraseMode'',''none'');',
   'hold off;' };
  slide(9).text={
   ' Finally, we can fix the cycle length a little more precisely, as',
   ' shown in this plot.',
   ' >> hold on;',
   ' >> index=find(power==max(power));',
   ' >> mainPeriodStr=num2str(period(index));',
   ' >> plot(period(index),power(index),''r.'', ...',
   '         ''MarkerSize'',25,''EraseMode'',''none'');',
   ' >> text(period(index)+2,power(index), ...'
   '         [''Period = '',mainPeriodStr],''EraseMode'',''none'');',
   ' >> hold off;' };
end
