%Similar to red_noise_prob but now also plotting using logarithmic axes.

clear;

N=1500; 
t=1:N;
alpha=0.7;
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; 
subplot(311); hold on; 
[s P]=fftPH(x,1,[]);
P=P*N;
s=[0 s 1/(2*diff(t(1:2)))];
P=[(sum(x))^2/N P (sum(Nyq'.*x))^2/N];
h=plot(s,2*repmat(chi2confPH(0.95,2),(N/2+1),1),'k--'); set(h,'linewidth',2);
plot(s,P,'r'); 
h=ylabel('squared Fourier coefficients'); font(h,16); 
h=xlabel('frequency'); font(h,16); 
font(gca,16); 

subplot(312); hold on; 
plot(s,P,'r');  logPH; axis tight; h4=axis;
ci=chi2confPH(0.95,2);
h=plot([1/1000 1/1000],1*ci,'k'); set(h,'linewidth',2);
h=plot([1/1000 1/1000],1,'k.'); set(h,'markersize',25);
h=ylabel('squared Fourier coefficients'); font(h,16); 
h=xlabel('frequency'); font(h,16); 
font(gca,16); 
axis([s(1) h4(2:end)]); 

subplot(313); hold on; 
P2=[]; s2=[]; 
for ct=2:2:length(P)-1; 
  P2(end+1)=mean(P(ct:ct+1)); 
  s2(end+1)=mean(s(ct:ct+1)); 
end;
P=P2;
s=s2;
plot(s,P,'r');  logPH;
ci=chi2confPH(0.95,4);
h=plot([1/700 1/700],1*ci,'k'); set(h,'linewidth',2);
h=plot([1/700 1/700],1,'k.'); set(h,'markersize',25);
h=ylabel('averaged coefficients ( \nu=4) '); font(h,16); 
h=xlabel('frequency'); font(h,16); 
font(gca,16); 
axis([h4(1:end)]);