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