% eda12_10
% illustration of variance of power spectral density
% case of random timeseries tapered with Hamming window

% create an uncorrelated random timeseries
N=2^10;
sd2true=64.0;
sdtrue=sqrt(sd2true);
d = random('Normal',0,sdtrue,N,1);

% generic time/frequency set up
Dt=0.025;
T = N*Dt;
t=Dt*[0:N-1]'; 
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]'; 
Nf=N/2+1;
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]'; 
Nw=Nf;

% Hamming window function
W = 0.54 - 0.46*cos(2*pi*[0:N-1]'/(N-1));
d = W.*d;

% window changes variance of signal by ff
ff = sum(W.*W)/N;

% true and estimated variance of d
sd2est=std(d)^2;
disp(sprintf('true data variance %f', ff*sd2true));
disp(sprintf('estimated data variance %f', sd2est));

% compute power spectral density
fftraw = fft(d);                       % raw DFT output
dbar = Dt*fftraw(1:Nf);                % Fourier transform, + freqs only
s2 = (2/T) * (abs(dbar).^2);           % power spectral density
f=f(1:Nf);
w=w(1:Nf);

% statistics of power spectral density
s2meanest = mean(s2);
s2varest=std(s2)^2;
s2sigmaest=sqrt(s2varest);
sd2psd=Nf*Df*s2meanest;
disp(sprintf('variance of data estimated from PSD %f', sd2psd));
disp(' ');


% true mean and variance
p=2; % degrees of freedom
c = ( ff*sd2true ) / (2*Nf*Df);  % normalization factor
s2meantrue = p*c;
s2vartrue = 2*p*c^2;
disp('power spectral density');
disp(sprintf('true mean %f variance %f',s2meantrue,s2vartrue));
disp(sprintf('estimated mean %f variance %f',s2meanest,s2varest));
disp(' ');

% plot data
figure(1)
clf;
set(gca,'LineWidth',2);
hold on;
plot(t,d,'k-','LineWidth',2);
plot([0, T]', [-2*sdtrue*sqrt(ff), -2*sdtrue*sqrt(ff)], 'k-', 'LineWidth', 2 );
plot([0, T]', [0, 0], 'k-', 'LineWidth', 1 );
plot([0, T]', [2*sdtrue*sqrt(ff), 2*sdtrue*sqrt(ff)], 'k-', 'LineWidth', 2 );

% plot PSD
figure(2)
clf;
set(gca,'LineWidth',2);
hold on;
plot(f,s2,'k-','LineWidth',2);

% plot mean on spectra
plot([0, fmax], [s2meantrue, s2meantrue],'k-','LineWidth',2);

% histogram of values
Nhist=100;
[s2hist, bins] = hist(s2,Nhist);

% plot histogram
figure(3)
clf;
set(gca,'LineWidth',2);
hold on;
axis( [min(bins), max(bins), 0, max(s2hist)] );
plot(bins,s2hist,'k-','LineWidth',2);
plot([s2meantrue, s2meantrue], [0, max(s2hist)],'k-','LineWidth',2);

% plot true chi-squared distribution
Db=(bins(2)-bins(1))/c;
s2histtrue=Nf*Db*chi2pdf(bins/c,p);
plot(bins,s2histtrue,'k-','LineWidth',2);

% plot 95% confidence level on spectra
figure(2);
pv=0.95;
cv = chi2inv(0.95,p);
disp(sprintf('%f confidence at %f',pv,cv*c));
percent=100 * sum(s2hist(find(bins<=cv*c))) / sum(s2hist);
disp(sprintf('%f percent of points below %f confidence',percent,pv));
disp(' ');
plot([f(1), f(Nf)], [cv*c, cv*c],'k-','LineWidth',2);
figure(3);
plot([cv*c, cv*c], [0, max(s2hist)],'k-','LineWidth',2);

% probability that PSD exceeds maximum observed value
speak = max(s2);
ppeak = chi2cdf(speak/c,p);
disp(sprintf('probability of PSD exceeding maximum observed value %f', 1-ppeak));
disp(sprintf('Nf probability of PSD exceeding maximum observed value %f', 1-(ppeak^Nf)));


