Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/crbound.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 % CRBOUND computes the inverse of the Fisher Matrix
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: CRBOUND computes the Cramer-Rao lowe bound for parametric
5 % model given the input signals and the noise. The method
6 % accepts 2D (2 input/2 output) models and 1D models
7 % (1 input/1 output) and (2 input/1 output)
8 %
9 % CALL: bs = crbound(in,noise,pl)
10 %
11 % INPUTS: in - analysis objects with input signals to the system
12 % (x1) if 1 input
13 % (x2) if 2 input
14 % noise - analysis objects with measured noise
15 % (x1) if 1 input: S11
16 % (x4) if 2 input: (S11,S12,S21,S22)
17 % model - symbolic model (smodel) containing the evaluated
18 % transfer function models
19 % (x1) if 1 input/1 output
20 % (x2) if 2 input/1 output
21 % (x4) if 2 input/2 output
22 % pl - parameter list
23 %
24 % OUTPUTS: bs - covariance matrix AO
25 %
26 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'crbound')">Parameters Description</a>
27 %
28 % VERSION: $Id: crbound.m,v 1.13 2011/04/08 08:56:17 hewitson Exp $
29 %
30 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31
32 function varargout = crbound(varargin)
33
34 % Check if this is a call for parameters
35 if utils.helper.isinfocall(varargin{:})
36 varargout{1} = getInfo(varargin{3});
37 return
38 end
39
40 import utils.const.*
41 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
42
43 % Method can not be used as a modifier
44 if nargout == 0
45 error('### crbound cannot be used as a modifier. Please give an output variable.');
46 end
47
48 % Collect input variable names
49 in_names = cell(size(varargin));
50 for ii = 1:nargin,in_names{ii} = inputname(ii);end
51
52 % Collect all AOs smodels and plists
53 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
54 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
55 % Combine plists
56 pl = parse(pl, getDefaultPlist);
57
58 % get params
59 params = find(pl,'params');
60 numparams = find(pl,'paramsValues');
61 mdl2 = find(pl,'models');
62 bmdl = find(pl,'built-in');
63 f1 = find(pl,'f1');
64 f2 = find(pl,'f2');
65 freqs = find(pl,'frequencies');
66 pseudoinv = find(pl,'pinv');
67 tol = find(pl,'tol');
68 % ninputs = find(pl,'Ninputs');
69
70
71 % Decide on a deep copy or a modify
72 bs = copy(as, nargout);
73 if ~isempty(mdl2)
74 mdl = copy(mdl2, nargout);
75 end
76
77 if numel(bs) == 2
78 % implementing here for 1D
79 in(1) = bs(1);
80 n(1) = bs(2);
81
82
83
84 % fft
85 i1 = fft(in(1));
86 Nsamples = length(i1.x);
87 if ~isempty(f1) && ~isempty(f2)
88 i1 = split(i1,plist('frequencies',[f1 f2]));
89 elseif ~isempty(freqs)
90 i1 = split(i1,plist('frequencies',freqs));
91 i1 = join(i1);
92 end
93
94 % get frequency vector
95 f = i1.x;
96
97 if ~isempty(mdl)
98 % get model
99 h(1) = mdl;
100 h(1).setXvals(f);
101 else
102 error('### One model should be introduced, please check the help.')
103 end
104
105 elseif numel(bs) == 3
106 if numel(mdl) == 2
107 in(1) = bs(1);
108 in(2) = bs(2);
109 n(1) = bs(3);
110
111 % fft
112 i1 = fft(in(1));
113 Nsamples = length(i1.x);
114 i2 = fft(in(2));
115
116 % Get rid of fft f =0, reduce frequency range if needed
117 if ~isempty(f1) && ~isempty(f2)
118 i1 = split(i1,plist('frequencies',[f1 f2]));
119 i2 = split(i2,plist('frequencies',[f1 f2]));
120 elseif ~isempty(freqs)
121 i1 = split(i1,plist('frequencies',freqs));
122 i1 = join(i1);
123 i2 = split(i2,plist('frequencies',freqs));
124 i2 = join(i2);
125 end
126
127 % get frequency vector
128 f = i1.x;
129
130 if ~isempty(mdl)
131 % compute built-in matrix
132 h(1) = mdl(1);
133 h(1).setXvals(f);
134 h(2) = mdl(2);
135 h(2).setXvals(f);
136 else
137 error('### One model should be introduced, please check the help.')
138 end
139 elseif numel(mdl) == 1
140 error('Check the model or use crbound with the matrix class')
141 end
142 % Get input objects
143 elseif numel(bs) == 4
144 in(1) = bs(1);
145 in(2) = bs(2);
146 n(1) = bs(3);
147 n(2) = bs(4);
148
149 % fft
150 i1 = fft(in(1));
151 Nsamples = length(i1.x);
152 i2 = fft(in(2));
153
154 % Get rid of fft f =0, reduce frequency range if needed
155 if ~isempty(f1) && ~isempty(f2)
156 i1 = split(i1,plist('frequencies',[f1 f2]));
157 i2 = split(i2,plist('frequencies',[f1 f2]));
158 elseif ~isempty(freqs)
159 i1 = split(i1,plist('frequencies',freqs));
160 i1 = join(i1);
161 i2 = split(i2,plist('frequencies',freqs));
162 i2 = join(i2);
163 end
164
165 % get frequency vector
166 f = i1.x;
167
168 if ~isempty(bmdl)
169 % compute built-in smodels
170 for i = 1:4
171 if strcmp(bmdl{i},'0');
172 h(i) = smodel('0');
173 h(i).setXvar('f');
174 h(i).setXvals(f);
175 else
176 h(i) = smodel(plist('built-in',bmdl{i},'f',f));
177 % set all params to all models. It is not true but harmless
178 for k = 1:numel(params)
179 vecparams(k) = {numparams(k)*ones(size(f))};
180 end
181 h(i).setParams(params,vecparams);
182 end
183 end
184 elseif ~isempty(mdl)
185 % compute built-in matrix
186 h(1) = mdl{1};
187 h(1).setXvals(f);
188 h(2) = mdl{2};
189 h(2).setXvals(f);
190 h(3) = mdl{3};
191 h(3).setXvals(f);
192 h(4) = mdl{4};
193 h(4).setXvals(f);
194 end
195
196 else
197 error('### The number of input in crbound objects is not correct, please check the help.')
198 end
199
200 % Set params
201 for i = 1:numel(h)
202 h(i).setParams(params,numparams);
203 end
204
205 if numel(bs) == 2
206 % Compute psd
207 n1 = psd(n(1), pl);
208
209 % interpolate to fft frequencies
210 S11 = interp(n1,plist('vertices',f));
211
212 for i = 1:length(params)
213 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
214 % differentiate symbolically
215 dH11 = diff(h(1),params{i});
216 % evaluate
217 d11(i) = eval(dH11);
218 end
219
220 % get some parameters used below
221 fs = S11.fs;
222 % N = len(S11);
223 Navs = find(pl,'Navs');
224
225 % scaling of PSD
226 % PSD = 2/(N*fs) * FFT *conj(FFT)
227 C11 = Nsamples*fs/2.*S11.y;
228
229 InvS11 = 1./C11;
230
231 % compute Fisher Matrix
232 for i =1:length(params)
233 for j =1:length(params)
234 v1v1 = conj(d11(i).y.*i1.y).*(d11(j).y.*i1.y);
235
236 FisMat(i,j) = sum(real(InvS11.*v1v1));
237 end
238 end
239
240 elseif numel(bs) == 3
241
242 if numel(mdl) == 2
243 % Compute psd
244 n1 = psd(n(1), pl);
245
246 % interpolate to fft frequencies
247 S11 = interp(n1,plist('vertices',f));
248
249 for i = 1:length(params)
250 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
251 % differentiate symbolically
252 dH11 = diff(h(1),params{i});
253 dH12 = diff(h(2),params{i});
254
255 % evaluate
256 d11(i) = eval(dH11);
257 d12(i) = eval(dH12);
258 end
259
260 % get some parameters used below
261 fs = S11.fs;
262 % N = len(S11);
263 Navs = find(pl,'Navs');
264
265 % scaling of PSD
266 % PSD = 2/(N*fs) * FFT *conj(FFT)
267 C11 = Nsamples*fs/2.*S11.y;
268
269 % compute elements of inverse cross-spectrum matrix
270 InvS11 = 1./C11;
271
272 % compute Fisher Matrix
273 for i =1:length(params)
274 for j =1:length(params)
275 v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
276
277 FisMat(i,j) = sum(real(InvS11.*v1v1));
278 end
279 end
280
281 elseif numel(mdl) == 1
282 error('Please check model sizes')
283 end
284
285
286 elseif numel(bs) == 4
287
288 % Compute psd
289 n1 = psd(n(1), pl);
290 n2 = psd(n(2), pl);
291 n12 = cpsd(n(1),n(2), pl);
292
293 % interpolate to fft frequencies
294 S11 = interp(n1,plist('vertices',f));
295 S12 = interp(n12,plist('vertices',f));
296 S22 = interp(n2,plist('vertices',f));
297 S21 = conj(S12);
298
299 for i = 1:length(params)
300 utils.helper.msg(msg.IMPORTANT, sprintf('computing symbolic differentiation with respect %s',params{i}), mfilename('class'), mfilename);
301 % differentiate symbolically
302 dH11 = diff(h(1),params{i});
303 dH12 = diff(h(2),params{i});
304 dH21 = diff(h(3),params{i});
305 dH22 = diff(h(4),params{i});
306 % evaluate
307 d11(i) = eval(dH11);
308 d12(i) = eval(dH12);
309 d21(i) = eval(dH21);
310 d22(i) = eval(dH22);
311 end
312
313 % get some parameters used below
314 fs = S11.fs;
315 % N = len(S11);
316 Navs = find(pl,'Navs');
317
318 % scaling of PSD
319 % PSD = 2/(N*fs) * FFT *conj(FFT)
320 C11 = Nsamples*fs/2.*S11.y;
321 C22 = Nsamples*fs/2.*S22.y;
322 C12 = Nsamples*fs/2.*S12.y;
323 C21 = Nsamples*fs/2.*S21.y;
324
325 % compute elements of inverse cross-spectrum matrix
326 InvS11 = (C22./(C11.*C22 - C12.*C21));
327 InvS22 = (C11./(C11.*C22 - C12.*C21));
328 InvS12 = (C21./(C11.*C22 - C12.*C21));
329 InvS21 = (C12./(C11.*C22 - C12.*C21));
330
331 % compute Fisher Matrix
332 for i =1:length(params)
333 for j =1:length(params)
334
335 v1v1 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
336 v2v2 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
337 v1v2 = conj(d11(i).y.*i1.y + d12(i).y.*i2.y).*(d21(j).y.*i1.y + d22(j).y.*i2.y);
338 v2v1 = conj(d21(i).y.*i1.y + d22(i).y.*i2.y).*(d11(j).y.*i1.y + d12(j).y.*i2.y);
339
340 FisMat(i,j) = sum(real(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1));
341 end
342 end
343
344 end
345
346 % inverse is the optimal covariance matrix
347 if pseudoinv && isempty(tol)
348 cov = pinv(FisMat);
349 elseif pseudoinv
350 cov = pinv(FisMat,tol);
351 else
352 cov = FisMat\eye(size(FisMat));
353 end
354
355 % create AO
356 out = ao(cov);
357 % Fisher Matrix in the procinfo
358 out.setProcinfo(plist('FisMat',FisMat));
359
360 varargout{1} = out;
361 end
362
363
364 %--------------------------------------------------------------------------
365 % Get Info Object
366 %--------------------------------------------------------------------------
367 function ii = getInfo(varargin)
368 if nargin == 1 && strcmpi(varargin{1}, 'None')
369 sets = {};
370 pls = [];
371 else
372 sets = {'Default'};
373 pls = getDefaultPlist;
374 end
375 % Build info object
376 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: crbound.m,v 1.13 2011/04/08 08:56:17 hewitson Exp $', sets, pls);
377 end
378
379 %--------------------------------------------------------------------------
380 % Get Default Plist
381 %--------------------------------------------------------------------------
382
383 function plout = getDefaultPlist()
384 persistent pl;
385 if exist('pl', 'var')==0 || isempty(pl)
386 pl = buildplist();
387 end
388 plout = pl;
389 end
390
391 function pl = buildplist()
392 pl = plist.WELCH_PLIST;
393 pset(pl,'Navs',1)
394
395 p = plist({'params', 'Parameters of the model'}, paramValue.EMPTY_STRING);
396 pl.append(p);
397
398 p = plist({'models','Symbolic models of the system'}, paramValue.EMPTY_STRING);
399 pl.append(p);
400
401 p = plist({'frequencies','Array of start/sop frequencies where the analysis is performed'}, paramValue.EMPTY_STRING);
402 pl.append(p);
403
404 p = plist({'pinv','Use the Penrose-Moore pseudoinverse'}, paramValue.TRUE_FALSE);
405 pl.append(p);
406
407 p = plist({'tol','Tolerance for the Penrose-Moore pseudoinverse'}, paramValue.EMPTY_DOUBLE);
408 pl.append(p);
409
410 end