+
−%  test_ao_lpsd_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 05-08-09
+
−%
+
−% $Id: test_ao_lpsd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $
+
−
+
−% function test_ao_lpsd_variance_montecarlo()
+
−
+
−clear
+
−
+
−% data
+
−nsecs = 500;
+
−fs    = 5;
+
−pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
+
−
+
−%lpsd
+
−pl2 = plist('Kdes', 20, 'Jdes', 100);
+
−
+
−% loop
+
−for  i = 1:100
+
−  a(i) = ao(pl);
+
−  b1(i) = lpsd(a(i),pl2);
+
−end
+
−
+
−
+
−%% compare
+
−
+
−% select a realisation
+
−index = 10;
+
−
+
−% compare estimated std to the one using all realisations
+
−clear rel
+
−for i =1:len(b1(1))
+
−  mn(i) = [mean(b1(:).y(i))];
+
−  rel(:,i) = [std(b1(:).y(i))/mn(i) b1(index).dy(i)/b1(index).y(i)];
+
−end
+
−
+
−% k = find(b1(1).procinfo,'k')
+
−
+
−figure
+
−loglog(b1(1).x,rel')
+
−hold
+
−loglog(b1(1).x,abs(rel(2,:)-rel(1,:)),'-or')
+
−ylabel('difference (%)')
+
−% xlabel('Frequency [Hz]')
+
−% ylabel('Rel. error')
+
−% legend('montecarlo','estimated','difference')
+
−% print -deps2c test_ao_lpsd_variance_montecarlo_1.eps
+
−
+
−% end
+
−
+
−
+
−