Global Index (short | long) | Local contents | Local Index (short | long)
Look at power spectra of successive 1000 year chunks
| 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);
[b, a] = butter(9, 2/9);
lpct = filtfilt(b, a, ct);
for i = 1:10;
subplot(10,1,i);
ind = 1000*(i-1)+[1:1000];
h = plot(yr(ind), ct(ind), 'k', yr(ind), lpct(ind), 'k');
set(h(2), 'linewidth', 2);
axis([yr(ind(1)) yr(ind(1000)) -1.1 1.1])
grid on;
set(gca, 'fontsize', 8, 'YTick', -1:.5:1);
end
subplot(10,1,1);
t1 = title('10K CSIRO: CT index');
set(t1, 'fontsize', 10);
subplot(10,1,10);
xl = xlabel('Year');
set(xl, 'fontsize', 10);
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
print -dps2 10K_CTindex.ps
% AR1 Process
figure(1); fo(1); clf;
nfft = 128; noverlap = 3*nfft/4; n2 = nfft/2;
nx = 1000;
for i = 1:10;
ind = 1000*(i-1)+[1:1000];
[p, f] = spectrum(ct(ind), nfft, noverlap);
f2 = f/2;
% AR1 process
rho = 0.5*(corr(ct(ind), ct(ind), 1) + ...
sqrt(corr(ct(ind), ct(ind), 2)));
% rho = corr(ct(ind), ct(ind), 1);
rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
rn = rn / mean(rn(3:(n2+1)));
rn = rn * mean(p(3:(n2+1),1));
maxrn(i) = find((f.*rn')==max(f.*rn'));
% AR2 process
rho1 = corr(ct(ind), ct(ind), 1);
rho2 = corr(ct(ind), ct(ind), 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2(3:(n2+1)));
rn2 = rn2 * mean(p(3:(n2+1),1));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
dofx = nx / n2;
ind2 = 2:(n2+1);
subplot(5,2,i)
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(ind(1))) '-' ...
num2str(yr(ind(1000))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ' yrs']);
set(t1, 'fontsize', 9);
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
print -dps2 10K_AR1_spect.ps
% AR2 Process
figure(1); fo(1); clf;
nfft = 128; noverlap = 3*nfft/4; n2 = nfft/2;
nx = 1000;
for i = 1:10;
ind = 1000*(i-1)+[1:1000];
[p, f] = spectrum(ct(ind), nfft, noverlap);
f2 = f/2;
% AR1 process
rho = 0.5*(corr(ct(ind), ct(ind), 1) + ...
sqrt(corr(ct(ind), ct(ind), 2)));
% rho = corr(ct(ind), ct(ind), 1);
rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
rn = rn / mean(rn(3:(n2+1)));
rn = rn * mean(p(3:(n2+1),1));
maxrn(i) = find((f.*rn')==max(f.*rn'));
% AR2 process
rho1 = corr(ct(ind), ct(ind), 1);
rho2 = corr(ct(ind), ct(ind), 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2(3:(n2+1)));
rn2 = rn2 * mean(p(3:(n2+1),1));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
dofx = nx / n2;
ind2 = 2:(n2+1);
subplot(5,2,i)
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn2(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn2(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn2(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(ind(1))) '--' ...
num2str(yr(ind(1000))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ' yrs']);
set(t1, 'fontsize', 9);
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
print -dps2 10K_AR2_spect.ps
% Both, for first and second half of data set
figure(1); fo(1); clf;
nfft = 256; noverlap = 3*nfft/4; n2 = nfft/2;
for i = 1:1;
if 0;
if i == 1; ind = 1:2000; lev = 1.67;
elseif i == 2; ind = 2001:4000; lev = 1.67;
elseif i == 3; ind = 4001:10000; lev = 1.36;
elseif i == 4; ind = 1:5000; lev = 1.39;
elseif i == 5; ind = 5001:10000; lev = 1.39;
end
else; ct = tem2; ind = 1:100000; end
nx = length(ind);
[p, f] = spectrum(ct(ind), nfft, noverlap);
f2 = f/2;
% AR1 process
% rho = 0.5*(corr(ct(ind), ct(ind), 1) + ...
% sqrt(corr(ct(ind), ct(ind), 2)));
rho = corr(ct(ind), ct(ind), 1);
rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
rn = rn / mean(rn(3:(n2+1)));
rn = rn * mean(p(3:(n2+1),1));
maxrn(i) = find((f.*rn')==max(f.*rn'));
% AR2 process
rho1 = corr(ct(ind), ct(ind), 1);
rho2 = corr(ct(ind), ct(ind), 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2(3:(n2+1)));
rn2 = rn2 * mean(p(3:(n2+1),1));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
dofx = nx / n2;
ind2 = 2:(n2+1);
subplot(5,2,(2*i)-1)
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ' yrs']);
set(t1, 'fontsize', 9);
subplot(5,2,(2*i))
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn2(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn2(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn2(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ' yrs']);
set(t1, 'fontsize', 9);
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Longrun/Figs
print -dps2 10K_AR1_AR2_spect.ps
% Look at random AR1 process:
tem = randn(1, 100000);
tem2 = zeros(1, 100000);
for i = 2:100000;
tem2(i) = 0.6*tem2(i-1)+.8*tem(i);
end
% AR2 - entire record
% AR2 Process
figure(1); fo(1); clf;
nfft = 256; noverlap = 3*nfft/4; n2 = nfft/2;
nx = 5000;
for j = 1:2;
for i = 1:2;
ind = nx*(i-1)+[1:nx];
[p, f] = spectrum(ct(ind), nfft, noverlap);
f2 = f/2;
% AR1 process
rho = 0.5*(corr(ct(ind), ct(ind), 1) + ...
sqrt(corr(ct(ind), ct(ind), 2)));
% rho = corr(ct(ind), ct(ind), 1);
rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
rn = rn / mean(rn(3:(n2+1)));
rn = rn * mean(p(3:(n2+1),1));
maxrn(i) = find((f.*rn')==max(f.*rn'));
% AR2 process
rho1 = corr(ct(ind), ct(ind), 1);
rho2 = corr(ct(ind), ct(ind), 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2(3:(n2+1)));
rn2 = rn2 * mean(p(3:(n2+1),1));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
dofx = nx / n2;
ind2 = 2:(n2+1);
subplot(5,2,i+2*(j-1))
% New - f*p vs log(f)
if j == 1;
hh = plot(log(f2(ind2)), f2(ind2)'.*rn(ind2), '-k', ...
log(f2(ind2)), 1.39*f2(ind2)'.*rn(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn(ind2)/1.39, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Period = ' ...
num2str(round(100./(f2(maxrn(i))))/100) ' yrs']);
else
hh = plot(log(f2(ind2)), f2(ind2)'.*rn2(ind2), '-k', ...
log(f2(ind2)), 1.39*f2(ind2)'.*rn2(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn2(ind2)/1.39, '--k', ...
log(f2(ind2)), f2(ind2).*p(ind2,1), '-k');
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/100) ' yrs']);
end
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
set(t1, 'fontsize', 9);
end
end
% Average over frequency bands
figure(1); fo(1); clf;
for i = 1:5;
if 1;
if i == 1; ind = 1:2000; lev = 1.67;
elseif i == 2; ind = 2001:4000; lev = 1.67;
elseif i == 3; ind = 4001:10000; lev = 1.36;
elseif i == 4; ind = 1:5000; lev = 1.39;
elseif i == 5; ind = 5001:10000; lev = 1.39;
end
else; ct = tem2; ind = 1:100000; end
nx = length(ind); nfft = nx; noverlap = 3*nfft/4; n2 = nfft/2;
[p, f] = spectrum(ct(ind), nfft, noverlap);
f2 = f/2;
% AR1 process
% rho = 0.5*(corr(ct(ind), ct(ind), 1) + ...
% sqrt(corr(ct(ind), ct(ind), 2)));
rho = corr(ct(ind), ct(ind), 1);
rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
rn = rn / mean(rn(3:(n2+1)));
rn = rn * mean(p(3:(n2+1),1));
maxrn(i) = find((f.*rn')==max(f.*rn'));
% AR2 process
rho1 = corr(ct(ind), ct(ind), 1);
rho2 = corr(ct(ind), ct(ind), 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);
freq = (pi/n2)*(0:n2);
rn2 = 1./(1+phi1.^2+phi2.^2-2*phi1*cos(freq)-...
2*phi2*cos(2*freq)+2*phi1*phi2*cos(freq));
rn2 = rn2 / mean(rn2(3:(n2+1)));
rn2 = rn2 * mean(p(3:(n2+1),1));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
maxrn2(i) = find((f.*rn2')==max(f.*rn2'));
% rner1 = rn * 2.05; rner5 = rn * 1.68;
dofx = nx / n2;
ind2 = 2:(n2+1);
subplot(5,2,(2*i)-1)
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*rave(p(ind2,1), nfft/100), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Efold = ' ...
num2str(round(100./(2*pi*f2(maxrn(i))))/100) ' yrs']);
set(t1, 'fontsize', 9);
subplot(5,2,(2*i))
% New - f*p vs log(f)
hh = plot(log(f2(ind2)), f2(ind2)'.*rn2(ind2), '-k', ...
log(f2(ind2)), 1.67*f2(ind2)'.*rn2(ind2), '--k', ...
log(f2(ind2)), f2(ind2)'.*rn2(ind2)/1.67, '--k', ...
log(f2(ind2)), f2(ind2).*rave(p(ind2,1), nfft/200), '-k');
set(hh(4), 'linewidth', 2);
axis([-log(140) log(0.55) -0.001 0.0376]);
xtl = ['128'; ' 64'; ' 32'; ' 16'; ' 8 '; ' 4 '; ...
' 2 '; ' 1 '; '0.5'];
set(gca, 'XTick', -log([128 64 32 16 8 4 2 1 0.5]), ...
'XTickLabel', xtl, 'YTick', [0:.025:.15]);
set(gca, 'YTick', 0:.0125:0.05, 'YTickLabel', ['']);
grid on
set(gca, 'fontsize', 8);
t1 = title(['Years ' num2str(yr(min(ind))) '--' ...
num2str(yr(max(ind))) '; Period = ' ...
num2str(round(100./(f2(maxrn2(i))))/200) ' yrs']);
set(t1, 'fontsize', 9);
end