Mercurial > hg > ltpda

comparison m-toolbox/test/MDC2/test_mdc2_noisegen_accuracy.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 % A test script for to check the accuracy of the noise generation procedure
2 % for the MDC2
3 %
4 %
5 % L. Ferraioli 04-02-2009
6 %
7 % $Id: test_mdc2_noisegen_accuracy.m,v 1.4 2009/04/20 14:34:12 luigi Exp $
8 %
9
10 %% General use variables and vectors
11
12 fs = 1;
13 f = logspace(-6,log10(fs/2),300);
14 f = f.';
15 Nsecs = 1e6; % number of seconds
16
17 %% MDC2 Models
18
19 b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f));
20
21 %%
22
23 iplot(b)
24
25 %%
26
27 CSD = [b(1) b(2);conj(b(2)) b(3)];
28 MSC = b(2).*conj(b(2))./(b(1).*b(3));
29
30 %% some ploting
31
32 iplot(CSD(1,1),CSD(2,2),MSC)
33 iplot(CSD(1,1),CSD(2,2),CSD(1,2),CSD(2,1))
34
35 %% Make white noise
36
37 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
38 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
39
40 %% Noise generation
41
42 plng = plist(...
43 'csd11', CSD(1,1), ...
44 'csd12', CSD(1,2), ...
45 'csd21', CSD(2,1), ...
46 'csd22', CSD(2,2), ...
47 'MaxIter', 80, ...
48 'PoleType', 2, ...
49 'MinOrder', 35, ...
50 'MaxOrder', 40, ...
51 'Weights', 2, ...
52 'Plot', false,...
53 'MSEVARTOL', 1e-2,...
54 'FITTOL', 1e-5,...
55 'UseSym', 0,...
56 'Disp', false);
57
58 [ac1,ac2] = noisegen2D(a1,a2, plng);
59
60 % extracting filters
61 h11 = ac1.procinfo.find('Filt11');
62 h12 = ac1.procinfo.find('Filt12');
63 h21 = ac2.procinfo.find('Filt21');
64 h22 = ac2.procinfo.find('Filt22');
65
66
67 %% spectra calculation
68
69 acxx1 = lpsd(ac1);
70 acxx2 = lpsd(ac2);
71 acxx12 = lcohere(ac1,ac2);
72
73 % get dofs
74 dof1 = getdof(acxx1,plist('method','lpsd'));
75 dof2 = getdof(acxx2,plist('method','lpsd'));
76 dof12 = getdof(acxx12,plist('method','mslcohere'));
77
78 % Interpolate theorethical data
79 CSD11 = interp(CSD(1,1),plist('vertices',acxx1.data.x));
80 CSD22 = interp(CSD(2,2),plist('vertices',acxx2.data.x));
81 MSC = interp(MSC,plist('vertices',acxx12.data.x));
82
83 % getting variance
84 [lw1,up1,sig1] = confint(acxx1,plist('method','lpsd'));
85 [lw2,up2,sig2] = confint(acxx2,plist('method','lpsd'));
86 [lw12,up12,sig12] = confint(acxx12,plist('method','mslcohere'));
87
88 % chi square
89 k1 = acxx1 - CSD11;
90 k1 = k1.^2;
91 k1 = k1./sig1;
92 chi1 = sum(k1)./length(k1.data.x);
93
94 k2 = acxx2 - CSD22;
95 k2 = k2.^2;
96 k2 = k2./sig2;
97 chi2 = sum(k2)./length(k2.data.x);
98
99 k12 = acxx12 - MSC;
100 k12= k12.^2;
101 k12 = k12./sig12;
102 k12 = split(k12,plist('split_type','samples','samples',[3 342]));
103 chi12 = sum(k12)./length(k12.data.x);
104
105 proc1 = acxx1.procinfo;
106 proc2 = acxx2.procinfo;
107 proc12 = acxx12.procinfo;
108
109 %%
110
111 iplot(CSD11,acxx1,lw1,up1)
112 iplot(CSD22,acxx2,lw2,up2)
113 iplot(MSC,acxx12,lw12,up12)
114
115 %% Running the loop
116
117 N = 100;
118
119 for ii = 2:N
120 disp(num2str(ii))
121
122 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
123 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
124
125 x11 = filter(a1,h11);
126 x12 = filter(a2,h12);
127 x21 = filter(a1,h21);
128 x22 = filter(a2,h22);
129
130 x1 = x11 + x12;
131 x2 = x21 + x22;
132
133 % sectra
134 xx1 = lpsd(x1);
135 xx2 = lpsd(x2);
136 xx12 = lcohere(x1,x2);
137
138 % sum spectra
139 acxx1 = acxx1 + xx1;
140 acxx2 = acxx2 + xx2;
141 acxx12 = acxx12 + xx12;
142
143 % getting variance
144 [lw,up,sig1] = confint(xx1,plist('method','lpsd'));
145 [lw,up,sig2] = confint(xx2,plist('method','lpsd'));
146 [lw,up,sig12] = confint(xx12,plist('method','mslcohere'));
147
148 % chi square
149 k1 = xx1 - CSD11;
150 k1 = k1.^2;
151 k1 = k1./sig1;
152 xc1 = sum(k1)./length(k1.data.x);
153
154 k2 = xx2 - CSD22;
155 k2 = k2.^2;
156 k2 = k2./sig2;
157 xc2 = sum(k2)./length(k2.data.x);
158
159 k12 = xx12 - MSC;
160 k12= k12.^2;
161 k12 = k12./sig12;
162 k12 = split(k12,plist('split_type','samples','samples',[3 342]));
163 xc12 = sum(k12)./length(k12.data.x);
164
165 % adding chisquare
166 chi1 = chi1 + xc1;
167 chi2 = chi2 + xc2;
168 chi12 = chi12 + xc12;
169
170 end
171
172 %% Do Averages
173
174 % do average
175 acxx1 = acxx1./N;
176 acxx1.setName;
177 acxx2 = acxx2./N;
178 acxx2.setName;
179 acxx12 = acxx12./N;
180 acxx12.setName;
181
182 % average chi square
183 chi1 = chi1./N;
184 chi2 = chi2./N;
185 chi12 = chi12./N;
186
187 % dofs for the averages
188 adof1 = dof1.*N;
189 adof2 = dof2.*N;
190 adof12 = dof12.*N;
191
192 % get confidence bounds for the 'true' spactra
193 [lw1,up1,sig1] = confint(acxx1,plist('method','lpsd','dof',adof1));
194 [lw2,up2,sig2] = confint(acxx2,plist('method','lpsd','dof',adof2));
195 [lw12,up12,sig12] = confint(acxx12,plist('method','mslcohere','dof',adof12));
196
197 %% plotting
198
199 iplot(CSD(1,1),acxx1)
200 iplot(CSD(2,2),acxx2)
201 iplot(MSC,acxx12)
202
203 iplot(CSD(1,1),lw1,up1)
204 iplot(CSD(2,2),lw2,up2)
205 iplot(MSC,lw12,up12)
206
207 %% plotting conflevels 1
208
209 x = lw1.data.x;
210 y1 = lw1.data.y;
211 y2 = up1.data.y;
212
213 figure
214 y = [y1 (y2-y1)]; % y1 and y2 are columns
215 ha = area(x, y);
216 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
217 set(ha(2), 'FaceColor', 'r')
218 set(ha, 'LineStyle', 'none')
219 grid on
220
221 % plot the line edges
222 hold on
223 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
224 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
225 hd = plot(CSD(1,1).data.x, CSD(1,1).data.y);
226
227 set(gca,'xscale','log','yscale','log');
228 % set(gca,'Layer','top')
229 ylabel('Power Spectrum [m^{2} / Hz]');
230 xlabel('Frequency [Hz]');
231 legend([hd hb],{'Model CSD(1,1)','95% Conf. level'});
232
233 %% plotting conflevels 2
234
235 x = lw2.data.x;
236 y1 = lw2.data.y;
237 y2 = up2.data.y;
238
239 figure
240 y = [y1 (y2-y1)]; % y1 and y2 are columns
241 ha = area(x, y);
242 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
243 set(ha(2), 'FaceColor', 'r')
244 set(ha, 'LineStyle', 'none')
245 grid on
246
247 % plot the line edges
248 hold on
249 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
250 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
251 hd = plot(CSD(2,2).data.x, CSD(2,2).data.y);
252
253 set(gca,'xscale','log','yscale','log');
254 % set(gca,'Layer','top')
255 ylabel('Power Spectrum [m^{2} / Hz]');
256 xlabel('Frequency [Hz]');
257 legend([hd hb],{'Model CSD(2,2)','95% Conf. level'});
258
259 %% plotting conflevels 3
260
261 x = lw12.data.x;
262 y1 = lw12.data.y;
263 y2 = up12.data.y;
264
265 figure
266 y = [y1 (y2-y1)]; % y1 and y2 are columns
267 ha = area(x, y);
268 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
269 set(ha(2), 'FaceColor', 'r')
270 set(ha, 'LineStyle', 'none')
271 grid on
272
273 % plot the line edges
274 hold on
275 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
276 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
277 hd = plot(MSC.data.x, MSC.data.y);
278
279 set(gca,'xscale','log','yscale','log');
280 % set(gca,'Layer','top')
281 ylabel('Magnitude Square Coherence');
282 xlabel('Frequency [Hz]');
283 legend([hd hb],{'Model MSC','95% Conf. level'});