Mercurial > hg > ltpda
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 |