--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_cohere_variance_montecarlo.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,58 @@
+%  test_ao_cohere_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_cohere_variance_montecarlo.m,v 1.1 2009/08/11 14:20:10 miquel Exp $
+
+% function test_ao_cohere_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 = 100;
+  win  = specwin('Hanning', Nfft);
+  pl2 = plist('Nfft',Nfft, 'win',win,'Olap',0)
+  
+  % loop
+  for  i = 1:100
+    a1 = ao(pl);
+    a2 = filter(a1,plist('filter', f1));
+    c1(i) = cohere(a1,a2,plist('Nfft',100,'type','MS'));
+    % Do with MATLAB
+    [cxy, f] = mscohere(a1.data.y, a2.data.y, win.win, Nfft/2, Nfft, a1.data.fs);
+    c2(i) = ao(fsdata(f, cxy));
+  end
+
+  %% mean
+  index = 6;
+  
+  % compare mean
+  mn = [mean(c1(:).y(index)) mean(c2(:).y(index))]
+  %  error
+  err = std(c1(:).y(index))
+  % compare standard deviation
+ clear rel
+  for i =1:len(c1(1))
+     mn(i) = [mean(c1(:).y(i))]; % both means are equal
+  rel(:,i) = [std(c1(:).y(i)) mean(c1(:).dy(i))]/abs(mn(i));
+  end
+
+   figure
+  loglog(c1(1).x,rel')
+  figure
+  loglog(c1(1).x(:),100*abs(rel(2,:)-rel(1,:)))
+ ylabel('difference (%)')
+