%Uses pmtmPH to estimate the spectra of a red noise plus periodic component process

clear;

N=1000; 
t=1:N;
alpha=0.9;
x=0;
for ct=1:N;
  x(ct+1)=alpha*x(ct)+randn;
  Nyq(ct)=(-1)^ct; 
end;
x=x(1:end-1)'+0.3*sin(2*pi/10*t)';

figure(1); clf; 

for ct=[1:5 10 15 20 25]; 
  cla; hold on;
  pmtmPH(x,1,ct,1); 
  h=title(['multitaper spectral estimate with  \nu=',num2str(ct*2),' degrees of freedom']); 
  font(h,16);
  if ct==1; h4=axis; else; axis(h4); end;
  pause;
end;