%Generates a red noise time series and the associated periodogram.
%Illustrates why a consant error bar is not practical.  

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(2,3,[1 2]); hold on; 
plot(t,x,'r'); 
h=ylabel('white noise realizations'); font(h,16); 
h=xlabel('time'); font(h,16); 
font(gca,16); 
axis tight;

subplot(233); hold on;
[occ vals]=hist(x,20);
bar(vals,occ/(diff(vals(1:2))*N),'b'); 
h=ylabel('normalized occurences'); font(h,16); 
h=xlabel('time value'); font(h,16); 
font(gca,16); 
h=plot(vals,normpdf(vals,0,1),'r'); 
set(h,'linewidth',3); 

subplot(2,3,[4 5]); 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(236); hold on;
vals=0:1:20; 
[occ vals]=hist(P,20);
bar(vals,2*occ/(N*diff(vals(1:2))),'b'); 
h=ylabel('normalized occurences'); font(h,16); 
h=xlabel('periodogram values'); font(h,16); 
font(gca,16);
h=plot(vals,chi2pdf(vals,2),'r'); 
set(h,'linewidth',3); 

[2 mean(P)],