--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_psd_variance_montecarlo.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,52 @@
+%  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 (%)')
\ No newline at end of file