%Class for CO907.
%__________________________________________
%Author: Jiarui Cao
%Date: 19 Nov 2013
%Purpose: Problem Sheet 2. Q2.6 Compute power spectrum density estimate from detreded real data. 


close all;
clear all;
clc;
clf;


%% Load data
%set file path
%filepath =['/Users/jiaruicao/Dropbox/Teaching/TACO907/Sheet2/data/'];
%load data
%filename1=[filepath 'ftse_detrended.dat'];
% filename1=['ftse_detrended.dat'];
% Data=load(filename1);
% Day = Data(:, 1);
% FTSE =Data(:, 2);

filename2=['temperatureanomaly_detrended.dat'];
Data=load(filename2);
Temp=Data;

% filename3=['xray_detrended.dat'];
% Data=load(filename3);
% t = Data(:, 1);
% xray =Data(:, 2);


% %% PSD: xray
% %Fs=length(xray)/(max(t)-min(t)); %Frequency. Event/unit Time.
% Fs=1;
% hold on
% [PSDE_p,F_p] = periodogram(xray,rectwin(length(xray)),length(xray),Fs);
% plot(F_p,10*log10(PSDE_p));
% %plot(F_p,PSDE_p);
% 
% M=2;
% [PSDE_c,F_c] = pcov(xray,M,length(xray),Fs);
% plot(F_c,10*log10(PSDE_c),'r','linewidth',2);
% %plot(F_c,PSDE_c,'r','linewidth',2);
% grid on;
% hold off
% 
% title('PSD Estimate. X-Ray.','fontsize',24);
% xlabel('Frequency(Sample/UnitTime)','fontsize',24);
% ylabel('Power Frequency(dB/Sample/UnitTime)','fontsize',24);
% lg2=['Pcov(M=' num2str(M) ')'];
% lg=legend('Periodgram',lg2);
% set(lg,'fontsize',24);
% %ylim([-30,40])


% model =arima(2,0,0); %guess this is a AR(2) model. 
% 
% estimate(model,xray) %estimate the data


% PSD: Temperature
N=length(Temp);
%T=N/12;
Fs=12; % Frequency. Event/Month

hold on
[PSDE_p,F_p] = periodogram(Temp,rectwin(N),N,Fs);
plot(F_p,10*log10(PSDE_p));

M=10;
[PSDE_c,F_c] = pcov(Temp,M,N,Fs);
plot(F_c,10*log10(PSDE_c),'r','linewidth',2);

title('PSD Estimate. Temperature.','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);
set(lg,'fontsize',24);
grid on;
hold off
ylim([-10,40]);

% % PSD: FTSE
% N=length(FTSE);
% Fs=1; % Frequency. Event/Trading Day
% 
% hold on
% [PSDE_p,F_p] = periodogram(FTSE,rectwin(N),N,Fs);
% plot(F_p,10*log10(PSDE_p));
% 
% M=3;
% [PSDE_c,F_c] = pcov(FTSE,M,N,Fs);
% plot(F_c,10*log10(PSDE_c),'r','linewidth',2);
% 
% 
% title('PSD Estimate. FTSE.','fontsize',24);
% xlabel('Frequency(Sample/TradingDay)','fontsize',24);
% ylabel('Power Frequency(dB/Sample/TradingDay)','fontsize',24);
% lg2=['Pcov(M=' num2str(M) ')'];
% lg=legend('Periodgram',lg2);
% set(lg,'fontsize',24);
% grid on;
% hold off
% ylim([-10,80]);