%  test_ao_cpsd_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_cpsd_variance_montecarlo.m,v 1.1 2009/08/11 14:20:11 miquel Exp $


 
  clear

  % data
  nsecs = 200;
  fs    = 10;
  pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
  
  
  % Make a filter
  f1 = miir(plist('type', 'highpass', 'fc', 4, 'fs', fs));
  
  % loop
  for  i = 1:100
    a1 = ao(pl);
    a2 = filter(a1,plist('filter', f1));
    b(i) = cpsd(a1,a2,plist('Nfft',100,'order',-1));
  end
  
 %% mean

 clear rel
  for i =1:len(b(1))
     mn(i) = [mean(b(:).y(i))]; % both means are equal
     rel(:,i) = [std(b(:).y(i)) mean(b(:).dy(i))]/abs(mn(i));
  end

  figure
  loglog(b(1).x,rel')
  figure
  loglog(b(1).x(:),abs(rel(2,:)-rel(1,:)))