% 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)))