--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_lcpsd_variance_montecarlo.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,52 @@
+%  test_ao_ltfe_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_lcpsd_variance_montecarlo.m,v 1.1 2009/08/11 14:20:10 miquel Exp $
+
+% function test_ao_ltfe_variance_montecarlo()
+
+clear
+
+% data
+nsecs = 200;
+fs    = 10;
+pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'sin(2*pi*2*t) + randn(size(t))');
+
+
+% Make a filter
+f1 = miir(plist('type', 'highpass', 'fc', 4, 'fs', fs));
+
+% Window
+Nfft = 1000;
+win  = specwin('Hanning', Nfft);
+pl2 = plist('Lmin',200,'win',win,'Olap',-1)
+
+% loop
+for  i = 1:100
+  a1 = ao(pl);
+  a2 = filter(a1,plist('filter', f1));
+  t1(i) = lcpsd(a1,a2,pl2);
+end
+
+%% mean
+index = 6;
+
+clear rel
+for i =1:len(t1(1))
+  mn(i) = [mean(t1(:).y(i))]; % both means are equal
+  rel(:,i) = [std(t1(:).y(i)) mean(t1(:).dy(i))]/abs(mn(i));
+end
+
+k = find(t1(1).procinfo,'k')
+
+
+figure
+loglog(t1(1).x,rel')
+figure
+loglog(t1(1).x,abs(rel(2,:)-rel(1,:)))
+