Global Index (short | long) | Local contents | Local Index (short | long)
Plot the following: 1. Total time series spectrum 2. Yr 1:2000 3. Yr 2001:4000; 4. Yr 4001:10000;
| This script calls | |
|---|---|
clean
cd ~/model/csiro
load tsu_10k_ann_anomCT.ts
ct = tsu_10k_ann_anomCT(:,2);
yr = tsu_10k_ann_anomCT(:,1);
load ct_yr_1_1000_anom.mat
ct = [ct; detrend(ct2(101:1000))]; yr = [yr; [101:1000]'];
cd ~/matlab/CSIRO/New_calcs/data
load heat_pcs.mat
ct = [ct; std(ct)*rpcs(:,1)]; yr = [yr; [101:1000]'];
ind1 = [1:1000:9001];
ind2 = ind1+999;
lev1 = [];
lev2 = [];
nfft2 = [256 256*[1 1 1] 128];
nfft2 = 128*ones(1,10);
figure(1); fo(1); clf;
for i = 1:10;
ctann = detrend(ct(ind1(i):ind2(i)));
nfft = nfft2(i);
nx = length(ctann); noverlap = 0.75*nfft; n2 = nfft/2;
dofx(i) = round(nx/n2);
lev5 = finv(0.95, dofx(i), 100000);
lev1 = finv(0.99, dofx(i), 100000);
[p, f] = spectrum(ctann, nfft, noverlap, 'none');
f2 = f/2;
% AR1 process
rho = corr(ctann, ctann, 1);
freq = (pi/n2)*(0:n2);
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn(i) = find((freq.*rn)==max(freq.*rn));
% AR2 process
rho1 = corr(ctann, ctann, 1);
rho2 = corr(ctann, ctann, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq2 = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq2)-...
2*phi2*cos(2*freq2)+2*phi1*phi2*cos(freq2));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(:,1));
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
for j = 2:2;
if j == 1;
subplot(5,2,2*i-1);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn, '-k', ...
log(f2), lev5*f2'.*rn, '--k', ...
log(f2), lev1*f2'.*rn, '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
elseif j == 2
subplot(5,2,i);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn2, '-k', ...
log(f2), lev5*f2'.*rn2, '--k', ...
log(f2), lev1*f2'.*rn2, '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
hold on;
ll = line(log(f2(maxrn2(i)))*[1 1], [0 0.005]);
pp = plot(log(f2(maxrn2(i))), 0.005, '^k');
set(ll, 'color', 'k');
hold off;
end
set(hh(1), 'linewidth', 2);
axis([-log(140) log(0.5) -0.001 0.0385]);
set(gca, 'XTick', -log(2.^[7:-1:-1]), ...
'XTickLabel', 2.^[7:-1:-1], 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05);
grid on
set(gca, 'fontsize', 8);
end
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
%print -dps2 AR1_AR2_spect.ps
% Look at original time series
% Look at RAW pcs:
clean
cd /home/disk/tao/dvimont/matlab/CSIRO/New_calcs2/data
load heat_pcs.mat; ct2 = rpcs(:,1);
% Check if removing the 25yr periodicity does anything
[b,a]=butter(11, 2/25);
[b2,a2]=butter(11,2/32);
%ct3 = ct2 - (filtfilt(b,a,ct2)-filtfilt(b2,a2,ct2));
%ct2 = ct2 - filtfilt(b,a,ct2);
ct = [zeros(10000, 1); ct2];
yr = [zeros(10000, 1); [101:1000]'];
% Check if removing the 25yr periodicity does anything
%cd /home/disk/tao/dvimont/model/csiro
%load ct_yr_1_1000_anom.mat
%ct = detrend(ct2(101:1000))./std(detrend(ct2(101:1000)));
%ct = [zeros(10000, 1); ct];
%yr = [zeros(10000, 1); [101:1000]'];
ind1 = [10001 10001 10451]; ind2 = [10900 10450 10900];
nfft2 = [64*[1 1 1]];
figure(2); fo(1); clf;
for i = 1:3;
ctann = detrend(ct(ind1(i):ind2(i)));
nfft = nfft2(i);
nx = length(ctann); noverlap = 0.875*nfft; n2 = nfft/2;
dofx(i) = round(nx/n2);
lev5 = finv(0.95, dofx(i), 100000);
lev1 = finv(0.99, dofx(i), 100000);
[p, f2] = spectrum(ctann, nfft, noverlap, hanning(nfft), 1, 'none');
f = f2;
% AR1 process
rho = (corr(ctann, ctann, 1)+sqrt(corr(ctann, ctann, 2)))/2;
freq = (pi/n2)*(0:n2);
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn(i) = find((freq.*rn)==max(freq.*rn));
sig = find(p(:,1) > lev5*rn')
% AR2 process
rho1 = corr(ctann, ctann, 1);
rho2 = corr(ctann, ctann, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq2 = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq2)-...
2*phi2*cos(2*freq2)+2*phi1*phi2*cos(freq2));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(:,1));
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
sig2 = find(p(:,1) > lev5*rn2')
for j = 1:2;
if j == 1;
subplot(3,2,2*i-1);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn, '-k', ...
log(f2), lev5*f2'.*rn, '--k', ...
log(f2), lev1*f2'.*rn, '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
elseif j == 2
subplot(3,2,2*i);
hh = plot(log(f2), f2.*p(:,1), 'k', ...
log(f2), f2'.*rn2, '-k', ...
log(f2), lev5*f2'.*rn2, '--k', ...
log(f2), lev1*f2'.*rn2, '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
hold on;
ll = line(log(f2(maxrn2(i)))*[1 1], [0 0.1]);
pp = plot(log(f2(maxrn2(i))), 0.1, '^k')
set(ll, 'color', 'k');
hold off;
end
set(hh(1), 'linewidth', 2);
axis([-log(140) log(0.5) -0.05 .65]);
set(gca, 'XTick', -log(2.^[7:-1:-1]), 'XTickLabel', 2.^[7:-1:-1]);
set(gca, 'YTick', 0:.1:0.6);
% axis([-log(140) log(0.5) -0.001 0.0385]);
% set(gca, 'XTick', -log(2.^[7:-1:-1]), 'XTickLabel', 2.^[7:-1:-1]);
% set(gca, 'YTick', 0:.0125:0.05);
grid on
set(gca, 'fontsize', 8);
end
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
%print -dps2 AR1_AR2_spect_101-1000.ps
% Plot 10K run in loglog coordinates
ind1 = [1 1 2001 4001 10001];
ind2 = [10000 2000 4000 10000 10900];
lev1 = [];
lev2 = [];
nfft2 = [512*[1 1 1 1] 128];
figure(1); fo(1); clf;
for i = 1:4;
ctann = detrend(ct(ind1(i):ind2(i)));
nfft = nfft2(i);
nx = length(ctann); noverlap = 0.75*nfft; n2 = nfft/2;
dofx(i) = round(nx/n2);
lev5 = finv(0.95, dofx(i), 100000);
lev1 = finv(0.99, dofx(i), 100000);
[p, f] = spectrum(ctann, nfft, noverlap, 'none');
f2 = f/2;
% AR1 process
rho = corr(ctann, ctann, 1);
freq = (pi/n2)*(0:n2);
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn(i) = find((freq.*rn)==max(freq.*rn));
% AR2 process
rho1 = corr(ctann, ctann, 1);
rho2 = corr(ctann, ctann, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq2 = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq2)-...
2*phi2*cos(2*freq2)+2*phi1*phi2*cos(freq2));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(:,1));
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
for j = 1:2;
if j == 1;
subplot(5,2,2*i-1);
hh = plot(log(f2), log(p(:,1)), 'k', ...
log(f2), log(rn), '-k', ...
log(f2), log(lev5*rn), '--k', ...
log(f2), log(rn/lev5), '--k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
elseif j == 2
subplot(5,2,2*i);
hh = plot(log(f2), log(p(:,1)), 'k', ...
log(f2), log(rn2), '-k', ...
log(f2), log(lev5*rn2), '--k', ...
log(f2), log(rn2/lev5), '--k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
end
set(hh(1), 'linewidth', 2);
axis([-log(550) log(0.5) -5.1 0.1]);
set(gca, 'XTick', -log(2.^[10:-1:-1]), ...
'XTickLabel', 2.^[10:-1:-1], 'YTick', [-5:0]);
grid on
set(gca, 'fontsize', 8);
end
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
%print -dps2 AR1_AR2_spect_loglog.ps
% Look at 1K run in loglog coordinates
ct = getct; ct = detrend(ct)./std(detrend(ct));
ct = [zeros(10000, 1); ct];
yr = [zeros(10000, 1); [101:1000]'];
ind1 = [10001 10001 10451]; ind2 = [10900 10450 10900];
nfft2 = [256*[1 1 1]];
figure(1); fo(1); clf;
for i = 1:3;
ctann = detrend(ct(ind1(i):ind2(i)));
nfft = nfft2(i);
nx = length(ctann); noverlap = 0.75*nfft; n2 = nfft/2;
dofx(i) = round(nx/n2);
lev5 = finv(0.95, dofx(i), 100000);
lev1 = finv(0.99, dofx(i), 100000);
[p, f] = spectrum(ctann, nfft, noverlap, 'none');
f2 = f/2;
% AR1 process
rho = corr(ctann, ctann, 1);
freq = (pi/n2)*(0:n2);
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn(i) = find((freq.*rn)==max(freq.*rn));
% AR2 process
rho1 = corr(ctann, ctann, 1);
rho2 = corr(ctann, ctann, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq2 = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq2)-...
2*phi2*cos(2*freq2)+2*phi1*phi2*cos(freq2));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(:,1));
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
for j = 1:2;
if j == 1;
subplot(5,2,2*i-1);
hh = plot(log(f2), log(p(:,1)), 'k', ...
log(f2), log(rn), '-k', ...
log(f2), log(lev5*rn), '--k', ...
log(f2), log(lev1*rn), '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
elseif j == 2
subplot(5,2,2*i);
hh = plot(log(f2), log(p(:,1)), 'k', ...
log(f2), log(rn2), '-k', ...
log(f2), log(lev5*rn2), '--k', ...
log(f2), log(lev1*rn2), '-.k');
t1 = title(['Years ' num2str(yr(ind1(i))) '--' ...
num2str(yr(ind2(i))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ...
' yrs']);
set(t1, 'fontsize', 9);
hold on;
ll = line(log(f2(6))*[1 1], [0 0.1]);
pp = plot(log(f2(6)), 0.1, '^k')
set(ll, 'color', 'k');
hold off;
end
set(hh(1), 'linewidth', 2);
axis([-log(140) log(0.5) -3 3]);
set(gca, 'XTick', -log(2.^[7:-1:-1]), 'XTickLabel', 2.^[7:-1:-1]);
% set(gca, 'YTick', 0:.1:0.6);
% axis([-log(140) log(0.5) -0.001 0.0385]);
% set(gca, 'XTick', -log(2.^[7:-1:-1]), 'XTickLabel', 2.^[7:-1:-1]);
% set(gca, 'YTick', 0:.0125:0.05);
grid on
set(gca, 'fontsize', 8);
end
end
% Wavelet - type analysis of 10K run
clean
cd ~/model/csiro
load tsu_10k_ann_anomCT.ts
ct = tsu_10k_ann_anomCT(:,2);
yr = tsu_10k_ann_anomCT(:,1);
%load ct_yr_1_1000_anom.mat
%ct = [ct; detrend(ct2(1:1000))]; yr = [yr; [1:1000]'];
ntim = length(ct);
nfft = 1500;
nest = (ntim-nfft)/20+1;
nfft2 = 256; nover = 0.875*nfft2;n2 = nfft2/2;
pp = zeros((nfft2/2+1), nest);
for i = 1:nest
ind = 20*(i-1)+[1:nfft];
[p,f] = spectrum(ct(ind), nfft2, nover, hanning(nfft2), 1, 'none');
pp(:,i) = p(:,1);
end
ctann = ct; nx = length(ctann);
rho1 = corr(ctann, ctann, 1);
rho2 = corr(ctann, ctann, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq2 = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq2)-...
2*phi2*cos(2*freq2)+2*phi1*phi2*cos(freq2));
rn2 = rn2 / mean(rn2);
rn2 = rn2 * mean(p(:,1));
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
rho = corr(ctann, ctann, 1);
freq = (pi/n2)*(0:n2);
rn = (1-rho^2) ./ (1-2*rho*cos(freq)+rho^2);
rn = rn / mean(rn);
rn = rn * mean(p(:,1));
maxrn(i) = find((freq.*rn)==max(freq.*rn));
tind = yr((nfft/2+1):20:(ntim-(nfft/2)+1));
dofx = round((nfft)/n2);
levtop = finv(0.995, dofx, 100000);
levbot = finv(0.995, 100000, dofx);
figure(1); fo(1); clf;
subplot(1,1,1);
surface(log(f), tind, (pp'./(ones(nest,1)*rn)));
shading interp
caxis([1./levbot levtop]);
axis([log(f(2)) log(0.5) 0 10000]);
set(gca, 'XTick', -log(2.^[10:-1:1]), 'XTickLabel', 2.^[10:-1:1])
% axis([f(2) (0.5) 0 10000]);
% set(gca, 'XTick', [0:.05:.5], 'XTickLabel', [0:.05:.5].^-1)