% A program to solve the first part of Exercise 7,
% analyzing the Fourier components of the sunspot data.
% DPO 11/2007
% Load the data into array x.

sunspot

% Let t be the vector of years and s be the vector of sunspot numbers.

t = x(:,1);
s = x(:,2);

n = length(s);

% Plot the data.

figure(1)
plot(t,s,'b')
xlabel('year')
ylabel('sunspot number')
title('Sunspot number; data from NOAA')

% Take the discrete Fourier transform of the data and plot it.
% (DCT is more interesting, if you have the function on your
% Matlab installation.)

Fs = fft(s);
figure(2)
title('FFT of sunspot data')
plot([0:n-1],abs(Fs))

% Which periods are significant in the data?
% Find the indices of the largest components.

aFs = abs(Fs);
[a,order] = sort(aFs);
disp(sprintf('The top 10 components in the FFT are:'))
disp(sprintf('Period      Magnitude                '))
disp([order(1:10),a(order(1:10))])

% Note that there are many small components.  Let's try to
% zero them out.  This has several applications:
%
%  1.  Seeing if the data supports the hypothesis that there
%      are only a few significant cyclical components in the data.
%
%  2.  Compressing the data, by only saving the nonzeros.
%      This is used, for example, in transmission of cellphone
%      data.
%
%  3.  Removing noise by zeroing components that correspond to noise.


tol = max(aFs)*1.e-2;
index = find(aFs < tol);

Fs(index) = 0;

disp(sprintf('I chose the drop tolerance for components to be' ))
disp(sprintf('.001 times the magnitude of the largest component.'))
disp(sprintf('This left %d nonzero entries.',n-length(index)))

% Transform back and see how well the modified data fits the original.

shat = ifft(Fs);
figure(3)
hold on
plot(t,s,'b',t,shat,'r')
legend('Original data','Modified data')
xlabel('year')
ylabel('sunspot number')
title('Sunspot number: original and truncated FFT approximation')
disp(sprintf('The relative error in the data after truncation is %f',...
           norm(s-shat)/norm(s)))