%Class for CO907.
%__________________________________________
%Author: Jiarui Cao
%Date: 19 Nov 2013
%Purpose: Problem Sheet 2. Q2.7a. From the spectral density we observe the
%11-year sunspot cycle.
%Read more about sunspot cycle at: http://solarscience.msfc.nasa.gov/SunspotCycle.shtml


close all;
clear all;
clc;
clf;


%% Load data
%set file path
% filepath =['/Users/jiaruicao/Dropbox/Teaching/TACO907/Sheet2/data/'];
filepath =['./'];


filename=[filepath 'sunspots.dat'];
delimiter = ' ';startRow = 4;
formatSpec = '%f%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
fclose(fileID);sun = dataArray{:, 1};
clearvars filename delimiter startRow formatSpec fileID dataArray ans;


N=length(sun);
Fs=12; % Frequency. Event/year
sun=sun-mean(sun); %center the data (since there's no obvious trend)
hold on
% PSD by periodogram
[PSDE_p,F_p] = periodogram(sun,rectwin(N),N,Fs);
plot(F_p,10*log10(PSDE_p));
%PSD by pcov
M=30;[PSDE_c,F_c] = pcov(sun,M,N,Fs);
plot(F_c,10*log10(PSDE_c),'-xr','linewidth',2);

%A line by 1/11, since we want to see the 11 year cicular.
line([1/11 1/11], [-30 50],'color','k','linewidth',2,'LineStyle','--');

% set labels & title
title('PSD Estimate. Sunspot.','fontsize',24);
xlabel('Frequency(Sample/Year)','fontsize',24);
ylabel('Power Frequency(dB/Sample/Year)','fontsize',24);
lg2=['Pcov(M=' num2str(M) ')'];
lg=legend('Periodgram',lg2,'1/11');
set(lg,'fontsize',24);
grid on;
hold off
ylim([-30,50]);