Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/dispersion.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 % DISPERSION computes the dispersion function
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: DISPERSION computes the dispersion function
5 %
6 % CALL: bs = dipersion(in,pl)
7 %
8 % INPUTS: in - matrix objects with input signals to the system
9 % model - symbolic models containing the transfer function model
10 %
11 % pl - parameter list
12 %
13 % OUTPUTS: bs - dispersion function AO
14 %
15 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'dispersion')">Parameters Description</a>
16 %
17 % VERSION: $Id: dispersion.m,v 1.1 2011/06/22 09:54:19 miquel Exp $
18 %
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21 function varargout = dispersion(varargin)
22
23 % Check if this is a call for parameters
24 if utils.helper.isinfocall(varargin{:})
25 varargout{1} = getInfo(varargin{3});
26 return
27 end
28
29 import utils.const.*
30 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
31
32 % Method can not be used as a modifier
33 if nargout == 0
34 error('### crb cannot be used as a modifier. Please give an output variable.');
35 end
36
37 % Collect input variable names
38 in_names = cell(size(varargin));
39 for ii = 1:nargin,in_names{ii} = inputname(ii);end
40
41 % Collect all AOs smodels and plists
42 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
43 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
44
45 % Combine plists
46 pl = parse(pl, getDefaultPlist);
47
48 % get params
49 params = find(pl,'FitParams');
50 numparams = find(pl,'paramsValues');
51 mmdl = find(pl,'model');
52 channel = find(pl,'channel');
53 mtxns = find(pl,'noise');
54 outModel = find(pl,'outModel');
55 bmdl = find(pl,'built-in');
56 f1 = find(pl,'f1');
57 f2 = find(pl,'f2');
58 pseudoinv = find(pl,'pinv');
59 tol = find(pl,'tol');
60 outNames = find(pl,'outNames');
61 inNames = find(pl,'inNames');
62
63 % Decide on a deep copy or a modify
64 fin = copy(mtxs, nargout);
65 n = copy(mtxns, nargout);
66 mdl = copy(mmdl,1);
67
68 % Get number of experiments
69 nexp = numel(fin);
70
71 % fft
72 % fin = fft(in);
73
74 % N should get before spliting, in order to convert correctly from psd to
75 % fft
76 N = length(fin(1).getObjectAtIndex(1).x);
77
78 % Get rid of fft f =0, reduce frequency range if needed
79 if ~isempty(f1) && ~isempty(f2)
80 fin = split(fin,plist('frequencies',[f1 f2]));
81 end
82
83 FMall = zeros(numel(params),numel(params));
84 % loop over experiments
85 for k = 1:nexp
86
87 utils.helper.msg(msg.IMPORTANT, sprintf('Analysis of experiment #%d',k), mfilename('class'), mfilename);
88
89 % if (((numel(n(1).objs)) == 2) && (numel(fin(1).objs) == 2))
90
91 % use signal fft to get frequency vector. Take into account signal
92 % could be empty or set to zero
93 % 1st channel
94 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
95 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
96 else
97 i1 = fin(k).getObjectAtIndex(1,1);
98 freqs = i1.x;
99 end
100 % 2nd channel
101 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
102 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
103 else
104 i2 = fin(k).getObjectAtIndex(2,1);
105 freqs = i2.x;
106 end
107
108 % Compute psd
109 n1 = lpsd(n(k).getObjectAtIndex(1,1));
110 n2 = lpsd(n(k).getObjectAtIndex(2,1));
111 n12 = lcpsd(n(k).getObjectAtIndex(1,1),n(k).getObjectAtIndex(2,1));
112
113 % interpolate to given frequencies
114 % noise
115 S11 = interp(n1,plist('vertices',freqs));
116 S12 = interp(n12,plist('vertices',freqs));
117 S22 = interp(n2,plist('vertices',freqs));
118 S21 = conj(S12);
119
120 % get some parameters used below
121 fs = S11.fs;
122
123
124 % if (isempty(outModel))
125
126 % if ~isempty(bmdl)
127 % compute built-in smodels
128 % for i = 1:4
129 % if strcmp(bmdl{i},'0');
130 % h(i) = smodel('0');
131 % h(i).setXvar('f');
132 % h(i).setXvals(freqs);
133 % h(i).setParams(params,numparams);
134 % else
135 % h(i) = smodel(plist('built-in',bmdl{i},'f',freqs));
136 % % set all params to all models. It is not true but harmless
137 % for ii = 1:numel(params)
138 % vecparams(ii) = {numparams(ii)*ones(size(freqs))};
139 % end
140 % h(i).setParams(params,vecparams);
141 % end
142 % end
143 % elseif ~isempty(mdl) && all(strcmp(class(mdl),'matrix'))
144
145 for i = 1:numel(mdl.objs)
146 % set Xvals
147 h(i) = mdl.getObjectAtIndex(i).setXvals(freqs);
148 % set alias
149 h(i).assignalias(mdl.objs(i),plist('xvals',freqs));
150 % set paramaters
151 h(i).setParams(params,numparams);
152 end
153 % differentiate and eval
154 for i = 1:length(params)
155 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
156 % differentiate symbolically
157 dH11 = diff(h(1),params{i});
158 dH12 = diff(h(3),params{i}); % taking into account matrix index convention h(2) > H(2,1)
159 dH21 = diff(h(2),params{i});
160 dH22 = diff(h(4),params{i});
161 % evaluate
162 d11(i) = eval(dH11,plist('output type','fsdata','output x',freqs));
163 d12(i) = eval(dH12,plist('output type','fsdata','output x',freqs));
164 d21(i) = eval(dH21,plist('output type','fsdata','output x',freqs));
165 d22(i) = eval(dH22,plist('output type','fsdata','output x',freqs));
166 end
167
168 % elseif ~isempty(mdl) && all(strcmp(class(mdl),'ssm'))
169 %
170 % meval = copy(mdl,1);
171 % % set parameter values
172 % meval.doSetParameters(params, numparams);
173 %
174 % % make numeric
175 % % meval.doSubsParameters(params, true);
176 %
177 % % get the differentiation step
178 % step = find(pl,'diffStep');
179 % if isempty(step)
180 % error('### Please input a step for the numerical differentiation')
181 % end
182 %
183 % % differentiate and eval
184 % for i = 1:length(params)
185 % utils.helper.msg(msg.IMPORTANT, sprintf('computing numerical differentiation with respect %s, Step:%4.2d ',params{i},step(i)), mfilename('class'), mfilename);
186 % % differentiate numerically
187 % dH = meval.parameterDiff(plist('names', params(i),'values',step(i)));
188 % % create plist with correct outNames (since parameterDiff change them)
189 % out1 = strrep(outNames{1},'.', sprintf('_DIFF_%s.',params{i})); % 2x2 case
190 % out2 =strrep(outNames{2},'.', sprintf('_DIFF_%s.',params{i}));
191 % spl = plist('set', 'for bode', ...
192 % 'outputs', {out1,out2}, ...
193 % 'inputs', inNames, ...
194 % 'reorganize', true,...
195 % 'f', freqs);
196 % % do bode
197 % d = bode(dH, spl);
198 % % assign according matlab's matrix notation: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
199 % d11(i) = d(1);
200 % d21(i) = d(2);
201 % d12(i) = d(3);
202 % d22(i) = d(4);
203 % end
204 %
205 % else
206 % error('### please introduce models for the transfer functions')
207 % end
208
209 % elseif (~isempty(outModel))
210 %
211 % % if(~isempty(outModel))
212 % for lll=1:size(outModel,1)
213 % for kkk=1:size(outModel,2)
214 % outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[f1 f2]));
215 % end
216 % end
217 % %end
218 %
219 % % Avoid numerical differentiation (faster for the magnetic case)
220 % Param{1} = [1 0;
221 % 0 0;];
222 % Param{2} = [0 0;
223 % 0 1;];
224 %
225 % for pp = 1:length(params)
226 % for ll = 1:size(outModel,1)
227 % for kk = 1:size(Param{pp},2)
228 % % index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
229 % tmp = 0;
230 % for innerIndex = 1:size(outModel,2)
231 % tmp = tmp + outModel(ll,innerIndex).y * Param{pp}(innerIndex,kk);
232 % end
233 % h{pp}(:,(kk-1)*size(outModel,1) + ll) = tmp;
234 % end
235 % end
236 %
237 % end
238 %
239 % for i = 1:length(params)
240 % d11(i).y = h{i}(:,1);
241 % d21(i).y = h{i}(:,2);
242 % d12(i).y = h{i}(:,3);
243 % d22(i).y = h{i}(:,4);
244 % end
245 %
246 % end
247
248 % scaling of PSD
249 % PSD = 2/(N*fs) * FFT *conj(FFT)
250 C11 = N*fs/2.*S11.y;
251 C22 = N*fs/2.*S22.y;
252 C12 = N*fs/2.*S12.y;
253 C21 = N*fs/2.*S21.y;
254
255 % compute elements of inverse cross-spectrum matrix
256 InvS11 = (C22./(C11.*C22 - C12.*C21));
257 InvS22 = (C11./(C11.*C22 - C12.*C21));
258 InvS12 = (C21./(C11.*C22 - C12.*C21));
259 InvS21 = (C12./(C11.*C22 - C12.*C21));
260
261
262 % compute Fisher Matrix
263 for i =1:length(params)
264 for j =1:length(params)
265
266 v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
267 v2v2 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
268 v1v2 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
269 v2v1 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
270
271 FisMat(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
272 end
273 end
274 utils.helper.msg(msg.IMPORTANT, sprintf('rank(FisMat) = %d', rank(FisMat)));
275
276 for kk = 1:numel(freqs)
277 % create input signal with power at single freq.
278 % depending on input channel
279 p = zeros(1,numel(freqs));
280 if channel == 1
281 p(kk) = sum(i1.y);
282 % create aos
283 i1single = ao(plist('Xvals',freqs,'Yvals',p));
284 i2single = ao(plist('Xvals',freqs,'Yvals',zeros(1,numel(freqs))));
285 elseif channel == 2
286 p(kk) = sum(i2.y);
287 % create aos
288 i1single = ao(plist('Xvals',freqs,'Yvals',zeros(1,numel(freqs))));
289 i2single = ao(plist('Xvals',freqs,'Yvals',p));
290 else
291 error('### wrong channel')
292 end
293
294 % compute Fisher Matrix for single frequencies
295 for i =1:length(params)
296 for j =1:length(params)
297
298 v1v1 = conj(d11(i).y.*i1single.y + d12(i).y.*i2single.y).*(d11(j).y.*i1single.y + d12(j).y.*i2single.y);
299 v2v2 = conj(d21(i).y.*i1single.y + d22(i).y.*i2single.y).*(d21(j).y.*i1single.y + d22(j).y.*i2single.y);
300 v1v2 = conj(d11(i).y.*i1single.y + d12(i).y.*i2single.y).*(d21(j).y.*i1single.y + d22(j).y.*i2single.y);
301 v2v1 = conj(d21(i).y.*i1single.y + d22(i).y.*i2single.y).*(d11(j).y.*i1single.y + d12(j).y.*i2single.y);
302
303 FisMatsingle(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
304 end
305 end
306 d(kk) = trace(FisMat\FisMatsingle)/numel(i1.x); % had to divide for num. freqs
307 % d(kk) = trace(pinv(FisMat)*FisMatsingle)/numel(i1.x); % had to divide for num. freqs
308
309 end
310
311 % % store Fisher Matrix for this run
312 % FM{k} = FisMat;
313 % % adding up
314 % FMall = FMall + FisMat;
315
316 % elseif ((numel(n(1).objs) == 3) && (numel(in.objs) == 4) && ~isempty(outModel))
317 % % this is only valid for the magnetic model, where we have 4 inputs
318 % % (corresponding to the 4 conformator waveforms) and 3 outputs
319 % % (corresponding to IFO.x12, IFO.eta1 and IFO.phi1). And there is a
320 % % contribution of an outModel converting the conformator waveforms
321 % % into forces and torques.
322 %
323 %
324 % % For other cases not implemented yet.
325 %
326 % % use signal fft to get frequency vector. Take into account signal
327 % % could be empty or set to zero
328 % % 1st channel
329 % freqs = fin.getObjectAtIndex(1,1).x;
330 %
331 % for ii = 1:numel(n.objs)
332 % for jj = ii:numel(n.objs)
333 % % Compute psd
334 % if (ii==jj)
335 % spec(ii,jj) = psd(n(k).getObjectAtIndex(ii), pl);
336 % S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
337 % else
338 % spec(ii,jj) = cpsd(n(k).getObjectAtIndex(ii),n(k).getObjectAtIndex(jj),pl);
339 % S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
340 % S2(jj,ii) = conj(S2(ii,jj));
341 % end
342 % end
343 % end
344 %
345 % S = matrix(S2,plist('shape',[numel(n.objs) numel(n.objs)]));
346 %
347 % % get some parameters used below
348 % fs = S.getObjectAtIndex(1,1).fs;
349 %
350 %
351 % if(~isempty(outModel))
352 % for lll=1:size(outModel,1)
353 % for kkk=1:size(outModel,2)
354 % outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[f1 f2]));
355 % end
356 % end
357 % end
358 %
359 % % Avoid numerical differentiation (faster for the magnetic case)
360 % Param{1} = [ 1 0 0 0;
361 % 0 0 0 0;
362 % 0 0 0 0;];
363 % Param{2} = [ 0 1 0 0;
364 % 0 0 0 0;
365 % 0 0 0 0;];
366 % Param{3} = [ 0 0 0 0;
367 % 0 0 1 0;
368 % 0 0 0 0;];
369 % Param{4} = [ 0 0 0 0;
370 % 0 0 0 0;
371 % 0 0 0 1;];
372 %
373 % % scaling of PSD
374 % % PSD = 2/(N*fs) * FFT *conj(FFT)
375 % for j = 1: numel(S.objs)
376 % % spectra to variance
377 % C(:,j) = (N*fs/2)*S.objs(j).data.getY;
378 % end
379 %
380 % detm = (C(:,1).*C(:,5).*C(:,9) + ...
381 % C(:,2).*C(:,6).*C(:,7) + ...
382 % C(:,3).*C(:,4).*C(:,8) -...
383 % C(:,7).*C(:,5).*C(:,3) -...
384 % C(:,8).*C(:,6).*C(:,1) -...
385 % C(:,9).*C(:,4).*C(:,2));
386 %
387 % InvS11 = (C(:,5).*C(:,9) - C(:,8).*C(:,6))./detm;
388 % InvS12 = -(C(:,4).*C(:,9) - C(:,7).*C(:,6))./detm;
389 % InvS13 = (C(:,4).*C(:,8) - C(:,7).*C(:,5))./detm;
390 % InvS21 = -(C(:,2).*C(:,9) - C(:,8).*C(:,3))./detm;
391 % InvS22 = (C(:,1).*C(:,9) - C(:,7).*C(:,3))./detm;
392 % InvS23 = -(C(:,1).*C(:,8) - C(:,7).*C(:,2))./detm;
393 % InvS31 = (C(:,2).*C(:,6) - C(:,5).*C(:,3))./detm;
394 % InvS32 = -(C(:,1).*C(:,6) - C(:,4).*C(:,3))./detm;
395 % InvS33 = (C(:,1).*C(:,5) - C(:,4).*C(:,2))./detm;
396 %
397 % for pp = 1:length(params)
398 % for ll = 1:size(outModel,1)
399 % for kk = 1:size(Param{pp},2)
400 % % index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
401 % tmp = 0;
402 % for innerIndex = 1:size(outModel,2)
403 % tmp = tmp + outModel(ll,innerIndex).y * Param{pp}(innerIndex,kk);
404 % end
405 % h{pp}(:,(kk-1)*size(outModel,1) + ll) = tmp;
406 % end
407 % end
408 %
409 % end
410 %
411 % for kk = 1:numel(in.objs)
412 % inV(:,kk) = fin.objs(kk).data.getY;
413 % end
414 %
415 %
416 % % compute Fisher Matrix
417 % for i =1:length(params)
418 % for j =1:length(params)
419 %
420 % for ll = 1:size(outModel,1)
421 % tmp = 0;
422 % for kk = 1:size(Param{1},2)
423 % tmp = tmp + h{i}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
424 % end
425 % v{i}(:,ll) = tmp;
426 % end
427 %
428 %
429 % for ll = 1:size(outModel,1)
430 % tmp = 0;
431 % for kk = 1:size(Param{1},2)
432 % tmp = tmp + h{j}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
433 % end
434 % v{j}(:,ll) = tmp;
435 % end
436 %
437 % v1v1 = conj(v{i}(:,1)).*v{j}(:,1);
438 % v1v2 = conj(v{i}(:,1)).*v{j}(:,2);
439 % v1v3 = conj(v{i}(:,1)).*v{j}(:,3);
440 % v2v1 = conj(v{i}(:,2)).*v{j}(:,1);
441 % v2v2 = conj(v{i}(:,2)).*v{j}(:,2);
442 % v2v3 = conj(v{i}(:,2)).*v{j}(:,3);
443 % v3v1 = conj(v{i}(:,3)).*v{j}(:,1);
444 % v3v2 = conj(v{i}(:,3)).*v{j}(:,2);
445 % v3v3 = conj(v{i}(:,3)).*v{j}(:,3);
446 %
447 % FisMat(i,j) = sum(real(InvS11.*v1v1 +...
448 % InvS12.*v1v2 +...
449 % InvS13.*v1v3 +...
450 % InvS21.*v2v1 +...
451 % InvS22.*v2v2 +...
452 % InvS23.*v2v3 +...
453 % InvS31.*v3v1 +...
454 % InvS32.*v3v2 +...
455 % InvS33.*v3v3));
456 % end
457 % end
458 % % store Fisher Matrix for this run
459 % FM{k} = FisMat;
460 % % adding up
461 % FMall = FMall + FisMat;
462 % else
463 % error('Implemented cases: 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
464 % end
465
466
467 % end
468
469 % inverse is the optimal covariance matrix
470 % if pseudoinv && isempty(tol)
471 % cov = pinv(FMall);
472 % elseif pseudoinv
473 % cov = pinv(FMall,tol);
474 % else
475 % cov = FMall\eye(size(FMall));
476 % end
477
478
479 % create AO
480 out = ao(plist('Xvals',freqs,'Yvals',d,'type','fsdata','fs',fs,'name',''));
481 % spectrum in the procinfo
482 if channel == 1
483 out.setProcinfo(plist('S11',S11));
484 elseif channel == 2
485 out.setProcinfo(plist('S22',S22));
486 end
487 % Fisher Matrix in the procinfo
488 out.setProcinfo(plist('FisMat',FisMat));
489
490 varargout{1} = out;
491 end
492
493 end
494
495
496 %--------------------------------------------------------------------------
497 % Get Info Object
498 %--------------------------------------------------------------------------
499 function ii = getInfo(varargin)
500 if nargin == 1 && strcmpi(varargin{1}, 'None')
501 sets = {};
502 pls = [];
503 else
504 sets = {'Default'};
505 pls = getDefaultPlist;
506 end
507 % Build info object
508 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: dispersion.m,v 1.1 2011/06/22 09:54:19 miquel Exp $', sets, pls);
509 end
510
511 %--------------------------------------------------------------------------
512 % Get Default Plist
513 %--------------------------------------------------------------------------
514 function plout = getDefaultPlist()
515 persistent pl;
516 if exist('pl', 'var')==0 || isempty(pl)
517 pl = buildplist();
518 end
519 plout = pl;
520 end
521
522 function pl = buildplist()
523 pl = plist.LPSD_PLIST;
524 pset(pl,'Navs',1)
525
526 p = plist({'f1', 'Initial frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
527 pl.append(p);
528
529 p = plist({'f2', 'Final frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
530 pl.append(p);
531
532 p = plist({'FitParamas', 'Parameters of the model'}, paramValue.EMPTY_STRING);
533 pl.append(p);
534
535 p = plist({'model','An array of matrix models'}, paramValue.EMPTY_STRING);
536 pl.append(p);
537
538 p = plist({'noise','An array of matrices with the cross-spectrum matrices'}, paramValue.EMPTY_STRING);
539 pl.append(p);
540
541 p = plist({'built-in','Symbolic models of the system as a string of built-in models'}, paramValue.EMPTY_STRING);
542 pl.append(p);
543
544 p = plist({'frequencies','Array of start/sop frequencies where the analysis is performed'}, paramValue.EMPTY_STRING);
545 pl.append(p);
546
547 p = plist({'pinv','Use the Penrose-Moore pseudoinverse'}, paramValue.TRUE_FALSE);
548 pl.append(p);
549
550 p = plist({'tol','Tolerance for the Penrose-Moore pseudoinverse'}, paramValue.EMPTY_DOUBLE);
551 pl.append(p);
552
553 p = plist({'diffStep','Numerical differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE);
554 pl.append(p);
555
556 end