load sunspot.dat
year=sunspot(:,1); sun=sunspot(:,2);
Y=fft(sun);
Y(1)=[];
plot(sun);
title('Original Data:  Sunspots');
keyboard

plot(abs(Y(1:143)))
title('Power Spectrum (Mean removed):  Sunspots');
keyboard

[a,b]=max(abs(Y(1:143)));  %Check that b is 26!  See notes

T=288./(1:143);
plot(T,abs(Y(1:143)));
axis([0 30 0 4300]);  %Only look at the relevant part
title('Power Spectrum versus Years');
keyboard

Y=fft(sun);
[vals,idx]=sort(abs(Y));
idx=idx(end:-1:1);
Ynew=zeros(size(Y));
Ynew(idx(1:10))=Y(idx(1:10));
sun2=ifft(Ynew);
sun2=real(sun2);

plot(1:288,sun,1:288,sun2);
legend('Original','Filtered');