Mercurial > hg > ltpda

comparison m-toolbox/test/test_matrix_MultiChannelNoise_2.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST matrix/MultiChannelNoise
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test_matrix_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $';
11
12 %% Loading spectra
13
14 fs = 10;
15 f = [logspace(-6,log10(fs/2),2000)]';
16
17 % currPath = cd;
18 % cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
19 % [TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
20 % cd(currPath)
21
22 % get noise psd shapes
23 % parameters names
24 parnames = {'DH','S11','S1D','SD1','SDD','Dh1','Dh2','w1','w2',...
25 'TH','Th1','Th2'};
26
27 % nominal params values
28 parvalues = {1,1,0,-1e-4,1,1,1,-1.3e-6,-2.0e-6,0.35,0.25,0.28};
29
30 % Noise model
31 Ns = matrix(plist('built-in','mdc3_ifo2ifo_noise'));
32
33 % evalueat model
34 [rw,cl]=size(Ns.objs);
35 for ii=1:rw
36 for jj=1:cl
37 Ns.objs(ii,jj).setParams(parnames,parvalues);
38 Ns.objs(ii,jj).setXvals(f);
39 CSD(ii,jj)=eval(Ns.objs(ii,jj));
40 if ii==jj
41 CSD(ii,jj)= abs(CSD(ii,jj));
42 end
43 end
44 end
45
46 %% Plot
47
48 plpl = plist('Linecolors', {'k', 'r', 'b'},...
49 'LineStyles', {'-', '--', ':'},...
50 'LineWidths', {3, 3, 3},...
51 'Legends', {'S_{o_1 o_1}', 'S_{o_1 o_\Delta}', 'S_{o_\Delta o_\Delta}'});
52
53 % iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))
54 iplot(CSD(1,1),CSD(1,2),CSD(2,2),plpl)
55
56 %%
57 %--------------------------------------------------------------------------
58 % 2 Dimensional Noisegen Filter
59 %--------------------------------------------------------------------------
60
61
62 %% 2 dim noisegen filter in z domain - symbolic precision All pass
63
64 % plist for noise generation
65 plns = plist(...
66 'csd', CSD, ...
67 'targetobj', 'miir', ...
68 'fs', fs, ...
69 'MaxIter', 70, ...
70 'PoleType', 3, ...
71 'MinOrder', 25, ...
72 'MaxOrder', 55, ...
73 'Weights', 2, ...
74 'Plot', false,...
75 'MSEVARTOL', 1e-1,...
76 'FITTOL', 1e-4,...
77 'UseSym', 'on',...
78 'iunits', 'm',...
79 'ounits', 'm');
80
81 na = matrix(plns);
82
83 %% Check response of CSD
84
85 na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
86 na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
87 na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
88 na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
89
90 mna = matrix([na11 na12;na21 na22]);
91
92 ECSD = mna*conj(transpose(mna));
93
94 for kk = 1:numel(CSD)
95 CSD(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
96 CSD(kk).setXunits(unit('Hz'));
97 ECSD.objs(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
98 ECSD.objs(kk).setXunits(unit('Hz'));
99 end
100
101 % calculate cross-coherence
102 COH = CSD(1,2)./sqrt(CSD(1,1).*CSD(2,2));
103 COH.simplifyYunits;
104 RCOH = real(COH);
105 ICOH = imag(COH);
106
107 ECOH = ECSD.objs(1,2)./sqrt(ECSD.objs(1,1).*ECSD.objs(2,2));
108 ECOH.simplifyYunits;
109 RECOH = real(ECOH);
110 IECOH = imag(ECOH);
111
112 % set plist for plotting
113
114 plpl = plist('Linecolors', {'k', 'r', 'b'},...
115 'LineStyles', {'-', '--', ':'},...
116 'LineWidths', {3, 3, 3},...
117 'Legends', {'Model', 'Fit', 'Residuals'});
118
119 iplot(CSD(1,1),ECSD.objs(1,1),abs(CSD(1,1)-ECSD.objs(1,1)),plpl)
120 iplot(CSD(2,2),ECSD.objs(2,2),abs(CSD(2,2)-ECSD.objs(2,2)),plpl)
121
122 plpl = plist('Linecolors', {'k', 'r', 'b'},...
123 'LineStyles', {'-', '--', ':'},...
124 'LineWidths', {3, 3, 3},...
125 'Legends', {'Model', 'Fit', 'Residuals'},...
126 'Yscales', {'All', 'lin'},...
127 'Ylabels', {'All', 'Re(coh)'});
128
129 iplot(RCOH,RECOH,abs(RCOH-RECOH),plpl)
130
131 plpl = plist('Linecolors', {'k', 'r', 'b'},...
132 'LineStyles', {'-', '--', ':'},...
133 'LineWidths', {3, 3, 3},...
134 'Legends', {'Model', 'Fit', 'Residuals'},...
135 'Yscales', {'All', 'lin'},...
136 'Ylabels', {'All', 'Im(coh)'});
137
138 iplot(ICOH,IECOH,abs(ICOH-IECOH),plpl)
139
140 %% Check response of Filters
141
142 % na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
143 % na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
144 % na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
145 % na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
146
147
148 % iplot(TF(1,1),na11,abs(TF(1,1) - na11))
149 % iplot(TF(1,2),na12,abs(TF(1,2) - na12))
150 % iplot(TF(2,1),na21,abs(TF(2,1) - na21))
151 % iplot(TF(2,2),na22,abs(TF(2,2) - na22))
152
153 iplot(na11,na12,na21,na22)
154
155 %%
156 %--------------------------------------------------------------------------
157 % 2 Dimensional Noise Generation
158 %--------------------------------------------------------------------------
159 nsecs = 3e5;
160 plng = plist(...
161 'model', na,...
162 'Nsecs', nsecs,...
163 'fs', fs,...
164 'yunits', 'm');
165
166 no = matrix(plng);
167
168 %% Check psd
169
170 no1 = no.objs(1);
171 no2 = no.objs(2);
172
173 no1xx = no1.psd;
174 no2xx = no2.psd;
175
176 iplot(no1xx,csd11,no2xx,csd22)
177
178 %% Check cross-coherence
179
180 no1 = no.objs(1);
181 no2 = no.objs(2);
182 nfft = 6e4;
183
184 no12xx = cpsd(no1,no2,plist('nfft',nfft));
185 no1xx = no1.psd(plist('nfft',nfft));
186 no2xx = no2.psd(plist('nfft',nfft));
187 %%
188
189 ecoh = no12xx./(sqrt(no1xx.*no2xx));
190 recoh = real(ecoh);
191 iecoh = imag(ecoh);
192
193 iplot(recoh,RECOH,plist('Yscales',{'All','lin'}))
194 iplot(iecoh,IECOH,plist('Yscales',{'All','lin'}))
195
196 %% Get averaged noise spectral curves
197
198 N = 500;
199 nsecs = 3e5;
200 nfft = 6e4;
201
202 Mnxx1 = 0; Mn2xx1 = 0;
203 Mnxx2 = 0; Mn2xx2 = 0;
204
205 Mnxy1 = 0; Mn2xy1 = 0;
206 Mnxy2 = 0; Mn2xy2 = 0;
207
208 tic
209 for ii = 1:N
210
211 fprintf('===== iter %s =====\n',num2str(ii))
212
213 % noise generation
214 plng = plist(...
215 'model', na,...
216 'Nsecs', nsecs,...
217 'fs', fs,...
218 'yunits', 'm');
219
220 no = matrix(plng);
221
222 % periodogram calculation
223 no1 = no.objs(1);
224 no2 = no.objs(2);
225
226 no1xx = no1.psd;
227 no2xx = no2.psd;
228 f1 = no1xx.x;
229
230 % running average and variance for periodogram
231 if ii == 1
232 Mnxx1 = no1xx.y;
233 Mnxx2 = no2xx.y;
234 else
235 % Qxx1 = no1xx.y - Mnxx1;
236 Mnxx1 = Mnxx1 + (no1xx.y - Mnxx1)./ii;
237 Mn2xx1 = Mn2xx1 + (no1xx.y - Mnxx1).*(no1xx.y - Mnxx1);
238
239 % Qxx2 = no2xx.y - Mnxx2;
240 Mnxx2 = Mnxx2 + (no2xx.y - Mnxx2)./ii;
241 Mn2xx2 = Mn2xx2 + (no2xx.y - Mnxx2).*(no2xx.y - Mnxx2);
242 end
243 clear no1xx no2xx
244
245
246 % cross coherence calculation
247 no12xx = cpsd(no1,no2,plist('nfft',nfft));
248 no1xx2 = no1.psd(plist('nfft',nfft));
249 no2xx2 = no2.psd(plist('nfft',nfft));
250 f2 = no1xx2.x;
251
252 ecoh = (no12xx.y)./(sqrt((no1xx2.y).*(no2xx2.y)));
253 recoh = real(ecoh);
254 iecoh = imag(ecoh);
255
256 clear no12xx no1xx2 no2xx2
257 % running average and variance for cross coherence
258 if ii == 1
259 Mnxy1 = recoh;
260 Mnxy2 = iecoh;
261 else
262 % Qxy1 = recoh - Mnxy1;
263 Mnxy1 = Mnxy1 + (recoh - Mnxy1)./ii;
264 Mn2xy1 = Mn2xy1 + (recoh - Mnxy1).*(recoh - Mnxy1);
265
266 % Qxy2 = iecoh - Mnxy2;
267 Mnxy2 = Mnxy2 + (iecoh - Mnxy2)./ii;
268 Mn2xy2 = Mn2xy2 + (iecoh - Mnxy2).*(iecoh - Mnxy2);
269 end
270
271 end
272 toc
273
274 % get mean and variance for psd 1
275 Sxx1 = Mnxx1; % mean
276 if N == 1
277 Svxx1 = Sxx1;
278 else
279 Svxx1 = (Mn2xx1/(N-1))/N; % variance
280 end
281 % get mean and variance for psd 2
282 Sxx2 = Mnxx2; % mean
283 if N == 1
284 Svxx2 = Sxx2;
285 else
286 Svxx2 = (Mn2xx2/(N-1))/N; % variance
287 end
288 % get mean and variance for real part of cross coherence
289 Sxy1 = Mnxy1; % mean
290 if N == 1
291 Svxy1 = Sxy1;
292 else
293 Svxy1 = (Mn2xy1/(N-1))/N; % variance
294 end
295 % get mean and variance for imaginary part of cross coherence
296 Sxy2 = Mnxy2; % mean
297 if N == 1
298 Svxy2 = Sxy2;
299 else
300 Svxy2 = (Mn2xy2/(N-1))/N; % variance
301 end
302
303 % save results
304
305 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat','Sxx1')
306 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat','Svxx1')
307 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat','Sxx2')
308 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat','Svxx2')
309 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy1.mat','Sxy1')
310 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy1.mat','Svxy1')
311 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy2.mat','Sxy2')
312 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy2.mat','Svxy2')
313 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f1.mat','f1')
314 save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f2.mat','f2')
315 save(na,'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\na.mat')
316
317 %% Load saved data
318
319 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat'
320 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat'
321 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat'
322 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat'
323 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy1.mat'
324 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy1.mat'
325 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy2.mat'
326 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy2.mat'
327 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f1.mat'
328 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f2.mat'
329
330 %% Plot spectra
331
332 figure()
333
334 p22 = loglog(f1(1:end-1),Sxx1(1:end-1) + sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
335 hold on
336 grid on
337 p23 = loglog(f1(1:end-1),Sxx1(1:end-1) - sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
338 p21 = loglog(f1(1:end-1),Sxx1(1:end-1),'r','LineWidth',3);
339
340 p11 = loglog(ECSD.objs(1,1).x,ECSD.objs(1,1).y,'k','LineWidth',3);
341
342
343 xlabel('Frequency [Hz]','fontsize',18)
344 ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
345 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
346
347 figure()
348
349 p22 = loglog(f1(1:end-1),Sxx2(1:end-1) + sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
350 hold on
351 grid on
352 p23 = loglog(f1(1:end-1),Sxx2(1:end-1) - sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
353 p21 = loglog(f1(1:end-1),Sxx2(1:end-1),'r','LineWidth',3);
354
355 p11 = loglog(ECSD.objs(2,2).x,ECSD.objs(2,2).y,'k','LineWidth',3);
356
357
358 xlabel('Frequency [Hz]','fontsize',18)
359 ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
360 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
361
362
363 %% Plot coherence
364
365 figure()
366
367 p22 = semilogx(f2,Sxy1 + sqrt(Svxy1),'b','LineWidth',3);
368 hold on
369 grid on
370 p23 = semilogx(f2,Sxy1 - sqrt(Svxy1),'b','LineWidth',3);
371 p21 = semilogx(f2,Sxy1,'r','LineWidth',3);
372
373 p11 = semilogx(RECOH.x,RECOH.y,'k','LineWidth',3);
374
375
376 xlabel('Frequency [Hz]','fontsize',18)
377 ylabel('Re(coh) []','fontsize',18)
378 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
379
380 figure()
381
382 p22 = semilogx(f2,Sxy2 + sqrt(Svxy2),'b','LineWidth',3);
383 hold on
384 grid on
385 p23 = semilogx(f2,Sxy2 - sqrt(Svxy2),'b','LineWidth',3);
386 p21 = semilogx(f2,Sxy2,'r','LineWidth',3);
387
388 p11 = semilogx(IECOH.x,IECOH.y,'k','LineWidth',3);
389
390
391 xlabel('Frequency [Hz]','fontsize',18)
392 ylabel('Im(coh) []','fontsize',18)
393 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
394
395 %% Check response of CSD
396
397 currPath = cd;
398 cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
399 [TF,CSD] = mdc1_tf_models(plist('f',f1,'fs',fs));
400 cd(currPath)
401
402 na11 = resp(na.objs(1,1).filters,plist('f',f1,'bank','parallel'));
403 na12 = resp(na.objs(1,2).filters,plist('f',f1,'bank','parallel'));
404 na21 = resp(na.objs(2,1).filters,plist('f',f1,'bank','parallel'));
405 na22 = resp(na.objs(2,2).filters,plist('f',f1,'bank','parallel'));
406
407 mna = matrix([na11 na12;na21 na22]);
408
409 ECSD = mna*conj(transpose(mna));
410
411 for kk = 1:numel(CSD)
412 CSD(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
413 CSD(kk).setXunits(unit('Hz'));
414 ECSD.objs(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
415 ECSD.objs(kk).setXunits(unit('Hz'));
416 end
417
418 plpl = plist('Linecolors', {'k', 'r', 'b'},...
419 'LineStyles', {'-', '--', ':'},...
420 'LineWidths', {3, 3, 3},...
421 'Legends', {'Model', 'Fit', 'Residuals'});
422
423 iplot(CSD(1,1),ECSD.objs(1,1),abs(CSD(1,1)-ECSD.objs(1,1)),plpl)
424 iplot(CSD(2,2),ECSD.objs(2,2),abs(CSD(2,2)-ECSD.objs(2,2)),plpl)
425
426 %% get expected response + window
427
428 % we should considere that the window has a spectral response that can add
429 % some systematics inside our spectral estimations. Since the window is a
430 % product in time domain with the time series, that it corresponds to a
431 % convolution in frequency domain. We can ifft the frequency response of
432 % the filters, multiply in time domain for the window and then come back to
433 % freq domain by fft.
434 na = matrix('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\na.mat');
435 %%
436 % filter responses
437 fs = 10;
438 nsecs = 3e5;
439 nfft = nsecs.*fs;
440 ff1 = (0:floor(nfft/2)).*fs./nfft;
441
442 rsp11 = resp(na.objs(1,1).filters,plist('f',ff1,'bank','parallel'));
443 rsp11 = rsp11.y;
444 % rsp11 = [rsp11(1:end-1);flipud(conj(rsp11(2:end)))];
445 rsp12 = resp(na.objs(1,2).filters,plist('f',ff1,'bank','parallel'));
446 rsp12 = rsp12.y;
447 % rsp12 = [rsp12(1:end-1);flipud(conj(rsp12(2:end)))];
448 rsp21 = resp(na.objs(2,1).filters,plist('f',ff1,'bank','parallel'));
449 rsp21 = rsp21.y;
450 % rsp21 = [rsp21(1:end-1);flipud(conj(rsp21(2:end)))];
451 rsp22 = resp(na.objs(2,2).filters,plist('f',ff1,'bank','parallel'));
452 rsp22 = rsp22.y;
453 % rsp22 = [rsp22(1:end-1);flipud(conj(rsp22(2:end)))];
454 clear na
455
456 % psd calc (npsd should contain nfft samples)
457 npsd11 = rsp11.*conj(rsp11) + rsp12.*conj(rsp12);
458 % npsd11 = [npsd11;flipud(npsd11(2:end-1))];
459 npsd22 = rsp21.*conj(rsp21) + rsp22.*conj(rsp22);
460 % npsd22 = [npsd22;flipud(npsd22(2:end-1))];
461
462 %%
463
464 % get window frequency response
465
466 % BH92 in frequency domain
467 W = BH92(2.*pi.*ff1./fs,nfft);
468 W = W./max(abs(W));
469 W = W.*conj(W);
470
471 %% ifft dello spettro
472
473 % npsd11 = [npsd11;flipud(npsd11(2:end-1))]./2;
474 % npsd22 = [npsd22;flipud(npsd22(2:end-1))]./2;
475
476 % g11 = ifft(npsd11);
477 % g22 = ifft(npsd22);
478
479 W = blackmanharris(nfft);
480 if size(W,2)>1
481 W = W.';
482 end
483 W = W.^2;
484 W = ifftshift(W);
485
486 WS11 = fft(W.*g11);
487 WS22 = fft(W.*g22);
488
489 WS11 = 2.*WS11(1:numel(ff1));
490 WS22 = 2.*WS22(1:numel(ff1));
491
492 %% get convolution
493
494 % WS11 = conv(W,npsd11);
495 % WS11 = WS11(1:numel(ff1));
496
497
498 %% Plot spectra
499
500 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat'
501 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat'
502 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat'
503 load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat'
504 %%
505 figure()
506
507 p22 = loglog(ff1(1:end-1),Sxx1(1:end-1) + sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
508 hold on
509 grid on
510 p23 = loglog(ff1(1:end-1),Sxx1(1:end-1) - sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
511 p21 = loglog(ff1(1:end-1),Sxx1(1:end-1),'r','LineWidth',3);
512
513 p11 = loglog(ff1(1:end),WS11(1:end),'k','LineWidth',3);
514
515
516 xlabel('Frequency [Hz]','fontsize',18)
517 ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
518 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
519
520 %%
521 figure()
522
523 p22 = loglog(ff1(1:end-1),Sxx2(1:end-1) + sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
524 hold on
525 grid on
526 p23 = loglog(ff1(1:end-1),Sxx2(1:end-1) - sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
527 p21 = loglog(ff1(1:end-1),Sxx2(1:end-1),'r','LineWidth',3);
528
529 p11 = loglog(ff1(1:end-1),WS22(1:end-1),'k','LineWidth',3);
530
531
532 xlabel('Frequency [Hz]','fontsize',18)
533 ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
534 legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
535
536 %% get statistics
537
538 Sn11 = (npsd11(14:end-1)-Sxx1(14:end-1))./npsd11(14:end-1);
539 Sn11 = ao(cdata(Sn11));
540 Sn22 = (npsd22(14:end-1)-Sxx2(14:end-1))./npsd22(14:end-1);
541 Sn22 = ao(cdata(Sn22));
542
543 % get equivalent normal distribution and histogram
544 hpl1 = plist('N', 50);
545 hpl2 = plist('N', 200);
546
547 dS11h = hist(Sn11,hpl1);
548 dS11h = dS11h./max(dS11h);
549 dS11n = normdist(Sn11,hpl2);
550 dS11n = dS11n./max(dS11n);
551
552 dS22h = hist(Sn22,hpl1);
553 dS22h = dS22h./max(dS22h);
554 dS22n = normdist(Sn22,hpl2);
555 dS22n = dS22n./max(dS22n);
556
557 %% check
558
559 iplot(dS11h,dS11n)
560 iplot(dS22h,dS22n)
561
562 %% Plot results of the statistics analysis
563
564 figure()
565
566 p1 = bar(dS11h.x,dS11h.y,'b');
567 hold on
568 grid on
569 p2 = plot(dS11n.x,dS11n.y,'r','LineWidth',3);
570
571 xlabel('Normalized Deviation [\Deltax / x]','fontsize',18)
572 ylabel('Normalized Counts','fontsize',18)
573 legend([p1(1) p2(1)],'Histogram', 'Equivalent Norm. Distr.')
574
575 figure()
576
577 p1 = bar(dS22h.x,dS22h.y,'b');
578 hold on
579 grid on
580 p2 = plot(dS22n.x,dS22n.y,'r','LineWidth',3);
581
582 xlabel('Normalized Deviation [\Deltax / x]','fontsize',18)
583 ylabel('Normalized Counts','fontsize',18)
584 legend([p1(1) p2(1)],'Histogram', 'Equivalent Norm. Distr.')
585
586
587
588 %%
589 m1 = mean(Sn11.y);
590 m2 = mean(Sn22.y);
591 s1 = std(Sn11.y);
592 s2 = std(Sn22.y);
593
594
595
596
597
598 %%
599
600 nfft = 1e2;
601 fs = 10;
602 frq = linspace(0,pi,10000);
603 % D=DirichletKernel(frq,nfft);
604 % D=D./max(D);
605 % figure;
606 % plot(frq,20.*log10(D))
607 % grid on
608 W=BH92(frq,nfft);
609 % W=W./max(W);
610 figure;
611 plot(frq./pi,20.*log10(W))
612 grid on
613