0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % test_ao_lpsd_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 05-08-09
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % $Id: test_ao_lpsd_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_lpsd_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 %lpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 pl2 = plist('Kdes', 20, 'Jdes', 100);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 % loop
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 for i = 1:100
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 a(i) = ao(pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 b1(i) = lpsd(a(i),pl2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 %% compare
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 % select a realisation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 index = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 % compare estimated std to the one using all realisations
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 clear rel
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 for i =1:len(b1(1))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 mn(i) = [mean(b1(:).y(i))];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 rel(:,i) = [std(b1(:).y(i))/mn(i) b1(index).dy(i)/b1(index).y(i)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 % k = find(b1(1).procinfo,'k')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 loglog(b1(1).x,rel')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 hold
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 loglog(b1(1).x,abs(rel(2,:)-rel(1,:)),'-or')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 ylabel('difference (%)')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 % xlabel('Frequency [Hz]')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 % ylabel('Rel. error')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 % legend('montecarlo','estimated','difference')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 % print -deps2c test_ao_lpsd_variance_montecarlo_1.eps
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57
|