%Comparing the periodogrm with a power spectral estimate formed by averaging
%neighboring periodogram estimates.  Also shows how the chi-squared PDF
%changes going from two to four degrees of freedom.

clear;

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

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


[s P]=fftPH(x,1,[]);
s=[0 s 1/(2*diff(t(1:2)))];
P=N*[(sum(x))^2/N^2 P (sum(Nyq'.*x))^2/N^2]; 


figure(1); clf; 
subplot(2,3,[1 2]); hold on; 
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 ( \nu=2) '); font(h,16); 
h=xlabel('frequency'); font(h,16); 
font(gca,16); 

subplot(233); 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); 

subplot(2,3,[4 5]); hold on; 
P2=[]; s2=[]; 
for ct=1: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;
N=N/2;
h=plot(s,2*repmat(chi2confPH(0.95,4),(N/2),1),'k--'); set(h,'linewidth',2);
%plot(s,2*ones(N/2,1),'k--'); 
plot(s2,P2,'r'); 
h=ylabel('averaged coefficients ( \nu=4) '); font(h,16); 
h=xlabel('frequency'); font(h,16); 
font(gca,16); 

subplot(236); hold on;
vals=0:1:20; 
[occ vals]=hist(P2,20);
bar(vals,2*occ/(N*diff(vals(1:2))),'b'); 
h=ylabel('normalized occurences'); font(h,16); 
h=xlabel('power density values'); font(h,16); 
font(gca,16);
%h=plot(vals,chi2pdf(vals,4),'r'); 
h=plot(vals,gampdf(vals,4/2,1),'r'); 
set(h,'linewidth',3); 

[2 mean(P2)],