+
−%  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 (%)')
+
−