Documentation of test_AR2_spect

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

cd ~/model/enso
load tfield1;
yr = tfield1(:,1)./360;
ct = tfield1(:,2);
ct = annave(ct);
ct2 = ct(1:120000);

Cross-Reference Information

Listing of script test_AR2_spect

clean
cd ~/model/csiro
load tsu_10k_ann_anomCT.ts
ct = tsu_10k_ann_anomCT(:,2);
yr = tsu_10k_ann_anomCT(:,1);
ct3 = detrend(ct);
load ct_yr_1_1000_anom.mat
ct2 = detrend(ct2);

ntim = length(ct2);
ct3 = NaN*ones(ntim/3-1, 1);
for i = 1:(ntim/3-1);
  ind = 3*(i-1)+[3:5];
  ct3(i) = mean(ct2(ind));
end

rho1 = corr(ct3, ct3, 1);
rho2 = corr(ct3, ct3, 2);
phi1 = (rho1 - rho2*rho1)./(1-rho1.^2);
phi2 = (rho2 - rho1.^2)./(1-rho1.^2);

x = randn(1, 120000);
y = zeros(size(x));
for i = 3:length(y);
  y(i) = phi1*y(i-1)+phi2*y(i-2)+x(i);
end

nfft = 12*256;
nx = length(y); noverlap = 0.75*nfft; n2 = nfft/2;
dofx = round(nx/n2);
lev5 = finv(0.95, dofx, 100000);
lev1 = finv(0.99, dofx, 100000);

[p,f] = spectrum(y, nfft, noverlap);
f2 = 0.5*f;

%  AR1 process 
  rho = corr(y, y, 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 = find((f.*rn')==max(f.*rn'));

%  AR2 process
  rho1 = corr(y, y, 1);
  rho2 = corr(y, y, 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);
  rn2 = rn2 * mean(p(:,1));
  maxrn2 = find((f.*rn2')==max(f.*rn2'));

%  Calculate theoretical damping and such
dampfact = sqrt(-phi2);
fnot = acos(phi1/(2*sqrt(-phi2)));
phinot = (1-dampfact^2)/(1+dampfact^2)*tan(fnot);

damptime = -1/log(dampfact);
period = 2*pi/fnot;

%  Plot the data
figure(2); fo(1); clf;

subplot(4,2,6); 
  pos = get(gca, 'Position'); pos(1) = pos(1)+0.3*pos(3);
  pos(3) = 0.7*pos(3); set(gca, 'Position', pos);
  set(gca, 'Visible', 'off');
  axis([0 0.5 0 1]);
  text(0.1,0.65,'AR1');
  text(0.1,0.45,['Efold: ' ...
		 num2str(round(100./(2*pi*f2(maxrn)))/100) ...
		    ' yrs']);

subplot(4,2,5); cla;
  pos = get(gca, 'Position');
  pos(3) = 1.5*pos(3); set(gca, 'Position', pos);
  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');
  axis([-log(1000) log(6) -0.1 5.1]);
  set(gca, 'XTick', -log(2.^[10:-1:-6]), ...
           'XTickLabel', 2.^[10:-1:-6], 'YTick', [0:1:5]);
  grid on
  set(gca, 'fontsize', 8);

%  Plot spectra
subplot(4,2,8); 
  pos = get(gca, 'Position'); pos(1) = pos(1)+0.3*pos(3);
  pos(3) = 0.7*pos(3); set(gca, 'Position', pos);
  set(gca, 'Visible', 'off');
  axis([0 0.5 0 1]);
  text(0.1,0.65,'AR2');
  text(0.1,0.45,['Period: ' ...
		 num2str(round(100./(f2(maxrn2)))/100) ...
		    ' yrs']);

subplot(4,2,7); cla;
  pos = get(gca, 'Position');
  pos(3) = 1.5*pos(3); set(gca, 'Position', pos);
  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');
  axis([-log(1000) log(6) -0.1 5.1]);
  set(gca, 'XTick', -log(2.^[10:-1:-6]), ...
           'XTickLabel', 2.^[10:-1:-6], 'YTick', [0:1:5]);
  grid on
  set(gca, 'fontsize', 8);