Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/noisegen1D.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 % NOISEGEN1D generates colored noise from white noise.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: noisegen1D can work in two different modes:
5 %
6 % ------------------------------------------------------------------------
7 % 1) Generates colored noise from white noise with a given spectrum. The
8 % function constructs a coloring filter through a fitting procedure to the
9 % model provided. If no model is provided an error is prompted. The colored
10 % noise provided has one-sided psd corresponding to the input model.
11 %
12 % This mode correspond to the 'Default' set for the method (see the list of
13 % parameters).
14 %
15 % ALGORITHM:
16 % 1) Fit a set of partial fraction z-domain filters using
17 % utils.math.psd2tf
18 % 2) Convert to array of MIIR filters
19 % 3) Filter time-series in parallel
20 %
21 % CALL: b = noisegen1D(a, pl)
22 %
23 % INPUT:
24 % - a is a white noise time-series analysis object or a
25 % vector of analysis objects
26 % - pl is a plist with the input parameters
27 %
28 % OUTPUT:
29 %
30 % - b Colored time-series AOs. The coloring filters used
31 % are stored in the objects procinfo field under the
32 % parameter 'Filt'.
33 %
34 % ------------------------------------------------------------------------
35 %
36 % 2) Generates noise coloring filters for given input psd models.
37 %
38 % This mode correspond to the 'Filter' set for the method (see the list of
39 % parameters).
40 %
41 % ALGORITHM:
42 % 1) Fit a set of partial fraction z-domain filters
43 % 2) Convert to array of MIIR filters
44 %
45 % CALL: fil = noisegen1D(psd, pl)
46 %
47 % INPUT:
48 % - psd is a fsdata analysis object representing the desired
49 % model for the power spectral density of the colored noise
50 % - pl is a plist with the input parameters
51 %
52 % OUTPUT:
53 %
54 % - fil is a filterbank parallel object which elements are
55 % miir filters. Filters are initialized to avoid startup
56 % transients.
57 %
58 % ------------------------------------------------------------------------
59 %
60 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'noisegen1D')">Parameters Description</a>
61 %
62 % VERSION: $Id: noisegen1D.m,v 1.30 2011/04/08 08:56:15 hewitson Exp $
63 %
64 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65
66 function varargout = noisegen1D(varargin)
67
68 % Check if this is a call for parameters
69 if utils.helper.isinfocall(varargin{:})
70 varargout{1} = getInfo(varargin{3});
71 return
72 end
73
74 import utils.const.*
75 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
76
77 % Collect input variable names
78 in_names = cell(size(varargin));
79 for ii = 1:nargin,in_names{ii} = inputname(ii);end
80
81 % Collect all AOs and plists
82 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
83 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
84
85 % Decide on a deep copy or a modify
86 bs = copy(as, nargout);
87 inhists = [as.hist];
88
89 % combine plists
90 % define input/output combinations. Different combination are
91 % 1) input tsdata and csd# into the plist, output are colored tsdata
92 % 2) input fsdata, output is a coloring filter (in a matrix)
93 if isempty(pl) % no model input, get model from input
94 setpar = 'Filter';
95 elseif isa(as(1).data,'fsdata') % get model from input
96 setpar = 'Filter';
97 else % get model from plist, output tsdata, back compatibility mode
98 setpar = 'Default';
99 end
100
101 pl = parse(pl, getDefaultPlist(setpar));
102 pl.getSetRandState();
103
104 switch lower(setpar)
105 case 'default'
106 % Extract necessary parameters
107 model = find(pl, 'model');
108 case 'filter'
109 fs = find(pl,'fs');
110 Iunits = find(pl, 'Iunits');
111 Ounits = find(pl, 'Ounits');
112
113 % init filter objects
114 fil = filterbank.initObjectWithSize(size(bs,1),size(bs,2));
115 end
116
117 % start building input structure for psd2tf
118 params.idtp = 1;
119 params.Nmaxiter = find(pl, 'MaxIter');
120 params.minorder = find(pl, 'MinOrder');
121 params.maxorder = find(pl, 'MaxOrder');
122 params.spolesopt = find(pl, 'PoleType');
123 params.weightparam = find(pl, 'Weights');
124 params.spy = find(pl, 'Disp');
125
126 % Tolerance for MSE Value
127 lrscond = find(pl, 'FITTOL');
128 % give an error for strange values of lrscond
129 if lrscond<0
130 error('!!! Negative values for FITTOL are not allowed !!!')
131 end
132 % handling data
133 lrscond = -1*log10(lrscond);
134 % give a warning for strange values of lrscond
135 if lrscond<0
136 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
137 end
138 params.lrscond = lrscond;
139
140 % Tolerance for the MSE relative variation
141 msevar = find(pl, 'MSEVARTOL');
142 % handling data
143 msevar = -1*log10(msevar);
144 % give a warning for strange values of msevar
145 if msevar<0
146 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
147 end
148 params.msevar = msevar;
149
150 if isempty(params.msevar)
151 params.ctp = 'chival';
152 else
153 params.ctp = 'chivar';
154 end
155
156 if(find(pl, 'plot'))
157 params.plot = 1;
158 else
159 params.plot = 0;
160 end
161
162 params.usesym = 0;
163 params.dterm = 0; % it is better to fit without dterm
164
165 % Loop over input AOs
166 for jj=1:numel(bs)
167
168 switch lower(setpar)
169 case 'default'
170 if ~isa(bs(jj).data, 'tsdata')
171 warning('!!! %s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
172 else
173 %-------------- Colored noise from white noise
174
175 % 1) If we have no model gives an error
176 if(isempty(model))
177 error('!!! Input a model for the desired PSD')
178 end
179
180 % 2) Noise Generation
181
182 params.fs = bs(jj).fs;
183
184 % call psd2tf
185 [res, poles, dterm, mresp, rdl] = ...
186 utils.math.psd2tf(model.y,[],[],[],model.x,params);
187
188 % get init state
189 Zi = utils.math.getinitstate(res,poles,1,'mtd','svd');
190
191 % 3) Convert to MIIR filters
192 % add yunits, taking them from plist or, if empty, from input object
193 Iunits = as(jj).yunits;
194 Ounits = find(pl, 'yunits');
195 if eq(unit(Ounits), unit(''))
196 Ounits = as(jj).yunits;
197 end
198 pfilts = [];
199 for kk=1:numel(res)
200 ft = miir(res(kk), [ 1 -poles(kk)], bs(jj).fs);
201 ft.setIunits(Iunits);
202 ft.setOunits(Ounits);
203 ft.setHistout(Zi(kk));
204 % build parallel bank of filters
205 pfilts = [pfilts ft];
206 end
207
208 % 4) Filter data
209 bs(jj).filter(pfilts);
210
211 % 5) Output data
212 % add history
213 bs(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
214 % name for this object
215 bs(jj).setName(sprintf('noisegen1D(%s)', ao_invars{jj}));
216 % Collect the filters into procinfo
217 bs(jj).procinfo = plist('filter', pfilts);
218
219 end
220 case 'filter'
221
222 params.fs = fs;
223
224 % call psd2tf
225 [res, poles, dterm, mresp, rdl] = ...
226 utils.math.psd2tf(bs(jj).y,[],[],[],bs(jj).x,params);
227
228 % get init state
229 Zi = utils.math.getinitstate(res,poles,1,'mtd','svd');
230
231 % build miir
232 pfilts = [];
233 for kk=1:numel(res)
234 ft = miir(res(kk), [ 1 -poles(kk)], fs);
235 ft.setIunits(Iunits);
236 ft.setOunits(Ounits);
237 ft.setHistout(Zi(kk));
238 % build parallel bank of filters
239 pfilts = [pfilts ft];
240 end
241
242 % build output filterbanks
243 fil(jj) = filterbank(plist('filters',pfilts,'type','parallel'));
244 fil(jj).setName(sprintf('noisegen1D(%s)',ao_invars{jj}));
245 fil(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
246 end
247 end
248
249 % switch output
250 switch lower(setpar)
251
252 case 'default'
253 % Set output
254 if nargout > 1 && nargout == numel(bs)
255 % List of outputs
256 for ii = 1:numel(bs)
257 varargout{ii} = bs.index(ii);
258 end
259 else
260 % Single output
261 varargout{1} = bs;
262 end
263
264 case 'filter'
265
266 % Set output
267 if nargout > 1 && nargout == numel(bs)
268 % List of outputs
269 for ii = 1:numel(fil)
270 varargout{ii} = fil.index(ii);
271 end
272 else
273 % Single output
274 varargout{1} = fil;
275 end
276
277 end
278 end
279
280 %--------------------------------------------------------------------------
281 % Get Info Object
282 %--------------------------------------------------------------------------
283 function ii = getInfo(varargin)
284 if nargin == 1 && strcmpi(varargin{1}, 'None')
285 sets = {};
286 pl = [];
287 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
288 sets{1} = varargin{1};
289 pl = getDefaultPlist(sets{1});
290 else
291 sets = SETS();
292 % get plists
293 pl(size(sets)) = plist;
294 for kk = 1:numel(sets)
295 pl(kk) = getDefaultPlist(sets{kk});
296 end
297 end
298 % Build info object
299 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: noisegen1D.m,v 1.30 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
300
301 end
302
303
304 %--------------------------------------------------------------------------
305 % Defintion of Sets
306 %--------------------------------------------------------------------------
307
308 function out = SETS()
309 out = {...
310 'Default', ...
311 'Filter' ...
312 };
313 end
314
315 %--------------------------------------------------------------------------
316 % Get Default Plist
317 %--------------------------------------------------------------------------
318 function plout = getDefaultPlist(set)
319 persistent pl;
320 persistent lastset;
321 if exist('pl', 'var')==0 || isempty(pl) || ~strcmp(lastset, set)
322 pl = buildplist(set);
323 lastset = set;
324 end
325 plout = pl;
326 end
327
328 function pl = buildplist(set)
329
330 pl = plist();
331
332 switch lower(set)
333 case 'default'
334 % Yunits
335 p = param({'yunits',['Unit on Y axis. <br>' ...
336 'If left empty, it will take the y-units from the input object']}, paramValue.STRING_VALUE(''));
337 pl.append(p);
338
339 % Model
340 p = param({'model', 'A frequency-series AO describing the model psd'}, paramValue.EMPTY_DOUBLE);
341 pl.append(p);
342
343 case 'filter'
344 % Fs
345 p = param({'fs','The sampling frequency to design for.'}, paramValue.DOUBLE_VALUE(1));
346 pl.append(p);
347
348 % Iunits
349 p = param({'iunits','The input units of the filter.'}, paramValue.EMPTY_STRING);
350 pl.append(p);
351
352 % Ounits
353 p = param({'ounits','The output units of the filter.'}, paramValue.EMPTY_STRING);
354 pl.append(p);
355 end
356
357 % MaxIter
358 p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
359 pl.append(p);
360
361 % PoleType
362 p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
363 '<li>use real starting poles</li>' ...
364 '<li>generates complex conjugate poles of the<br>'...
365 'type <tt>a.*exp(theta*pi*j)</tt><br>'...
366 'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
367 '<li>generates complex conjugate poles of the type<br>'...
368 '<tt>a.*exp(theta*pi*j)</tt><br>'...
369 'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {3, {1,2,3}, paramValue.SINGLE});
370 pl.append(p);
371
372 % MinOrder
373 p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
374 pl.append(p);
375
376 % MaxOrder
377 p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
378 pl.append(p);
379
380 % Weights
381 p = param({'Weights', ['Choose weighting for the fit:<ol>'...
382 '<li>equal weights for each point</li>'...
383 '<li>weight with <tt>1/abs(model)</tt></li>'...
384 '<li>weight with <tt>1/abs(model).^2</tt></li>'...
385 '<li>weight with inverse of the square mean spread<br>'...
386 'of the model</li></ol>']}, paramValue.DOUBLE_VALUE(3));
387 pl.append(p);
388
389 % Plot
390 p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);
391 pl.append(p);
392
393 % Disp
394 p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);
395 pl.append(p);
396
397 % MSEVARTOL
398 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
399 'realtive variation of the mean squared error is<br>'...
400 'smaller than the value specified. This<br>'...
401 'option is useful for finding the minimum of Chi-squared.']}, ...
402 paramValue.DOUBLE_VALUE(1e-2));
403 pl.append(p);
404
405 % FITTOL
406 p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
407 'squared error value is lower than the value<br>'...
408 'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
409 pl.append(p);
410
411 end
412
413