0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % Test script for ao/confint
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % L Ferraioli 20-03-09
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % $Id: test_ao_confint.m,v 1.4 2011/04/29 14:03:08 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 %% params
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 fs = 1; % Hz
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 %% Set the model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 d = [1 -1.1 0.5];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 h11 = miir([1 -0.5],d,fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 h12 = miir([0 -0.5],d,fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 h21 = miir([0 0.4],d,fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 h22 = miir([1 -0.6],d,fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 %% Theoretical sectra
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 f = logspace(-5,log10(0.5),300);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 f = f.';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 plresp = plist('f',f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 rh11 = resp(h11,plresp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 rh12 = resp(h12,plresp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 rh21 = resp(h21,plresp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 rh22 = resp(h22,plresp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % psd11
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 G11.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 % psd 22
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 G22.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % cpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 G12.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 % mscohere
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 sK12 = (G12.y.*conj(G12.y))./(G11.y.*G22.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 sK12 = ao(plist('xvals', f, 'yvals', sK12, 'fs', fs, 'type', 'fsdata'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 sK12.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 %% Noise generation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 w1 = ao(plist('tsfcn','randn(size(t))','fs',1,'nsecs',1e5,'yunits','m'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 w2 = ao(plist('tsfcn','randn(size(t))','fs',1,'nsecs',1e5,'yunits','m'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 x11 = filter(w1,h11);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 x12 = filter(w2,h12);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 x21 = filter(w1,h21);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 x22 = filter(w2,h22);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 x1 = x11 + x12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 x2 = x21 + x22;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 %% mslcohere
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 sk12log = lcohere(x1,x2,plist('type','MS'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 if numel(sk12log)>1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 sk12log = sk12log(1,2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 % iplot(sK12,sk12log,plist('Yscales',{'All','linear'}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 out = confint(sk12log,plist('method','mslcohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 lc = out.getObjectAtIndex(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 uc = out.getObjectAtIndex(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 var = out.getObjectAtIndex(3);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 iplot(sK12,sk12log,lc,uc,plist('Yscales',{'All','linear'}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 dof = getdof(sk12log,plist('method','mslcohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 %% shaded plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 x = lc.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 y1 = lc.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 y2 = uc.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 mod = sK12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 hd = plot(mod.data.x, mod.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 hf = plot(sk12log.data.x,sk12log.data.y,'Color', 'g');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 set(gca,'xscale','log','yscale','lin');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 ylabel('Magnitude Squared Coherence');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 legend([hd hf ha(2)],{'Model MSC','Sample MSC','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 %% lpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 g11lpsd = lpsd(x1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 % iplot(2.*G11,g11lpsd)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 outlpsd = confint(g11lpsd,plist('method','lpsd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 lclpsd = outlpsd.getObjectAtIndex(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 uclpsd = outlpsd.getObjectAtIndex(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 varlpsd = outlpsd.getObjectAtIndex(3);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 iplot(2.*G11,g11lpsd,lclpsd,uclpsd)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 doflpsd = getdof(g11lpsd,plist('method','lpsd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 %% shaded plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 x = lclpsd.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 y1 = lclpsd.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 y2 = uclpsd.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 mod = 2.*G11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 hd = plot(mod.data.x, mod.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 hf = plot(g11lpsd.data.x,g11lpsd.data.y,'Color', 'g');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 set(gca,'xscale','log','yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 ylabel('Power Spectral Density [m^{2} / Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 legend([hd hf ha(2)],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 %% mscohere
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 sk12lin = cohere(x1,x2,plist('Nfft',1e3,'type','MS'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 if numel(sk12lin)>1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 sk12lin = sk12lin(1,2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160 iplot(sK12,sk12lin,plist('Yscales',{'All','linear'}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162 outlin = confint(sk12lin,plist('method','mscohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 lclin = outlin.getObjectAtIndex(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 uclin = outlin.getObjectAtIndex(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 varlin = outlin.getObjectAtIndex(3);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167 iplot(sK12,sk12lin,lclin,uclin,plist('Yscales',{'All','linear'}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169 doflin = getdof(sk12lin,plist('method','mscohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
171 %% shaded plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
172
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
173 x = lclin.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
174 y1 = lclin.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
175 y2 = uclin.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
176 mod = sK12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
177
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
178 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
179 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
180 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
181 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
182 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
183 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
184 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
185
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
186 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
187 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
188 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
189 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
190 hd = plot(mod.data.x, mod.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
191 hf = plot(sk12lin.data.x,sk12lin.data.y,'Color', 'g');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
192
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
193 set(gca,'xscale','log','yscale','lin');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
194 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
195 ylabel('Magnitude Squared Coherence');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
196 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
197 legend([hd hf ha(2)],{'Model MSC','Sample MSC','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
198
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
199 %% psd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
200
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
201 g11psd = psd(x1,plist('Nfft',1e4));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
202
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
203 iplot(2.*G11,g11psd)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
204
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
205 outpsd = confint(g11psd,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
206 lcpsd = outpsd.getObjectAtIndex(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
207 ucpsd = outpsd.getObjectAtIndex(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
208 varpsd = outpsd.getObjectAtIndex(3);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
209
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
210 iplot(2.*G11,g11psd,lcpsd,ucpsd)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
211
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
212 dofpsd = getdof(g11psd,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
213
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
214 %% shaded plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
215
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
216 x = lcpsd.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
217 y1 = lcpsd.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
218 y2 = ucpsd.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
219 mod = 2.*G11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
220
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
221 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
222 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
223 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
224 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
225 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
226 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
227 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
228
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
229 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
230 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
231 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
232 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
233 hd = plot(mod.data.x, mod.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
234 hf = plot(g11psd.data.x,g11psd.data.y,'Color', 'g');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
235
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
236 set(gca,'xscale','log','yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
237 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
238 ylabel('Power Spectral Density [m^{2} / Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
239 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
240 legend([hd hf ha(2)],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
241
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
242 %% Shaded plot 2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
243
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
244 x = [lcpsd.data.x; flipud(lcpsd.data.x)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
245 y = [ucpsd.data.y; flipud(lcpsd.data.y)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
246 mod = 2.*G11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
247
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
248 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
249 ha = fill(x,y,'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
250 set(ha,'EdgeColor','r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
251
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
252 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
253 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
254
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
255 hd = plot(mod.data.x, mod.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
256 hf = plot(g11psd.data.x,g11psd.data.y,'Color','g');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
257
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
258 set(gca,'xscale','log','yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
259 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
260 ylabel('Power Spectral Density [m^{2} / Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
261 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
262 legend([hd hf ha],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
263
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
264 %% Test it is working also with processed data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
265
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
266 g11psd = psd(x1,plist('Nfft',1e4));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
267
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
268 g11psd.setName('psd');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
269
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
270 plsp = plist('samples',[5 numel(g11psd.x)]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
271 g11psds = split(g11psd,plsp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
272
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
273 outpsd = confint(g11psds,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
274 lcpsd = outpsd.getObjectAtIndex(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
275 ucpsd = outpsd.getObjectAtIndex(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
276 varpsd = outpsd.getObjectAtIndex(3);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
277
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
278 iplot(2.*G11,g11psds,lcpsd,ucpsd)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
279
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
280 dofpsd = getdof(g11psd,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
281
|