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