Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/noisegen2D.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 % NOISEGEN2D generates cross correleted colored noise from white noise.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: noisegen2D can work in two different modes:
5 %
6 % ------------------------------------------------------------------------
7 %
8 % 1) Generates colored noise from white noise with a given cross spectrum.
9 % This mode correspond to the 'Default' set for the method (see the
10 % list of parameters).
11 %
12 % The coloring filter is constructed by a fitting procedure to the models
13 % provided. If no model is provided an error is prompted. The
14 % cross-spectral matrix is assumed to be frequency by frequency of the
15 % type:
16 %
17 % / csd11(f) csd12(f) \
18 % CSD(f) = | |
19 % \ csd21(f) csd22(f) /
20 %
21 % Note: The function output colored noise data with one-sided
22 % csd corresponding to the model provided.
23 %
24 % ALGORITHM:
25 % 1) Fit a set of partial fraction z-domain filters
26 % 2) Convert to array of MIIR filters
27 % 3) Filter time-series in parallel
28 % The filtering process is:
29 % b(1) = Filt11(a(1)) + Filt12(a(2))
30 % b(2) = Filt21(a(1)) + Filt22(a(2))
31 %
32 %
33 % CALL: b = noisegen2D(a, pl) % returns colored time-series AOs
34 % b = noisegen2D(a, pl)
35 % [b1,b2] = noisegen2D(a1, a2, pl)
36 % [b1,b2,...,bn] = noisegen2D(a1,a2,...,an, pl);
37 % Note: this method cannot be used as a modifier, the
38 % call a.noisegen2D(pl) is forbidden
39 %
40 % INPUT:
41 %
42 % - a is at least a couple of time series analysis objects
43 % - pl is a parameter list, see the list of accepted
44 % parameters below
45 %
46 % OUTPUT:
47 %
48 % - b are a couple of colored time-series AOs. The coloring
49 % filters used are stored in the objects procinfo field under
50 % the parameters:
51 % - b(1): 'Filt11' and 'Filt12'
52 % - b(2): 'Filt21' and 'Filt22'
53 % ------------------------------------------------------------------------
54 %
55 % 2) Generates coloring filter
56 % This mode correspond to the 'Filter' set for the method (see the
57 % list of parameters).
58 %
59 % The coloring filter is constructed by a fitting procedure to the models
60 % provided. The cross-spectral matrix is assumed to be frequency by
61 % frequency of the type:
62 %
63 % / csd11(f) csd12(f) \
64 % CSD(f) = | |
65 % \ csd21(f) csd22(f) /
66 %
67 % ALGORITHM:
68 % 1) Fit a set of partial fraction z-domain filters
69 % 2) Convert to array of MIIR filters
70 %
71 %
72 % CALL: fil = noisegen2D(csd11,csd12,csd21,csd22, pl)
73 % fil = noisegen2D(csd11,csd12,csd22, pl)
74 % Note: this method cannot be used as a modifier, the
75 % call a.noisegen2D(pl) is forbidden
76 %
77 % INPUT:
78 %
79 % - csd11, csd12, csd21,csd22 are the terms of the
80 % cross-spectral matrix. They must be frequency series
81 % analysis objects.
82 % - pl is a parameter list, see the list of accepted
83 % parameters below
84 %
85 % OUTPUT:
86 %
87 % - fil is a matrix object which represent a two dimensional
88 % filter. The elements of fil are filterbanks parallel
89 % objects of miir filters. Filters are initialized to
90 % avoid startup transients.
91 %
92 % ------------------------------------------------------------------------
93 %
94 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'noisegen2D')">Parameters Description</a>
95 %
96 % VERSION: $Id: noisegen2D.m,v 1.34 2011/04/08 08:56:13 hewitson Exp $
97 %
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99
100 function varargout = noisegen2D(varargin)
101
102 % Check if this is a call for parameters
103 if utils.helper.isinfocall(varargin{:})
104 varargout{1} = getInfo(varargin{3});
105 return
106 end
107
108 import utils.const.*
109 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
110
111 % Collect input variable names
112 in_names = cell(size(varargin));
113 for ii = 1:nargin,in_names{ii} = inputname(ii);end
114
115 % Collect all AOs and plists
116 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
117 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
118
119 % Decide on a deep copy or a modify
120 bs = copy(as, nargout);
121 inhists = [as.hist];
122
123 % combine plists
124 % define input/output combinations. Different combination are
125 % 1) input tsdata and csd# into the plist, output are colored tsdata
126 % 2) input fsdata, output is a coloring filter (in a matrix)
127 if isempty(pl) % no model input, get model from input
128 setpar = 'Filter';
129 elseif numel(as)==3 && isa(as(1).data,'fsdata') % get model from input
130 setpar = 'Filter';
131 elseif numel(as)==4 && isa(as(1).data,'fsdata') % get model from input
132 setpar = 'Filter';
133 else % get model from plist, output tsdata, back compatibility mode
134 setpar = 'Default';
135 end
136
137 pl = parse(pl, getDefaultPlist(setpar));
138 pl.getSetRandState();
139
140 if nargout == 0
141 error('### noisegen2D cannot be used as a modifier. Please give an output variable.');
142 end
143
144 % Check the number of input AOs
145 if numel(bs)==1
146 error('!!! One input AO! At least two independent white noise time series or three frequency series are needed');
147 end
148
149 switch lower(setpar)
150 case 'default' % back compatibility mode
151
152 % this copy will be used for data filtering
153 cs = copy(as, nargout);
154
155 % Extract necessary model parameters
156 csd11 = find(pl, 'csd11');
157 csd12 = find(pl, 'csd12');
158 csd21 = find(pl, 'csd21');
159 csd22 = find(pl, 'csd22');
160
161 if ~isempty(csd12) && isempty(csd21)
162 csd21 = conj(csd12);
163 elseif ~isempty(csd21) && isempty(csd12)
164 csd12 = conj(csd21);
165 elseif isempty(csd12) && isempty(csd21)
166 error('!!! One of the parameters ''csd12'' or ''csd21'' must be not empty!')
167 end
168
169 % Checks that the input AOs come in pairs
170 odc = 0;
171 if rem(numel(as),2)
172 warning('The input AOs must come in pairs! Skipping AO %s during calculation', ao_invars{end});
173 odc = 1;
174 end
175
176 % Loop over input AOs to check for non time series objects
177 fsv = zeros(numel(bs),1);
178 for jj=1:numel(bs)
179 if ~isa(bs(jj).data, 'tsdata')
180 error('!!! %s expects ao/tsdata objects. ', mfilename);
181 end
182 fsv(jj,1) = bs(jj).fs; % collecting sampling frequencies
183 end
184 % Check that input Aos have the same sampling frequency
185 if any(diff(fsv))
186 error('!!! Sampling frequency must be the same for all input objects')
187 end
188
189 case 'filter'
190
191 % get models from inputs
192 if numel(bs)==3
193 csd11 = bs(1);
194 csd12 = bs(2);
195 csd21 = conj(csd12);
196 csd22 = bs(3);
197 elseif numel(bs)==2
198 error('!!! A number of fsdata ao between 3 and 4 must be given as input')
199 else
200 csd11 = bs(1);
201 csd12 = bs(2);
202 csd21 = bs(3);
203 csd22 = bs(4);
204 end
205
206 fsv = find(pl,'fs');
207 Iunits = find(pl, 'Iunits');
208 Ounits = find(pl, 'Ounits');
209
210 otherwise
211
212 error('!!! Something with the input is going wrong! Check function help for details on how input data properly')
213
214 end
215
216
217 % ----------------------------------
218 % 1) - Fitting the models to identify the innovation filters
219
220 % Build input structure for psd2tf
221 params = struct();
222
223 params.idtp = 1;
224 params.Nmaxiter = find(pl, 'MaxIter');
225 params.minorder = find(pl, 'MinOrder');
226 params.maxorder = find(pl, 'MaxOrder');
227 params.spolesopt = find(pl, 'PoleType');
228 params.usesym = find(pl, 'UseSym');
229 params.spy = find(pl, 'Disp');
230
231 % check for weights
232 wobj = find(pl, 'Weights');
233 if isa(wobj,'ao')
234 warning('Using externally provided weights.')
235 params.weightparam = 0;
236 % check external weights dimensions
237 if numel(wobj)~= 4
238 erro('!!! Provide a weight for each CSD element')
239 end
240 for ii=1:4
241 twobj = wobj.index(ii).y;
242 % willing to work with columns
243 [aaw,bbw] = size(twobj);
244 if aaw<bbw
245 twobj = twobj.';
246 end
247 wobj2(:,ii) = twobj;
248 end
249 params.extweights = wobj2;
250 else
251 params.weightparam = wobj;
252 end
253
254 % Tolerance for MSE Value
255 lrscond = find(pl, 'FITTOL');
256 % give an error for strange values of lrscond
257 if lrscond<0
258 error('!!! Negative values for FITTOL are not allowed !!!')
259 end
260 % handling data
261 lrscond = -1*log10(lrscond);
262 % give a warning for strange values of lrscond
263 if lrscond<0
264 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
265 end
266 params.lrscond = lrscond;
267
268 % Tolerance for the MSE relative variation
269 msevar = find(pl, 'MSEVARTOL');
270 % handling data
271 msevar = -1*log10(msevar);
272 % give a warning for strange values of msevar
273 if msevar<0
274 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
275 end
276 params.msevar = msevar;
277
278 if isempty(params.msevar)
279 params.ctp = 'chival';
280 else
281 params.ctp = 'chivar';
282 end
283
284 if(find(pl, 'plot'))
285 params.plot = 1;
286 else
287 params.plot = 0;
288 end
289
290 params.fs = fsv(1,1);
291 params.dterm = 0; % it is better to fit without direct term
292
293 % call psd2tf
294 ostruct = utils.math.psd2tf(csd11.y,csd12.y,csd21.y,csd22.y,csd11.x,params);
295
296 % ----------------------------------
297 % 2) - Convert into MIIR filters
298
299 fs = fsv(1,1);
300
301 % get init states for coloring filters
302 mres13 = [ostruct(1).res ostruct(3).res];
303 mres24 = [ostruct(2).res ostruct(4).res];
304 mpoles13 = [ostruct(1).poles ostruct(3).poles];
305 mpoles24 = [ostruct(4).poles ostruct(4).poles];
306
307 % initialize filters to cope with starting transients
308 % Zi1 = zeros(numel(mres13(:,1)),1);
309 % Zi3 = Zi1;
310
311 Zi = utils.math.getinitstate(mres13,mpoles13,1,'mtd','svd');
312 Zi1 = Zi(1:numel(mres13(:,1)));
313 Zi3 = Zi(numel(mres13(:,1))+1:2*numel(mres13(:,1)));
314
315 % Zi2 = zeros(numel(mres24(:,1)),1);
316 % Zi4 = Zi2;
317
318 Zi = utils.math.getinitstate(mres24,mpoles24,1,'mtd','svd');
319 Zi2 = Zi(1:numel(mres24(:,1)));
320 Zi4 = Zi(numel(mres24(:,1))+1:2*numel(mres24(:,1)));
321
322 % --- filter 1 ---
323 res = mres13(:,1);
324 poles = mpoles13(:,1);
325 % construct a struct array of miir filters vectors
326 pfilts1 = [];
327 for kk=1:numel(res)
328 ft = miir(res(kk), [ 1 -poles(kk)], fs);
329 ft.setHistout(Zi1(kk));
330 pfilts1 = [pfilts1 ft];
331 end
332
333 % --- filter 2 ---
334 res = mres24(:,1);
335 poles = mpoles24(:,1);
336 % construct a struct array of miir filters vectors
337 pfilts2 = [];
338 for kk=1:numel(res)
339 ft = miir(res(kk), [ 1 -poles(kk)], fs);
340 ft.setHistout(Zi2(kk));
341 pfilts2 = [pfilts2 ft];
342 end
343
344 % --- filter 3 ---
345 res = mres13(:,2);
346 poles = mpoles13(:,2);
347 % construct a struct array of miir filters vectors
348 pfilts3 = [];
349 for kk=1:numel(res)
350 ft = miir(res(kk), [ 1 -poles(kk)], fs);
351 ft.setHistout(Zi3(kk));
352 pfilts3 = [pfilts3 ft];
353 end
354
355 % --- filter 4 ---
356 res = mres24(:,2);
357 poles = mpoles24(:,2);
358 % construct a struct array of miir filters vectors
359 pfilts4 = [];
360 for kk=1:numel(res)
361 ft = miir(res(kk), [ 1 -poles(kk)], fs);
362 ft.setHistout(Zi4(kk));
363 pfilts4 = [pfilts4 ft];
364 end
365
366 % switch between output options
367 switch lower(setpar)
368 case 'default'
369 % ----------------------------------
370 % 3) Filtering data
371
372 for jj = 1:2:numel(bs)-1-odc
373
374 % add yunits, taking them from plist or, if empty, from input objects
375 Iunits1 = bs(jj).yunits;
376 Iunits2 = bs(jj+1).yunits;
377 Ounits = find(pl, 'yunits');
378 switch class(Ounits)
379 case 'cell'
380 Ounits1 = Ounits{1};
381 Ounits2 = Ounits{2};
382 case 'unit'
383 Ounits1 = Ounits(1);
384 Ounits2 = Ounits(2);
385 otherwise
386 error('### Bad format for unit container. Supporting vector or cell array');
387 end
388 if eq(unit(Ounits1), unit('')) && eq(unit(Ounits2), unit(''))
389 Ounits1 = bs(jj).yunits;
390 Ounits2 = bs(jj+1).yunits;
391 end
392
393 pfilts1.setIunits(Iunits1);
394 pfilts1.setOunits(Ounits1);
395 pfilts2.setIunits(Iunits2);
396 pfilts2.setOunits(Ounits1);
397 pfilts3.setIunits(Iunits1);
398 pfilts3.setOunits(Ounits2);
399 pfilts4.setIunits(Iunits2);
400 pfilts4.setOunits(Ounits2);
401
402 bs(jj) = filter(cs(jj), pfilts1) + filter(cs(jj+1), pfilts2);
403 bs(jj+1) = filter(cs(jj), pfilts3) + filter(cs(jj+1), pfilts4);
404
405 % -----------------------------------
406 % 4) Output data
407 % name for this objects
408 bs(jj).name = sprintf('noisegen2D(%s)_c1', [ao_invars{jj} ao_invars{jj+1}]);
409 bs(jj+1).name = sprintf('noisegen2D(%s)_c2', [ao_invars{jj} ao_invars{jj+1}]);
410 % Collect the filters into procinfo
411 bs(jj).procinfo = plist('Filt11', pfilts1,'Filt12', pfilts2);
412 bs(jj+1).procinfo = plist('Filt21', pfilts3,'Filt22', pfilts4);
413 % add history
414 bs(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
415 bs(jj+1).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
416 end
417
418
419 % Set output
420 if nargout == numel(bs)
421 % List of outputs
422 for ii = 1:numel(bs)
423 varargout{ii} = bs.index(ii);
424 end
425 else
426 % Single output
427 varargout{1} = bs;
428 end
429
430 case 'filter'
431
432 pfilts1.setIunits(Iunits);
433 pfilts1.setOunits(Ounits);
434 pfilts2.setIunits(Iunits);
435 pfilts2.setOunits(Ounits);
436 pfilts3.setIunits(Iunits);
437 pfilts3.setOunits(Ounits);
438 pfilts4.setIunits(Iunits);
439 pfilts4.setOunits(Ounits);
440
441 fil11 = filterbank(plist('filters',pfilts1,'type','parallel'));
442 fil11.setName('filter 11');
443 fil12 = filterbank(plist('filters',pfilts2,'type','parallel'));
444 fil12.setName('filter 12');
445 fil21 = filterbank(plist('filters',pfilts3,'type','parallel'));
446 fil21.setName('filter 21');
447 fil22 = filterbank(plist('filters',pfilts4,'type','parallel'));
448 fil22.setName('filter 22');
449
450 fil = matrix(plist('objs',[fil11,fil21,fil12,fil22],'shape',[2 2]));
451 fil.setName(sprintf('noisegen2D(%s)',ao_invars{:}));
452 fil.addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
453
454 % Single output
455 varargout{1} = fil;
456
457 end
458
459 end
460
461
462 %--------------------------------------------------------------------------
463 % Get Info Object
464 %--------------------------------------------------------------------------
465 function ii = getInfo(varargin)
466 if nargin == 1 && strcmpi(varargin{1}, 'None')
467 sets = {};
468 pl = [];
469 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
470 sets{1} = varargin{1};
471 pl = getDefaultPlist(sets{1});
472 else
473 sets = SETS();
474 % get plists
475 pl(size(sets)) = plist;
476 for kk = 1:numel(sets)
477 pl(kk) = getDefaultPlist(sets{kk});
478 end
479 end
480 % Build info object
481 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: noisegen2D.m,v 1.34 2011/04/08 08:56:13 hewitson Exp $', sets, pl);
482 ii.setArgsmin(2);
483 ii.setOutmin(1);
484 ii.setModifier(false);
485 end
486
487
488 %--------------------------------------------------------------------------
489 % Defintion of Sets
490 %--------------------------------------------------------------------------
491
492 function out = SETS()
493 out = {...
494 'Default', ...
495 'Filter' ...
496 };
497 end
498
499 %--------------------------------------------------------------------------
500 % Get Default Plist
501 %--------------------------------------------------------------------------
502 function plout = getDefaultPlist(set)
503 persistent pl;
504 persistent lastset;
505 if exist('pl', 'var')==0 || isempty(pl) || ~strcmp(lastset, set)
506 pl = buildplist(set);
507 lastset = set;
508 end
509 plout = pl;
510 end
511
512 function pl = buildplist(set)
513
514 pl = plist();
515
516 switch lower(set)
517 case 'default'
518 % CSD11
519 p = param({'csd11', 'A frequency-series AO describing the model csd11'}, paramValue.EMPTY_DOUBLE);
520 pl.append(p);
521
522 % CSD12
523 p = param({'csd12', 'A frequency-series AO describing the model csd12'}, paramValue.EMPTY_DOUBLE);
524 pl.append(p);
525
526 % CSD21
527 p = param({'csd21', 'A frequency-series AO describing the model csd21'}, paramValue.EMPTY_DOUBLE);
528 pl.append(p);
529
530 % CSD22
531 p = param({'csd22', 'A frequency-series AO describing the model csd22'}, paramValue.EMPTY_DOUBLE);
532 pl.append(p);
533
534 % Yunits
535 p = param({'yunits',['Unit on Y axis. <br>' ...
536 'If left empty, it will take the y-units from the input objects']}, {'',''});
537 pl.append(p);
538
539 case 'filter'
540
541 % Fs
542 p = param({'fs','The sampling frequency to design for.'}, paramValue.DOUBLE_VALUE(1));
543 pl.append(p);
544
545 % Iunits
546 p = param({'iunits','The input units of the filter.'}, paramValue.EMPTY_STRING);
547 pl.append(p);
548
549 % Ounits
550 p = param({'ounits','The output units of the filter.'}, paramValue.EMPTY_STRING);
551 pl.append(p);
552
553 end
554
555 % MaxIter
556 p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
557 pl.append(p);
558
559 % PoleType
560 p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
561 '<li>use real starting poles</li>' ...
562 '<li>generates complex conjugate poles of the<br>'...
563 'type <tt>a.*exp(theta*pi*j)</tt><br>'...
564 'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
565 '<li>generates complex conjugate poles of the type<br>'...
566 '<tt>a.*exp(theta*pi*j)</tt><br>'...
567 'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {3, {1,2,3}, paramValue.SINGLE});
568 pl.append(p);
569
570 % MinOrder
571 p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
572 pl.append(p);
573
574 % MaxOrder
575 p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
576 pl.append(p);
577
578 % Weights
579 p = param({'Weights', ['Choose weighting for the fit:<ol>'...
580 '<li>equal weights for each point</li>'...
581 '<li>weight with <tt>1/abs(model)</tt></li>'...
582 '<li>weight with <tt>1/abs(model).^2</tt></li>'...
583 '<li>weight with inverse of the square mean spread<br>'...
584 'of the model</li></ol>']}, paramValue.DOUBLE_VALUE(3));
585 pl.append(p);
586
587 % Plot
588 p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);
589 pl.append(p);
590
591 % Disp
592 p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);
593 pl.append(p);
594
595 % MSEVARTOL
596 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
597 'realtive variation of the mean squared error is<br>'...
598 'smaller than the value specified. This<br>'...
599 'option is useful for finding the minimum of Chi-squared.']}, ...
600 paramValue.DOUBLE_VALUE(1e-2));
601 pl.append(p);
602
603 % FITTOL
604 p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
605 'squared error value is lower than the value<br>'...
606 'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
607 pl.append(p);
608
609 % UseSym
610 p = param({'UseSym', ['Use symbolic calculation in eigen-decomposition.<ul>'...
611 '<li>0 - perform double-precision calculation in the<br>'...
612 'eigendecomposition procedure to identify 2-Dim<br>'...
613 'systems and for poles stabilization</li>'...
614 '<li>1 - uses symbolic math toolbox variable precision<br>'...
615 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
616 'system identification and double-precison for<br>'...
617 'poles stabilization</li>'...
618 '<li>2 - uses symbolic math toolbox variable precision<br>'...
619 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
620 'system identification and for poles stabilization.']}, {1, {0,1,2}, paramValue.SINGLE});
621 pl.append(p);
622
623
624 end
625
626
627