Documentation of Fig_spect

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

Help text

    ctann = rpcs(1:450,num);

Cross-Reference Information

Listing of script Fig_spect

clear
data
load seas_ct2.mat
back
loadpcs;
num = 1;
for i = 1:3;
  if i == 2;
    ctann = ct2(1:1800)';
    tit = ['Years 101 to 550'];
  elseif i == 3;
%    ctann = rpcs(451:900,num);
    ctann = ct2(1801:3600)';
    tit = ['Years 551 to 1000'];
  elseif i == 1;
%    ctann = rpcs(1:900,num);
    ctann = ct2';
    tit = ['Years 101 to 1000'];
  end
  nx = length(ctann);
  nfft = 4*128; noverlap = 0.5*nfft;
  n2 = nfft/2;

  [p, f] = spectrum(ctann, nfft, noverlap); f2 = 2*f;
%  rho = (corr(ctann, ctann, 1));
  rho = (corr(ctann, ctann, 1)+sqrt(corr(ctann,ctann,2)))/2;
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn);
  rn = rn * mean(p(:,1));
  if nfft == 4*128;
    if i ~= 1;
      rner1 = rn * 2.64; rner5 = rn * 2.01;
    else
      rner1 = rn * 2.05; rner5 = rn * 1.68;
    end
  else nfft == 256;
    if i ~= 1;
      rner1 = rn * 3.55; rner5 = rn * 2.50;
    else
      rner1 = rn * 2.64; rner5 = rn * 2.01;
    end
  end
  dofx = nx / n2; 
  figure(1);% figure_orient;
%  pos = [1.75/8.5 10/11-i*(8.5/33) 5.5/8.5 2.25/11];
  sptalk(4,1,i);
     h = semilogy(f2, rn, '-r', ...
	      f2, rner5, '--r', f2, rner1, '-.r', f2, p(:,1), '-b' )
     set(gca, 'YTick', [.1 1 10 100], 'XTick', [0:.05:.5]);
     set(gca, 'YTickLabel', [0.1 1 10])
     set(h(4), 'linewidth', 2);
%     axis([f2(3) 0.5 .03 15])
%     axis([f2(3) 0.5 .01 4])
     axis([0 0.5 .01 4])
     grid on
     ylabel(tit);
     if i == 1;
       title(['Power Spectra for PC' num2str(num) ' of 0-270m HC, ' tit]);
     elseif i == 3;
       xlabel('Frequency (yrs^-^1)');
     end
end

cd ~/Thesis/Talk
%print -dpsc2 PC1_Power_Spectrum.ps
%print -dpsc2 CT_Power_Spectrum.ps

%  Do the same for the CT index

clear

data
load seas_ct2.mat
back
loadpcs;

num = 1;
for i = 1:2;
  if i == 2;
    ctann = ct2(1:3600);
    tit = ['CT index'];
    nfft = 4*128;
  elseif i == 1;
    ctann = rpcs(:,1);
    tit = ['PC1'];
    nfft = 128;
  end
  nx = length(ctann);
  noverlap = 0.5*nfft;
  n2 = nfft/2;

  [p, f] = spectrum(ctann, nfft, noverlap); f2 = f*nfft/256;
  if i == 1;
    rho = corr(ctann, ctann, 1);
  elseif i == 2;
    rho = (corr(ctann, ctann, 1)+sqrt(corr(ctann,ctann,2)))/2;
  end
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn);
  rn = rn * mean(p(:,1));
  rner1 = rn * 2.05; rner5 = rn * 1.68;
  dofx = nx / n2; 
  figure(2); fo(1);
%  pos = [1.75/8.5 10/11-i*(8.5/33) 5.5/8.5 2.25/11];
%  subplot('position', pos)
   sptalk(3,1,i);
     h = semilogy(f2, rn, 'r-', ...
	      f2, rner5, 'r--', f2, rner1, 'r-.', f2, p(:,1), '-b' )
     set(gca, 'YTick', [.1 1 10 100], 'XTick', [0:.05:.5]);
     set(gca, 'YTickLabel', [0.1 1 10])
     set(h(4), 'linewidth', 2);
     if i == 1; axis([f2(3) 0.5 .03 15]);
     else; axis([f2(3) 0.5 .01 5]); end;
     grid on
     ylabel(tit);
     if i == 1;
       title(['Power Spectra;  Years 101-1000']);
     end
end

[b1, a1] = butter(9, 2/9);
[h1, w1] = freqz(b1, a1, 128); h1 = abs(h1).^2;
[b2, a2] = butter(9, 2/10);
[h2, w2] = freqz(b2, a2, 128); h2 = 1-abs(h2).^2;
f2 = w1./(2*pi);
[p, f] = spectrum(rpcs(:,1), 128, 64);

sptalk(3,1,3);
  h = plot(f2, h1, '-r', f2, h2, '-b');
  set(h, 'linewidth', 2);
  grid on;
  axis([f(3)/2 0.5 -0.1 1.1]);
  set(gca, 'YTick', -1:.2:1, 'XTick', [0:.05:.5]);
  ylabel('Filter Response');
  legend(h, 'LP9 Filter', 'HP10 Filter');

cd ~/Thesis/Talk
print -dpsc PC_CT_Spect_Filtfunc.ps