--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_lpsd_variance_montecarlo.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,57 @@
+%  test_ao_lpsd_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 05-08-09
+%
+% $Id: test_ao_lpsd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $
+
+% function test_ao_lpsd_variance_montecarlo()
+
+clear
+
+% data
+nsecs = 500;
+fs    = 5;
+pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
+
+%lpsd
+pl2 = plist('Kdes', 20, 'Jdes', 100);
+
+% loop
+for  i = 1:100
+  a(i) = ao(pl);
+  b1(i) = lpsd(a(i),pl2);
+end
+
+
+%% compare
+
+% select a realisation
+index = 10;
+
+% compare estimated std to the one using all realisations
+clear rel
+for i =1:len(b1(1))
+  mn(i) = [mean(b1(:).y(i))];
+  rel(:,i) = [std(b1(:).y(i))/mn(i) b1(index).dy(i)/b1(index).y(i)];
+end
+
+% k = find(b1(1).procinfo,'k')
+
+figure
+loglog(b1(1).x,rel')
+hold
+loglog(b1(1).x,abs(rel(2,:)-rel(1,:)),'-or')
+ylabel('difference (%)')
+% xlabel('Frequency [Hz]')
+% ylabel('Rel. error')
+% legend('montecarlo','estimated','difference')
+% print -deps2c test_ao_lpsd_variance_montecarlo_1.eps
+
+% end
+
+
+