%  test_ao_psd_variance_montecarlo
%
% Tests that the standard deviation returned by ao.dy in one
% frequency bin is equivalent to the matlab's std taking
% considering all realisations
%
% M Nofrarias 22-07-09
%
% $Id: test_ao_psd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $

% function test_ao_psd_variance_montecarlo()

clear

% data
nsecs = 500;
fs    = 5;
pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');

% Window
Nfft = 100;
win  = specwin('Hanning', Nfft);
pl2 = plist('Nfft',Nfft, 'win',win,'Olap',-1,'scale','PSD')

% loop
for  i = 1:100
  a(i) = ao(pl);
  b1(i) = psd(a(i),pl2);
  % matlab's
  [txy, f] = pwelch(a(i).data.y,  win.win, Nfft/2, Nfft, a(i).data.fs);
  b2(i) = ao(fsdata(f.', txy.'));
end

%% mean
index = 6;

% compare mean
mn = [mean(b1(:).y(index)) mean(b2(:).y(index))]
%  error
err = std(b1(:).y(index))
% compare standard deviation
clear rel
for i =1:len(b1(1))
  mn(i) = [mean(b1(:).y(i))]; % both means are equal
  rel(:,i) = [std(b1(:).y(i)) mean(b1(:).dy(i))]/abs(mn(i));
end

figure
loglog(b1(1).x,rel')
figure
loglog(b1(1).x,rel(2,:)-rel(1,:))
ylabel('difference (%)')