Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/crb.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 % CRB computes the inverse of the Fisher Matrix
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: CRB computes the inverse of the Fisher Matrix
5 %
6 % CALL: bs = crb(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 - covariance matrix AO
14 %
15 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'crb')">Parameters Description</a>
16 %
17 % VERSION: $Id: crb.m,v 1.17 2011/10/07 08:19:55 miquel Exp $
18 %
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21 function varargout = crb(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 mdl = find(pl,'model');
52 mtxns = find(pl,'noise');
53 outModel = find(pl,'outModel');
54 bmdl = find(pl,'built-in');
55 f1 = find(pl,'f1');
56 f2 = find(pl,'f2');
57 pseudoinv = find(pl,'pinv');
58 tol = find(pl,'tol');
59 outNames = find(pl,'outNames');
60 inNames = find(pl,'inNames');
61
62 % Decide on a deep copy or a modify
63 in = copy(mtxs, nargout);
64 n = copy(mtxns, nargout);
65
66 % Get number of experiments
67 nexp = numel(in);
68
69 % fft
70 fin = fft(in);
71
72 % N should get before spliting, in order to convert correctly from psd to
73 % fft
74 N = length(fin(1).getObjectAtIndex(1).x);
75
76 % Get rid of fft f =0, reduce frequency range if needed
77 if ~isempty(f1) && ~isempty(f2)
78 fin = split(fin,plist('frequencies',[f1 f2]));
79 end
80
81 FMall = zeros(numel(params),numel(params));
82 % loop over experiments
83 for k = 1:nexp
84
85 utils.helper.msg(msg.IMPORTANT, sprintf('Analysis of experiment #%d',k), mfilename('class'), mfilename);
86
87 if (((numel(n(1).objs)) == 1) && (numel(in(1).objs) == 1))
88
89 % use signal fft to get frequency vector.
90 i1 = fin(k).getObjectAtIndex(1,1);
91 freqs = i1.x;
92
93 FisMat = utils.math.fisher_1x1(i1,n(k),mdl,params,numparams,freqs,N,pl,inNames,outNames);
94 % store Fisher Matrix for this run
95 FM{k} = FisMat;
96 % adding up
97 FMall = FMall + FisMat;
98
99 elseif (((numel(n(1).objs)) == 2) && (numel(in(1).objs) == 2))
100 % use signal fft to get frequency vector. Take into account signal
101 % could be empty or set to zero
102 % 1st channel
103 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
104 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
105 else
106 i1 = fin(k).getObjectAtIndex(1,1);
107 freqs = i1.x;
108 end
109 % 2nd channel
110 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
111 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
112 else
113 i2 = fin(k).getObjectAtIndex(2,1);
114 freqs = i2.x;
115 end
116
117 FisMat = utils.math.fisher_2x2(i1,i2,n(k),mdl,params,numparams,freqs,N,pl,inNames,outNames);
118 % store Fisher Matrix for this run
119 FM{k} = FisMat;
120 % adding up
121 FMall = FMall + FisMat;
122
123 elseif ((numel(n(1).objs) == 3) && (numel(in.objs) == 4) && ~isempty(outModel))
124 % this is only valid for the magnetic model, where we have 4 inputs
125 % (corresponding to the 4 conformator waveforms) and 3 outputs
126 % (corresponding to IFO.x12, IFO.eta1 and IFO.phi1). And there is a
127 % contribution of an outModel converting the conformator waveforms
128 % into forces and torques.
129
130
131 % For other cases not implemented yet.
132
133 % use signal fft to get frequency vector. Take into account signal
134 % could be empty or set to zero
135 % 1st channel
136 freqs = fin.getObjectAtIndex(1,1).x;
137
138 for ii = 1:numel(n.objs)
139 for jj = ii:numel(n.objs)
140 % Compute psd
141 if (ii==jj)
142 spec(ii,jj) = psd(n(k).getObjectAtIndex(ii), pl);
143 S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
144 else
145 spec(ii,jj) = cpsd(n(k).getObjectAtIndex(ii),n(k).getObjectAtIndex(jj),pl);
146 S2(ii,jj) = interp(spec(ii,jj),plist('vertices',freqs,'method','linear'));
147 S2(jj,ii) = conj(S2(ii,jj));
148 end
149 end
150 end
151
152 S = matrix(S2,plist('shape',[numel(n.objs) numel(n.objs)]));
153
154 % get some parameters used below
155 fs = S.getObjectAtIndex(1,1).fs;
156
157
158 if(~isempty(outModel))
159 for lll=1:size(outModel,1)
160 for kkk=1:size(outModel,2)
161 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[f1 f2]));
162 end
163 end
164 end
165
166 % Avoid numerical differentiation (faster for the magnetic case)
167 Param{1} = [ 1 0 0 0;
168 0 0 0 0;
169 0 0 0 0;];
170 Param{2} = [ 0 1 0 0;
171 0 0 0 0;
172 0 0 0 0;];
173 Param{3} = [ 0 0 0 0;
174 0 0 1 0;
175 0 0 0 0;];
176 Param{4} = [ 0 0 0 0;
177 0 0 0 0;
178 0 0 0 1;];
179
180 % scaling of PSD
181 % PSD = 2/(N*fs) * FFT *conj(FFT)
182 for j = 1: numel(S.objs)
183 % spectra to variance
184 C(:,j) = (N*fs/2)*S.objs(j).data.getY;
185 end
186
187 detm = (C(:,1).*C(:,5).*C(:,9) + ...
188 C(:,2).*C(:,6).*C(:,7) + ...
189 C(:,3).*C(:,4).*C(:,8) -...
190 C(:,7).*C(:,5).*C(:,3) -...
191 C(:,8).*C(:,6).*C(:,1) -...
192 C(:,9).*C(:,4).*C(:,2));
193
194 InvS11 = (C(:,5).*C(:,9) - C(:,8).*C(:,6))./detm;
195 InvS12 = -(C(:,4).*C(:,9) - C(:,7).*C(:,6))./detm;
196 InvS13 = (C(:,4).*C(:,8) - C(:,7).*C(:,5))./detm;
197 InvS21 = -(C(:,2).*C(:,9) - C(:,8).*C(:,3))./detm;
198 InvS22 = (C(:,1).*C(:,9) - C(:,7).*C(:,3))./detm;
199 InvS23 = -(C(:,1).*C(:,8) - C(:,7).*C(:,2))./detm;
200 InvS31 = (C(:,2).*C(:,6) - C(:,5).*C(:,3))./detm;
201 InvS32 = -(C(:,1).*C(:,6) - C(:,4).*C(:,3))./detm;
202 InvS33 = (C(:,1).*C(:,5) - C(:,4).*C(:,2))./detm;
203
204 for pp = 1:length(params)
205 for ll = 1:size(outModel,1)
206 for kk = 1:size(Param{pp},2)
207 % index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
208 tmp = 0;
209 for innerIndex = 1:size(outModel,2)
210 tmp = tmp + outModel(ll,innerIndex).y * Param{pp}(innerIndex,kk);
211 end
212 h{pp}(:,(kk-1)*size(outModel,1) + ll) = tmp;
213 end
214 end
215
216 end
217
218 for kk = 1:numel(in.objs)
219 inV(:,kk) = fin.objs(kk).data.getY;
220 end
221
222
223 % compute Fisher Matrix
224 for i =1:length(params)
225 for j =1:length(params)
226
227 for ll = 1:size(outModel,1)
228 tmp = 0;
229 for kk = 1:size(Param{1},2)
230 tmp = tmp + h{i}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
231 end
232 v{i}(:,ll) = tmp;
233 end
234
235
236 for ll = 1:size(outModel,1)
237 tmp = 0;
238 for kk = 1:size(Param{1},2)
239 tmp = tmp + h{j}(:,(kk-1)*size(outModel,1) + ll).*inV(:,kk);
240 end
241 v{j}(:,ll) = tmp;
242 end
243
244 v1v1 = conj(v{i}(:,1)).*v{j}(:,1);
245 v1v2 = conj(v{i}(:,1)).*v{j}(:,2);
246 v1v3 = conj(v{i}(:,1)).*v{j}(:,3);
247 v2v1 = conj(v{i}(:,2)).*v{j}(:,1);
248 v2v2 = conj(v{i}(:,2)).*v{j}(:,2);
249 v2v3 = conj(v{i}(:,2)).*v{j}(:,3);
250 v3v1 = conj(v{i}(:,3)).*v{j}(:,1);
251 v3v2 = conj(v{i}(:,3)).*v{j}(:,2);
252 v3v3 = conj(v{i}(:,3)).*v{j}(:,3);
253
254 FisMat(i,j) = sum(real(InvS11.*v1v1 +...
255 InvS12.*v1v2 +...
256 InvS13.*v1v3 +...
257 InvS21.*v2v1 +...
258 InvS22.*v2v2 +...
259 InvS23.*v2v3 +...
260 InvS31.*v3v1 +...
261 InvS32.*v3v2 +...
262 InvS33.*v3v3));
263 end
264 end
265 % store Fisher Matrix for this run
266 FM{k} = FisMat;
267 % adding up
268 FMall = FMall + FisMat;
269 else
270 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)');
271 end
272
273
274 end
275
276 % inverse is the optimal covariance matrix
277 if pseudoinv && isempty(tol)
278 cov = pinv(FMall);
279 elseif pseudoinv
280 cov = pinv(FMall,tol);
281 else
282 cov = FMall\eye(size(FMall));
283 end
284
285
286 % create AO
287 out = ao(cov);
288 % Fisher Matrix in the procinfo
289 out.setProcinfo(plist('FisMat',FM));
290
291 varargout{1} = out;
292 end
293
294
295 %--------------------------------------------------------------------------
296 % Get Info Object
297 %--------------------------------------------------------------------------
298 function ii = getInfo(varargin)
299 if nargin == 1 && strcmpi(varargin{1}, 'None')
300 sets = {};
301 pls = [];
302 else
303 sets = {'Default'};
304 pls = getDefaultPlist;
305 end
306 % Build info object
307 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: crb.m,v 1.17 2011/10/07 08:19:55 miquel Exp $', sets, pls);
308 end
309
310 %--------------------------------------------------------------------------
311 % Get Default Plist
312 %--------------------------------------------------------------------------
313 function plout = getDefaultPlist()
314 persistent pl;
315 if exist('pl', 'var')==0 || isempty(pl)
316 pl = buildplist();
317 end
318 plout = pl;
319 end
320
321 function pl = buildplist()
322 pl = plist.WELCH_PLIST;
323 pset(pl,'Navs',1)
324
325 p = plist({'f1', 'Initial frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
326 pl.append(p);
327
328 p = plist({'f2', 'Final frequency for the analysis'}, paramValue.EMPTY_DOUBLE);
329 pl.append(p);
330
331 p = plist({'FitParamas', 'Parameters of the model'}, paramValue.EMPTY_STRING);
332 pl.append(p);
333
334 p = plist({'model','An array of matrix models'}, paramValue.EMPTY_STRING);
335 pl.append(p);
336
337 p = plist({'noise','An array of matrices with the cross-spectrum matrices'}, paramValue.EMPTY_STRING);
338 pl.append(p);
339
340 p = plist({'built-in','Symbolic models of the system as a string of built-in models'}, paramValue.EMPTY_STRING);
341 pl.append(p);
342
343 p = plist({'frequencies','Array of start/sop frequencies where the analysis is performed'}, paramValue.EMPTY_STRING);
344 pl.append(p);
345
346 p = plist({'pinv','Use the Penrose-Moore pseudoinverse'}, paramValue.TRUE_FALSE);
347 pl.append(p);
348
349 p = plist({'tol','Tolerance for the Penrose-Moore pseudoinverse'}, paramValue.EMPTY_DOUBLE);
350 pl.append(p);
351
352 p = plist({'step','Numerical differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE);
353 pl.append(p);
354
355 p = plist({'ngrid','Number of points in the grid to compute the optimal differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE);
356 pl.append(p);
357
358 p = plist({'stepRanges','An array with upper and lower values for the parameters ranges. To be used to compute the optimal differentiation step for ssm models.'}, paramValue.EMPTY_DOUBLE);
359 pl.append(p);
360 end