% eda06_13
% comparison of variances
% between time series and its power spectral density

clear all;

% random time series with variance sigmad^2
N=256;
Dt=0.01;
t=Dt*[0:N-1];
sigmad1=10;
sigmad1sq=sigmad1^2;
d = random('Normal',0,sigmad1,N,1);
disp(sprintf('true variance %f', sigmad1sq));

% variance estimated from time series
sigmad2sq = (d'*d)/N;
disp(sprintf('variance estimated from time series %f', sigmad2sq));

% set up for fitting of sines and cosines type fourier series
Nf=N/2+1;
fmax = 1/(2*Dt);
Df = fmax/(N/2);
f = Df*[0:Nf-1]';
Nw=Nf;
wmax = 2*pi*fmax;
Dw = wmax/(N/2);
w = Dw*[0:Nw-1]';

% set up G
M=N;
G=zeros(N,M);
G(:,1)=1;
for i = [1:M/2-1]
    j = 2*i;
    k = j+1;
    G(:,j)=cos(w(i+1).*t);
    G(:,k)=sin(w(i+1).*t);
end
G(:,M)=cos(w(Nw).*t);
    
% solve least squares problem
gtgi = 2*ones(M,1)/N;
gtgi(1)=1/N;
gtgi(M)=1/N;
mest = gtgi .* (G'*d);

sigmad3sq = 0.5*mest'*mest;
disp(sprintf('variance estimated from sine/cosines %f', sigmad3sq));

% standard set for discrete fourier transform
N=256;
Dt=0.1;
tmax=Dt*N; 
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;

% discrete fourier transform
C=fft(d);


% comparison of coefficients
disp('check agreement between sines/cosine and fft at 2nd frequency:');
disp(sprintf('A2 B2 %f %f  (2/N)real(C) -(2/N)imag(C) %f %f', mest(2), mest(3), (2/N)*real(C(2)), -(2/N)*imag(C(2)) ));

% variance estimated from fft coefficients
sigmad4sq = (2/(N^2)) * C(1:Nf)'*C(1:Nf);
disp(sprintf('variance estimated from fft coefficients %f', sigmad4sq));

% apprximation to integral transform
dbar=Dt*C;

% variance estimated from fourier transform
sigmad4sq = (2/(Dt^2*N^2)) * dbar(1:Nf)'*dbar(1:Nf);
disp(sprintf('variance estimated from fourier transform %f', sigmad4sq));

% integral from 0 to fny of |dbar(f)|2
theintegral = Df * dbar(1:Nf)'*dbar(1:Nf);
sigmad5sq = (2/(Df*Dt^2*N^2)) * theintegral;
disp(sprintf('variance estimated from frequency integral %f', sigmad5sq));

% definition of power spectral density
thepsd = (2/tmax) * abs(dbar(1:Nf)).^2;

% integral of power spectral density
P = Df * sum(thepsd);
sigmad6sq = P;
disp(sprintf('variance estimated from power spectral density %f', sigmad6sq));
