clear;

N=1000; 
t=1:N;
x=randn(N,1)+0.15*sin(2*pi/10*t)'; 

for ct=1:N;
  Nyq(ct)=(-1)^ct; 
end;


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)],