Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/whiten2D.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 % WHITEN2D whiten the noise for two cross correlated time series.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: whiten2D whitens cross-correlated time-series. Whitening
5 % filters are constructed by a fitting procedure to the cross-spectrum
6 % models provided.
7 % Note: The function assumes that the input model corresponds to the
8 % one-sided csd of the data to be whitened.
9 %
10 % ALGORITHM:
11 % 1) Fit a set of partial fraction z-domain filters using
12 % utils.math.psd2wf
13 % 2) Convert to bank of mIIR filters
14 % 3) Filter time-series in parallel
15 % The filtering process is:
16 % b(1) = Filt11(a(1)) + Filt12(a(2))
17 % b(2) = Filt21(a(1)) + Filt22(a(2))
18 %
19 % CALL: b = whiten2D(a, pl) % returns whitened time-series AOs
20 % [b1,b2] = whiten2D(a1, a2, pl)
21 % [b1,b2,...,bn] = whiten2D(a1,a2,...,an, pl);
22 % Note: Input AOs must come in couples.
23 % Note: this method cannot be used as a modifier, the
24 % call a.whiten2D(pl) is forbidden.
25 %
26 % INPUT:
27 %
28 % - a is a couple of two colored noise time-series AOs
29 %
30 % OUTPUT:
31 %
32 % - b is a couple of "whitened" time-series AOs. The whitening
33 % filters used are stored in the objects procinfo field under
34 % the parameters:
35 % - b(1): 'Filt11' and 'Filt12'
36 % - b(2): 'Filt21' and 'Filt22'
37 %
38 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'whiten2D')">Parameters Description</a>
39 %
40 % VERSION: $Id: whiten2D.m,v 1.35 2011/04/08 08:56:16 hewitson Exp $
41 %
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43
44 function varargout = whiten2D(varargin)
45
46 % Check if this is a call for parameters
47 if utils.helper.isinfocall(varargin{:})
48 varargout{1} = getInfo(varargin{3});
49 return
50 end
51
52 import utils.const.*
53 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
54
55 % Collect input variable names
56 in_names = cell(size(varargin));
57 for ii = 1:nargin,in_names{ii} = inputname(ii);end
58
59 % Collect all AOs and plists
60 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
61 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
62
63 % Decide on a deep copy or a modify
64 bs = copy(as, nargout);
65 cs = copy(as, nargout);
66 inhists = [as.hist];
67
68 % combine plists
69 pl = parse(pl, getDefaultPlist());
70 pl.getSetRandState();
71
72
73 % Extract necessary model parameters
74 csd11 = find(pl, 'csd11');
75 csd12 = find(pl, 'csd12');
76 csd21 = find(pl, 'csd21');
77 csd22 = find(pl, 'csd22');
78
79 if nargout == 0
80 error('### noisegen2D cannot be used as a modifier. Please give an output variable.');
81 end
82
83 % Check the number of input AO
84 if numel(bs)==1
85 error('!!! One input AO! The input AOs must come in pairs!');
86 end
87 % Checks that the input AOs come in pairs
88 odc = 0;
89 if rem(numel(as),2)
90 warning('The input AOs must come in pairs! Skipping AO %s during calculation', ao_invars{end});
91 odc = 1;
92 end
93
94 % Loop over input AOs to check for non time series objects
95 fsv = zeros(numel(bs),1);
96 for jj=1:numel(bs)
97 if ~isa(bs(jj).data, 'tsdata')
98 error('!!! %s expects ao/tsdata objects. ', mfilename);
99 end
100 fsv(jj,1) = bs(jj).data.fs; % collecting sampling frequencies
101 end
102 % Check that input Aos have the same sampling frequency
103 if any(diff(fsv))
104 error('!!! Sampling frequency must be the same for all input objects')
105 end
106
107 % ------------------- Coloring Noise
108
109 % ----------------------------------
110 % 1) - Fitting the models to identify the innovation filters
111
112 % Build input structure for psd2tf
113 params = struct();
114
115 params.idtp = 1;
116 params.Nmaxiter = find(pl, 'MaxIter');
117 params.minorder = find(pl, 'MinOrder');
118 params.maxorder = find(pl, 'MaxOrder');
119 params.spolesopt = find(pl, 'PoleType');
120 params.weightparam = find(pl, 'Weights');
121 params.usesym = find(pl, 'UseSym');
122 params.spy = find(pl, 'Disp');
123 params.keepvar = find(pl, 'keepvar');
124
125 % Tolerance for MSE Value
126 lrscond = find(pl, 'FITTOL');
127 % give an error for strange values of lrscond
128 if lrscond<0
129 error('!!! Negative values for FITTOL are not allowed !!!')
130 end
131 % handling data
132 lrscond = -1*log10(lrscond);
133 % give a warning for strange values of lrscond
134 if lrscond<0
135 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
136 end
137 params.lrscond = lrscond;
138
139 % Tolerance for the MSE relative variation
140 msevar = find(pl, 'MSEVARTOL');
141 % handling data
142 msevar = -1*log10(msevar);
143 % give a warning for strange values of msevar
144 if msevar<0
145 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
146 end
147 params.msevar = msevar;
148
149 if isempty(params.msevar)
150 params.ctp = 'chival';
151 else
152 params.ctp = 'chivar';
153 end
154
155 if(find(pl, 'plot'))
156 params.plot = 1;
157 else
158 params.plot = 0;
159 end
160
161 params.fs = fsv(1,1);
162 params.dterm = 0;
163
164 % get data variance
165 vars = [1 1];
166 if numel(bs)==2
167 for ii=1:numel(bs)
168 b = bs(ii).y;
169 v = var(b);
170 vars(ii) = v;
171 end
172 end
173 params.vars = vars;
174
175 % call psd2wf
176 ostruct = utils.math.psd2wf(csd11.y,csd12.y,csd21.y,csd22.y,csd11.x,params);
177
178 % ----------------------------------
179 % 2) - Convert into MIIR filters
180
181 fs = fsv(1,1);
182
183 % --- filter 1 ---
184 res = ostruct(1).res;
185 poles = ostruct(1).poles;
186 dterm = ostruct(1).dterm;
187 % construct a struct array of miir filters vectors
188 pfilts1 = [];
189 for kk=1:numel(res)
190 ft = miir(res(kk), [ 1 -poles(kk)], fs);
191 pfilts1 = [pfilts1 ft];
192 end
193 % pfilts1 = [pfilts1 miir(dterm, [1 0], fs)];
194
195 % --- filter 2 ---
196 res = ostruct(2).res;
197 poles = ostruct(2).poles;
198 dterm = ostruct(2).dterm;
199 % construct a struct array of miir filters vectors
200 pfilts2 = [];
201 for kk=1:numel(res)
202 ft = miir(res(kk), [ 1 -poles(kk)], fs);
203 pfilts2 = [pfilts2 ft];
204 end
205 % pfilts2 = [pfilts2 miir(dterm, [1 0], fs)];
206
207 % --- filter 3 ---
208 res = ostruct(3).res;
209 poles = ostruct(3).poles;
210 dterm = ostruct(3).dterm;
211 % construct a struct array of miir filters vectors
212 pfilts3 = [];
213 for kk=1:numel(res)
214 ft = miir(res(kk), [ 1 -poles(kk)], fs);
215 pfilts3 = [pfilts3 ft];
216 end
217 % pfilts3 = [pfilts3 miir(dterm, [1 0], fs)];
218
219 % --- filter 4 ---
220 res = ostruct(4).res;
221 poles = ostruct(4).poles;
222 dterm = ostruct(4).dterm;
223 % construct a struct array of miir filters vectors
224 pfilts4 = [];
225 for kk=1:numel(res)
226 ft = miir(res(kk), [ 1 -poles(kk)], fs);
227 pfilts4 = [pfilts4 ft];
228 end
229 % pfilts4 = [pfilts4 miir(dterm, [1 0], fs)];
230
231 % ----------------------------------
232 % 3) Filtering data
233
234 for jj = 1:2:numel(bs)-1-odc
235
236 bs(jj) = filter(cs(jj), pfilts1) + filter(cs(jj+1), pfilts2);
237 bs(jj+1) = filter(cs(jj), pfilts3) + filter(cs(jj+1), pfilts4);
238
239 % -----------------------------------
240 % 4) Output data
241
242 % name for this object
243 bs(jj).name = sprintf('whiten2D(%s)_c1', [ao_invars{jj} ao_invars{jj+1}]);
244 % Collect the filters into procinfo
245 bs(jj).procinfo = combine(plist('Filt11', pfilts1,'Filt12', pfilts2),as(jj).procinfo);
246 % add history
247 bs(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
248
249 % name for this object
250 bs(jj+1).name = sprintf('whiten2D(%s)_c2', [ao_invars{jj} ao_invars{jj+1}]);
251 % Collect the filters into procinfo
252 bs(jj+1).procinfo = combine(plist('Filt21', pfilts3,'Filt22', pfilts4),as(jj+1).procinfo);
253 % add history
254 bs(jj+1).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
255
256 end
257
258 % clear errors
259 bs.clearErrors;
260
261 % Set output
262 if nargout == numel(bs)
263 % List of outputs
264 for ii = 1:numel(bs)
265 varargout{ii} = bs.index(ii);
266 end
267 else
268 % Single output
269 varargout{1} = bs;
270 end
271 end
272
273 %--------------------------------------------------------------------------
274 % Get Info Object
275 %--------------------------------------------------------------------------
276 function ii = getInfo(varargin)
277 if nargin == 1 && strcmpi(varargin{1}, 'None')
278 sets = {};
279 pl = [];
280 else
281 sets = {'Default'};
282 pl = getDefaultPlist;
283 end
284 % Build info object
285 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: whiten2D.m,v 1.35 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
286 ii.setArgsmin(2);
287 ii.setOutmin(2);
288 ii.setModifier(false);
289 end
290
291 %--------------------------------------------------------------------------
292 % Get Default Plist
293 %--------------------------------------------------------------------------
294
295 function plout = getDefaultPlist()
296 persistent pl;
297 if exist('pl', 'var')==0 || isempty(pl)
298 pl = buildplist();
299 end
300 plout = pl;
301 end
302
303 function pl = buildplist()
304
305 pl = plist();
306
307 % CSD11
308 p = param({'csd11', 'A frequency-series AO describing the model csd11'}, paramValue.EMPTY_DOUBLE);
309 pl.append(p);
310
311 % CSD12
312 p = param({'csd12', 'A frequency-series AO describing the model csd12'}, paramValue.EMPTY_DOUBLE);
313 pl.append(p);
314
315 % CSD21
316 p = param({'csd21', 'A frequency-series AO describing the model csd21'}, paramValue.EMPTY_DOUBLE);
317 pl.append(p);
318
319 % CSD22
320 p = param({'csd22', 'A frequency-series AO describing the model csd22'}, paramValue.EMPTY_DOUBLE);
321 pl.append(p);
322
323 % MaxIter
324 p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
325 pl.append(p);
326
327 % PoleType
328 p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
329 '<li>use real starting poles</li>' ...
330 '<li>generates complex conjugate poles of the<br>'...
331 'type <tt>a.*exp(theta*pi*j)</tt><br>'...
332 'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
333 '<li>generates complex conjugate poles of the type<br>'...
334 '<tt>a.*exp(theta*pi*j)</tt><br>'...
335 'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {3, {1,2,3}, paramValue.SINGLE});
336 pl.append(p);
337
338 % MinOrder
339 p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
340 pl.append(p);
341
342 % MaxOrder
343 p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
344 pl.append(p);
345
346 % Weights
347 p = param({'Weights', ['Choose weighting for the fit:<ol>'...
348 '<li>equal weights for each point</li>'...
349 '<li>weight with <tt>1/abs(model)</tt></li>'...
350 '<li>weight with <tt>1/abs(model).^2</tt></li>'...
351 '<li>weight with inverse of the square mean spread<br>'...
352 'of the model</li></ol>']}, paramValue.DOUBLE_VALUE(3));
353 pl.append(p);
354
355 % Plot
356 p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);
357 pl.append(p);
358
359 % Disp
360 p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);
361 pl.append(p);
362
363 % MSEVARTOL
364 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
365 'realtive variation of the mean squared error is<br>'...
366 'smaller than the value specified. This<br>'...
367 'option is useful for finding the minimum of Chi-squared.']}, ...
368 paramValue.DOUBLE_VALUE(1e-2));
369 pl.append(p);
370
371 % FITTOL
372 p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
373 'squared error value is lower than the value<br>'...
374 'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
375 pl.append(p);
376
377 % UseSym
378 p = param({'UseSym', ['Use symbolic calculation in eigen-decomposition.<ul>'...
379 '<li>0 - perform double-precision calculation in the<br>'...
380 'eigendecomposition procedure to identify 2-Dim<br>'...
381 'systems and for poles stabilization</li>'...
382 '<li>1 - uses symbolic math toolbox variable precision<br>'...
383 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
384 'system identification and double-precison for<br>'...
385 'poles stabilization</li>'...
386 '<li>2 - uses symbolic math toolbox variable precision<br>'...
387 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
388 'system identification and for poles stabilization.']}, {1, {0,1,2}, paramValue.SINGLE});
389 pl.append(p);
390
391 % Keep var
392 p = param({'keepvar', '???'}, paramValue.TRUE_FALSE);
393 p.val.setValIndex(2);
394 pl.append(p);
395
396 end
397
398
399 % PARAMETERS:
400 %
401 % 'csd11' - a frequency-series AO describing the model csd11
402 % 'csd12' - a frequency-series AO describing the model csd12
403 % 'csd21' - a frequency-series AO describing the model csd21
404 % 'csd22' - a frequency-series AO describing the model csd22
405 %
406 % 'MaxIter' - Maximum number of iterations in fit routine
407 % [default: 30]
408 %
409 % 'PoleType' - Choose the pole type for fitting:
410 % 1 - use real starting poles
411 % 2 - generates complex conjugate poles of the
412 % type a.*exp(theta*pi*j)
413 % with theta = linspace(0,pi,N/2+1)
414 % 3 - generates complex conjugate poles of the type
415 % a.*exp(theta*pi*j)
416 % with theta = linspace(0,pi,N/2+2) [default]
417 %
418 % 'MinOrder' - Minimum order to fit with. [default: 2]
419 %
420 % 'MaxOrder' - Maximum order to fit with. [default: 25]
421 %
422 % 'Weights' - choose weighting for the fit: [default: 2]
423 % 1 - equal weights for each point
424 % 2 - weight with 1/abs(model)
425 % 3 - weight with 1/sqrt(abs(model))
426 % 4 - weight with inverse of the square mean spread
427 % of the model
428 %
429 % 'Plot' - plot results of each fitting step. [default: false]
430 %
431 % 'Disp' - Display the progress of the fitting iteration.
432 % [default: false]
433 %
434 % 'FITTOL' - Mean Squared Error Value - Check if the mean
435 % squared error value is lower than the value
436 % specified. [defalut: 1e-2]
437 %
438 % 'MSEVARTOL' - Mean Squared Error Variation - Check if the
439 % realtive variation of the mean squared error is
440 % smaller than the value specified. This
441 % option is useful for finding the minimum of Chi
442 % squared. [default: 1e-2]]
443 %
444 % 'UseSym' - Use symbolic calculation in eigendecomposition.
445 % [default: 0]
446 % 0 - perform double-precision calculation in the
447 % eigendecomposition procedure to identify 2dim
448 % systems and for poles stabilization
449 % 1 - uses symbolic math toolbox variable precision
450 % arithmetic in the eigendecomposition for 2dim
451 % system identification and double-precison for
452 % poles stabilization
453 % 2 - uses symbolic math toolbox variable precision
454 % arithmetic in the eigendecomposition for 2dim
455 % system identification and for poles
456 % stabilization