Mercurial > hg > ltpda
comparison m-toolbox/test/test_ao_psd_variance_montecarlo.m @ 0:f0afece42f48
Import.
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Wed, 23 Nov 2011 19:22:13 +0100 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:f0afece42f48 |
|---|---|
| 1 % test_ao_psd_variance_montecarlo | |
| 2 % | |
| 3 % Tests that the standard deviation returned by ao.dy in one | |
| 4 % frequency bin is equivalent to the matlab's std taking | |
| 5 % considering all realisations | |
| 6 % | |
| 7 % M Nofrarias 22-07-09 | |
| 8 % | |
| 9 % $Id: test_ao_psd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $ | |
| 10 | |
| 11 % function test_ao_psd_variance_montecarlo() | |
| 12 | |
| 13 clear | |
| 14 | |
| 15 % data | |
| 16 nsecs = 500; | |
| 17 fs = 5; | |
| 18 pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))'); | |
| 19 | |
| 20 % Window | |
| 21 Nfft = 100; | |
| 22 win = specwin('Hanning', Nfft); | |
| 23 pl2 = plist('Nfft',Nfft, 'win',win,'Olap',-1,'scale','PSD') | |
| 24 | |
| 25 % loop | |
| 26 for i = 1:100 | |
| 27 a(i) = ao(pl); | |
| 28 b1(i) = psd(a(i),pl2); | |
| 29 % matlab's | |
| 30 [txy, f] = pwelch(a(i).data.y, win.win, Nfft/2, Nfft, a(i).data.fs); | |
| 31 b2(i) = ao(fsdata(f.', txy.')); | |
| 32 end | |
| 33 | |
| 34 %% mean | |
| 35 index = 6; | |
| 36 | |
| 37 % compare mean | |
| 38 mn = [mean(b1(:).y(index)) mean(b2(:).y(index))] | |
| 39 % error | |
| 40 err = std(b1(:).y(index)) | |
| 41 % compare standard deviation | |
| 42 clear rel | |
| 43 for i =1:len(b1(1)) | |
| 44 mn(i) = [mean(b1(:).y(i))]; % both means are equal | |
| 45 rel(:,i) = [std(b1(:).y(i)) mean(b1(:).dy(i))]/abs(mn(i)); | |
| 46 end | |
| 47 | |
| 48 figure | |
| 49 loglog(b1(1).x,rel') | |
| 50 figure | |
| 51 loglog(b1(1).x,rel(2,:)-rel(1,:)) | |
| 52 ylabel('difference (%)') |