%** create plots of power specra for paleoclimate data sets
%   with UNIFORM sample interval
%
%   this script uses quick and easy function f_powerspec 
%   [freq, pow]=f_powerspec(SDATA, SF)
%           
%   SDATA:  equally-spaced (in time) samples
%           
%   SF:     sample frequency
%   freq:   frequency spectrum
%   pow:    signal power at those frequencies
%
%   this script uses function f_autorcorr
%   R=f_autorcorr(SDATA, SF)
%
%   R:      level of correlation at a range of lags 
%           1: lag in time units
%           2: correlation
%
%  detrend the data to clean up the analysis
%  buit-in Matlab detrend takes flags 'constant' and 'linear'


%** SPECMAP data

load specmap_delO18_SST.mat

figure(100)
clf

%* spectral analysis
ts=detrend(SPEC_delO18(:,2),'linear');
sf=SPEC_delO18(2,1)-SPEC_delO18(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);

subplot(3,1,1)
plot(fSPEC(2:101), pSPEC(2:101))
grid on
ylabel('power'), xlabel('frequency (cycles/ka)')
title('spectral analysis of SPECMAP benthic \delta O^{18}')


%* correlation analysis
autoC=f_autocorr(ts,sf);
[fSPAC, pSPAC]=f_powerspec(autoC(:,2), sf);

subplot(3,1,2)
plot(autoC(:,1), autoC(:,2), 'm-')
grid on
ylabel('correlation'), xlabel('lag')

subplot(3,1,3)
plot(fSPAC(2:27), pSPAC(2:27))
grid on
ylabel('power')
xlabel('frequency spectrum of autocorrelation (cycles/ka)')



%** SPECMAP and INSOLOATION

figure(200)
clf

ts=detrend(INSOL_Jun(1:1000,3),'linear');
sf=INSOL_Jun(2,1)-INSOL_Jun(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
subplot(3,1,1)
plot(fSPEC(2:101), pSPEC(2:101))
grid on
ylabel('power'), %xlabel('frequency (cycles/ka)')
title('frequency spectrum of June insolation at 60 N')


ts=detrend(INSOL_Dec(1:1000,3),'linear');
sf=INSOL_Dec(2,1)-INSOL_Dec(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
subplot(3,1,2)
plot(fSPEC(2:101), pSPEC(2:101))
grid on
ylabel('power'), %xlabel('frequency (cycles/ka)')
title('frequency spectrum of December insolation at 60 N')


ts=detrend(SPEC_delO18(:,2),'linear');
sf=SPEC_delO18(2,1)-SPEC_delO18(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
subplot(3,1,3)
plot(fSPEC(2:101), pSPEC(2:101))
grid on
ylabel('power'), xlabel('frequency (cycles/ka)')
title('frequency spectrum of SPECMAP benthic \delta O^{18}')



%** Lisiecki deep sea record
figure(300)
clf

ts=detrend(PliPleHol_s_delO18(1:501,2),'linear');
sf=PliPleHol_s_delO18(2,1)-PliPleHol_s_delO18(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
nd=ceil(length(fSPEC)/4);       % don't need to plot a lot of these

subplot(3,1,1)
plot(fSPEC(2:nd)/sf, pSPEC(2:nd))
grid on
ylabel('power'), %xlabel('frequency (cycles/ka)')
title('frequency spectrum of Lisieki benthic \delta O^{18}')
legend('0 to 0.5 Ma')


ts=detrend(PliPleHol_s_delO18(1:2001,2),'linear');
sf=PliPleHol_s_delO18(2,1)-PliPleHol_s_delO18(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
nd=ceil(length(fSPEC)/4);

subplot(3,1,2)
plot(fSPEC(3:nd)/sf, pSPEC(3:nd))
grid on
ylabel('power'), %xlabel('frequency (cycles/ka)')
%title('power spectrum of Lisieki benthic \delta O^{18}')
legend('0 to 2 Ma')


ts=detrend(PliPleHol_s_delO18(1001:2001,2),'linear');
sf=PliPleHol_s_delO18(2,1)-PliPleHol_s_delO18(1,1);
[fSPEC, pSPEC]=f_powerspec(ts, sf);
nd=ceil(length(fSPEC)/4);

subplot(3,1,3)
plot(fSPEC(3:nd)/sf, pSPEC(3:nd))
grid on
ylabel('power'), xlabel('frequency (cycles/ka)')
%title('power spectrum of Lisieki benthic \delta O^{18}')
legend('1 to 2 Ma')