comparison m-toolbox/classes/@matrix/mchNoisegenFilter.m @ 0:f0afece42f48

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 % MCHNOISEGENFILTER Construct a matrix filter from cross-spectral density matrix
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % FUNCTION: mchNoisegenFilter
5 %
6 % DESCRIPTION: Construct matrix filter from cross-spectral density. Such a
7 % filter can be used for multichannel noise generation in combination with
8 % the mchNoisegen method of the matrix class.
9 %
10 %
11 % CALL: fil = mchNoisegenFilter(mod, pl)
12 %
13 % INPUT:
14 % mod: is a matrix object containing the model for the target
15 % cross-spectral density matrix. Elements of mod must be fsdata
16 % analysis objects.
17 %
18 % OUTPUT:
19 % fil: is a matrix object containing the noise generating filter.
20 % Such a filter can be used to generate colored noise from
21 % uncorrelated unitary variance white time series. Fil can be a
22 % matrix of filterbanks objects or of parfrac objects according to
23 % the chosen output options.
24 %
25 % NOTE:
26 %
27 % The cross-spectral matrix is assumed to be frequency by frequency
28 % of the type:
29 %
30 % / csd11(f) csd12(f) \
31 % CSD(f) = | |
32 % \ csd21(f) csd22(f) /
33 %
34 %
35 %
36 % HISTORY: 22-04-2009 L Ferraioli
37 % Creation
38 %
39 % ------------------------------------------------------------------------
40 %
41 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mchNoisegenFilter')">Parameters Description</a>
42 %
43 % VERSION: $Id: mchNoisegenFilter.m,v 1.9 2011/05/16 10:36:18 luigi Exp $
44 %
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 function varargout = mchNoisegenFilter(varargin)
47
48 % Check if this is a call for parameters
49 if utils.helper.isinfocall(varargin{:})
50 varargout{1} = getInfo(varargin{3});
51 return
52 end
53
54 import utils.const.*
55 utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
56
57 % Collect input variable names
58 in_names = cell(size(varargin));
59 for ii = 1:nargin,in_names{ii} = inputname(ii);end
60
61 % Collect all ltpdauoh objects
62 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
63 [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist');
64
65 inhists = mtxs.hist;
66
67 % combine plists
68 pl = parse(pl, getDefaultPlist());
69 pl.getSetRandState();
70
71 % get elements out of the input matrix
72 csdm = copy(mtxs,nargout);
73 csdao = csdm.objs;
74
75 % Get parameters and set params for fit
76 fs = find(pl, 'fs');
77
78 target = lower(find(pl, 'targetobj')); % decide to perform s domain or z domain identification
79 % if target is parfrac output a matrix of parfarc objs (s domain
80 % identification)
81 % if target is miir output a matrix of filterbank parallel miir objects
82 % (z domain identification)
83
84 usesym = lower(find(pl, 'UseSym'));
85
86 if (fs == 0) && strcmpi(target,'miir')
87 error('### Please provide a valid sampling frequency for CSD constructor.');
88 elseif isempty(fs) && strcmpi(target,'miir')
89 error('### Please provide a valid sampling frequency for CSD constructor.');
90 end
91
92 % get units for filters
93 tgiunit = find(pl,'iunits');
94 tgounit = find(pl,'ounits');
95
96 % require filter initialization
97 initfilter = utils.prog.yes2true(find(pl, 'InitFilter'));
98
99 params = struct();
100
101 params.Nmaxiter = find(pl, 'MaxIter');
102 params.minorder = find(pl, 'MinOrder');
103 params.maxorder = find(pl, 'MaxOrder');
104 params.spolesopt = find(pl, 'PoleType');
105 params.weightparam = find(pl, 'Weights');
106
107 % set the target output
108 if strcmpi(target,'miir')
109 params.TargetDomain = 'z';
110 elseif strcmpi(target,'parfrac')
111 params.TargetDomain = 's';
112 else
113 error('### Unknown option for ''targetobj''.');
114 end
115
116 % Tolerance for MSE Value
117 lrscond = find(pl, 'FITTOL');
118 % give an error for strange values of lrscond
119 if lrscond<0
120 error('### Negative values for FITTOL are not allowed. ')
121 end
122 % handling data
123 lrscond = -1*log10(lrscond);
124 % give a warning for strange values of lrscond
125 if lrscond<0
126 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
127 end
128 params.lrscond = lrscond;
129
130 % Tolerance for the MSE relative variation
131 msevar = find(pl, 'MSEVARTOL');
132 % handling data
133 msevar = -1*log10(msevar);
134 % give a warning for strange values of msevar
135 if msevar<0
136 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
137 end
138 params.msevar = msevar;
139
140 if isempty(params.msevar)
141 params.ctp = 'chival';
142 else
143 params.ctp = 'chivar';
144 end
145
146 if(find(pl, 'plot'))
147 params.plot = 1;
148 else
149 params.plot = 0;
150 end
151
152 params.fs = fs;
153 params.dterm = 0; % it is better to fit without direct term
154
155 % check if symbolic calculation is required
156 if strcmpi(usesym,'on')
157 params.usesym = 1;
158 elseif strcmpi(usesym,'off')
159 params.usesym = 0;
160 else
161 error('### Unknown option for ''UseSym''.');
162 end
163
164 % extracting csd
165 if numel(csdao)==1 % one dimensional psd
166 csd = csdao.y;
167 freq = csdao.x;
168 dim = 'one';
169 else % multichannel system
170 dim = 'multi';
171 [nn,mm] = size(csdao);
172 if nn~=mm
173 error('### CSD Matrix must be square. ')
174 end
175 freq = csdao(1).x;
176 for ii = 1:nn
177 for jj = 1:nn
178 tcsd = csdao(ii,jj).y;
179 % willing to work with columns
180 [aa,bb] = size(tcsd);
181 if aa<bb
182 tcsd = tcsd.';
183 end
184 csd(ii,jj,:) = tcsd;
185 end
186 end
187
188 end
189
190
191 % call csd2tf
192 % ostruct is a struct array whose fields contain the residues and poles
193 % of estimated TFs. Since the fit is porformed on the columns of the TF
194 % matrix, each element of the array contains the residues and poles
195 % corresponding to the functions on the given column of the TF matrix.
196
197 %ostruct = utils.math.csd2tf(csd,freq,params);
198
199 ostruct = utils.math.csd2tf2(csd,freq,params);
200
201
202 % the filter for each channel is implemented by the rows of the TF matrix
203
204 switch dim
205 case 'one'
206
207 switch target
208 case 'miir'
209 % --- filter ---
210 res = ostruct.res;
211 poles = ostruct.poles;
212 % check if filter init is required
213 if initfilter
214 Zi = utils.math.getinitstate(res,poles,1,'mtd','svd');
215 else
216 Zi = zeros(size(res));
217 end
218
219 % construct a struct array of miir filters vectors
220 pfilts(numel(res),1) = miir;
221 for kk=1:numel(res)
222 ft = miir(res(kk), [ 1 -poles(kk)], fs);
223 ft.setIunits(unit(tgiunit));
224 ft.setOunits(unit(tgounit));
225 ft.setHistout(Zi(kk));
226 pfilts(kk,1) = ft;
227 clear ft
228 end
229 filt = filterbank(pfilts,'parallel');
230
231 csdm = matrix(filt);
232
233 % Add history
234 csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
235
236 case 'parfrac'
237 res = ostruct.res;
238 poles = ostruct.poles;
239
240 fbk = parfrac(res,poles,0);
241 fbk.setIunits(unit(tgiunit));
242 fbk.setOunits(unit(tgounit));
243
244 csdm = matrix(fbk);
245
246 % Add history
247 csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
248 end
249
250
251 case 'multi'
252
253 switch target
254 case 'miir'
255 % init filters array
256 %fbk(nn*nn,1) = filterbank;
257 %fbk = filterbank.newarray([nn nn]);
258
259 for zz=1:nn*nn % run over system dimension
260 % --- get column filter coefficients ---
261 % each column of mres\mpoles are the coefficients of a given filter
262 clear res poles
263 res = ostruct(zz).res;
264 poles = ostruct(zz).poles;
265
266 % construct a struct array of miir filters vectors
267 %ft(numel(res),1) = miir;
268 for kk=1:numel(res)
269 ft(kk,1) = miir(res(kk), [1 -poles(kk)], fs);
270 ft(kk,1).setIunits(unit(tgiunit));
271 ft(kk,1).setOunits(unit(tgounit));
272 end
273
274 fbk(zz,1) = filterbank(ft,'parallel');
275 clear ft
276
277 end
278
279 mfbk = reshape(fbk,nn,nn);
280
281 % check if filter init is required
282 if initfilter
283 ckidx = 0;
284 while ckidx<nn*nn
285 resv = [];
286 plsv = [];
287 for ii=1+ckidx:nn+ckidx
288 resv = [resv; ostruct(ii).res];
289 plsv = [plsv; ostruct(ii).poles];
290 end
291 % get init states
292 Zi = utils.math.getinitstate(resv,plsv,1,'mtd','svd');
293 clear resv plsv
294 % unpdate into the filters
295 for ii=1+ckidx:nn+ckidx
296 for kk=1:numel(mfbk)
297 mfbk(ii).filters(kk).setHistout(Zi(kk));
298 end
299 end
300 % update ckidx
301 ckidx = ckidx + nn;
302 end
303 end
304
305 csdm = matrix(mfbk);
306
307
308 % Add history
309 csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
310
311 case 'parfrac'
312 % init filters array
313 %fbk(nn*nn,1) = parfrac;
314
315 for zz=1:nn*nn % run over system dimension
316 % --- get column filter coefficients ---
317 % each column of mres\mpoles are the coefficients of a given filter
318 clear res poles
319 res = ostruct(zz).res;
320 poles = ostruct(zz).poles;
321
322 fbk(zz,1) = parfrac(res,poles,0);
323 fbk(zz,1).setIunits(unit(tgiunit));
324 fbk(zz,1).setOunits(unit(tgounit));
325
326 end
327
328 mfbk = reshape(fbk,nn,nn);
329
330 csdm = matrix(mfbk);
331
332 % Add history
333 csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
334 end
335
336 end
337
338 % Set properties from the plist
339 warning('off', utils.const.warnings.METHOD_NOT_FOUND);
340 % remove parameters we already used
341 pl_set = copy(pl,1);
342 if pl_set.isparam('fs')
343 pl_set.remove('fs');
344 end
345 if pl_set.isparam('targetobj')
346 pl_set.remove('targetobj');
347 end
348 if pl_set.isparam('UseSym')
349 pl_set.remove('UseSym');
350 end
351 if pl_set.isparam('iunits')
352 pl_set.remove('iunits');
353 end
354 if pl_set.isparam('ounits')
355 pl_set.remove('ounits');
356 end
357 if pl_set.isparam('InitFilter')
358 pl_set.remove('InitFilter');
359 end
360 if pl_set.isparam('MaxIter')
361 pl_set.remove('MaxIter');
362 end
363 if pl_set.isparam('MinOrder')
364 pl_set.remove('MinOrder');
365 end
366 if pl_set.isparam('MaxOrder')
367 pl_set.remove('MaxOrder');
368 end
369 if pl_set.isparam('PoleType')
370 pl_set.remove('PoleType');
371 end
372 if pl_set.isparam('Weights')
373 pl_set.remove('Weights');
374 end
375 if pl_set.isparam('FITTOL')
376 pl_set.remove('FITTOL');
377 end
378 if pl_set.isparam('MSEVARTOL')
379 pl_set.remove('MSEVARTOL');
380 end
381 if pl_set.isparam('plot')
382 pl_set.remove('plot');
383 end
384 csdm.setProperties(pl_set);
385 warning('on', utils.const.warnings.METHOD_NOT_FOUND);
386
387 % output data
388 varargout{1} = csdm;
389
390 end
391
392 %--------------------------------------------------------------------------
393 % Get Info Object
394 %--------------------------------------------------------------------------
395 function ii = getInfo(varargin)
396 if nargin == 1 && strcmpi(varargin{1}, 'None')
397 sets = {};
398 pl = [];
399 else
400 sets = {'Default'};
401 pl = getDefaultPlist;
402 end
403 % Build info object
404 ii = minfo(mfilename, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: mchNoisegenFilter.m,v 1.9 2011/05/16 10:36:18 luigi Exp $', sets, pl);
405 ii.setArgsmin(1);
406 ii.setOutmin(1);
407 ii.setOutmax(1);
408 end
409
410 %--------------------------------------------------------------------------
411 % Get Default Plist
412 %--------------------------------------------------------------------------
413 function plout = getDefaultPlist()
414 persistent pl;
415 if exist('pl', 'var')==0 || isempty(pl)
416 pl = buildplist();
417 end
418 plout = pl;
419 end
420
421 function pl = buildplist()
422
423 pl = plist();
424
425 % Target Objects
426 p = param({'targetobj', ['Choose the type of output objects:<ul>',...
427 '<li>''miir'' output a matrix containing filterbanks of parallel miir filters</li>',...
428 '<li>''parfrac'' output a matrix containing parafracs objects</li>']}, ...
429 {1, {'miir','parfrac'}, paramValue.OPTIONAL});
430 pl.append(p);
431
432 % Fs
433 p = param({'fs', 'The sampling frequency of the discrete filters.'}, {1, {1}, paramValue.OPTIONAL});
434 pl.append(p);
435
436 % Iunits
437 p = param({'iunits', 'The unit to set as input unit for the output filters'}, paramValue.EMPTY_STRING);
438 pl.append(p);
439
440 % Ounits
441 p = param({'ounits', 'The unit to set as output unit for the output filters'}, paramValue.EMPTY_STRING);
442 pl.append(p);
443
444 % Plot
445 p = param({'InitFilter', 'Initialize filters (works only for miir objects) to cope with startup transients.'}, paramValue.TRUE_FALSE);
446 p.val.setValIndex(1);
447 pl.append(p);
448
449 % Max Iter
450 p = param({'MaxIter', 'Maximum number of fit iterations.'}, {1, {50}, paramValue.OPTIONAL});
451 pl.append(p);
452
453 % Pole type
454 p = param({'PoleType',['Choose the pole type for fitting initialization:<ul>',...
455 '<li>1 == use real starting poles</li>',...
456 '<li>2 == generates complex conjugate poles of the type <tt>a.*exp(theta*pi*j)</tt> with <tt>theta = linspace(0,pi,N/2+1)</tt></li>',...
457 '<li>3 == generates complex conjugate poles of the type <tt>a.*exp(theta*pi*j)</tt> with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ul>']}, ...
458 {1, {1, 2, 3}, paramValue.SINGLE});
459 pl.append(p);
460
461 % Min order
462 p = param({'MinOrder','Minimum order to fit with.'}, {1, {7}, paramValue.OPTIONAL});
463 pl.append(p);
464
465 % Max Order
466 p = param({'MaxOrder','Maximum order to fit with.'}, {1, {35}, paramValue.OPTIONAL});
467 pl.append(p);
468
469 % Weights
470 p = param({'Weights',['Choose weighting for the fit:<ul>',...
471 '<li> 1 == equal weights for each point</li>',...
472 '<li> 2 == weight with <tt>1/abs(model)</tt></li>',...
473 '<li> 3 == weight with <tt>1/abs(model).^2</tt></li>',...
474 '<li> 4 == weight with inverse of the square mean spread of the model</li></ul>']}, {3, {1 2 3 4}, paramValue.SINGLE});
475 pl.append(p);
476
477 % Plot
478 p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.TRUE_FALSE);
479 p.val.setValIndex(2);
480 pl.append(p);
481
482 % MSE Vartol
483 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the realtive variation of the mean squared error is<br>',...
484 'smaller than the value specified. This option is useful for finding the minimum of Chi squared.']}, ...
485 {1, {1e-1}, paramValue.OPTIONAL});
486 pl.append(p);
487
488 % FIT TOL
489 p = param({'FITTOL',['Mean Squared Error Value - Check if the mean squared error value <br>',...
490 ' is lower than the value specified.']}, {1, {1e-2}, paramValue.OPTIONAL});
491 pl.append(p);
492
493 % UseSym
494 p = param({'UseSym', ['Use symbolic calculation in eigen-decomposition.<ul>'...
495 '<li>''on'' - uses symbolic math toolbox calculation<br>'...
496 'for poles stabilization</li>'...
497 '<li>''off'' - perform double-precision calculation<br>'...
498 'for poles stabilization</li>']}, {1, {'on','off'}, paramValue.SINGLE});
499 pl.append(p);
500
501
502
503
504 end
505
506