Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/fromCSD.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 % FROMCSD Construct a matrix filter from cross-spectral density matrix
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % FUNCTION: fromCSD
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 MultiChannelNoise constructor of the matrix class.
9 %
10 %
11 % CALL: a = fromCSD(a, pl)
12 %
13 % PARAMETER: pl: plist containing 'csd'
14 %
15 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16 function a = fromCSD(a, pli)
17
18 VERSION = '$Id: fromCSD.m,v 1.8 2010/10/29 16:09:12 ingo Exp $';
19
20 % get AO info
21 iobj = matrix.getInfo('matrix', 'From CSD');
22
23 % Set the method version string in the minfo object
24 iobj.setMversion([VERSION '-->' iobj.mversion]);
25
26 % Add default values
27 pl = parse(pli, iobj.plists);
28
29 % Get parameters and set params for fit
30 csdm = find(pl, 'csd');
31 if isa(csdm,'matrix') % get elements out of the input matrix
32 csdao = csdm.objs;
33 else
34 csdao = csdm;
35 end
36 fs = find(pl, 'fs');
37
38 target = lower(find(pl, 'targetobj')); % decide to perform s domain or z domain identification
39 % if target is parfrac output a matrix of parfarc objs (s domain
40 % identification)
41 % if target is miir output a matrix of filterbank parallel miir objects
42 % (z domain identification)
43
44 usesym = lower(find(pl, 'UseSym'));
45
46 if (fs == 0) && strcmpi(target,'miir')
47 error('### Please provide a valid sampling frequency for CSD constructor.');
48 elseif isempty(fs) && strcmpi(target,'miir')
49 error('### Please provide a valid sampling frequency for CSD constructor.');
50 end
51
52 % get units for filters
53 tgiunit = find(pl,'iunits');
54 tgounit = find(pl,'ounits');
55
56 params = struct();
57
58 params.Nmaxiter = find(pl, 'MaxIter');
59 params.minorder = find(pl, 'MinOrder');
60 params.maxorder = find(pl, 'MaxOrder');
61 params.spolesopt = find(pl, 'PoleType');
62 params.weightparam = find(pl, 'Weights');
63
64 % set the target output
65 if strcmpi(target,'miir')
66 params.TargetDomain = 'z';
67 elseif strcmpi(target,'parfrac')
68 params.TargetDomain = 's';
69 else
70 error('### Unknown option for ''targetobj''.');
71 end
72
73 % Tolerance for MSE Value
74 lrscond = find(pl, 'FITTOL');
75 % give an error for strange values of lrscond
76 if lrscond<0
77 error('### Negative values for FITTOL are not allowed. ')
78 end
79 % handling data
80 lrscond = -1*log10(lrscond);
81 % give a warning for strange values of lrscond
82 if lrscond<0
83 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
84 end
85 params.lrscond = lrscond;
86
87 % Tolerance for the MSE relative variation
88 msevar = find(pl, 'MSEVARTOL');
89 % handling data
90 msevar = -1*log10(msevar);
91 % give a warning for strange values of msevar
92 if msevar<0
93 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
94 end
95 params.msevar = msevar;
96
97 if isempty(params.msevar)
98 params.ctp = 'chival';
99 else
100 params.ctp = 'chivar';
101 end
102
103 if(find(pl, 'plot'))
104 params.plot = 1;
105 else
106 params.plot = 0;
107 end
108
109 params.fs = fs;
110 params.dterm = 0; % it is better to fit without direct term
111
112 % check if symbolic calculation is required
113 if strcmpi(usesym,'on')
114 params.usesym = 1;
115 elseif strcmpi(usesym,'off')
116 params.usesym = 0;
117 else
118 error('### Unknown option for ''UseSym''.');
119 end
120
121 % extracting csd
122 if numel(csdao)==1 % one dimensional psd
123 csd = csdao.y;
124 freq = csdao.x;
125 dim = 'one';
126 else % multichannel system
127 dim = 'multi';
128 [nn,mm] = size(csdao);
129 if nn~=mm
130 error('### CSD Matrix must be square. ')
131 end
132 freq = csdao.index(1).x;
133 for ii = 1:nn
134 for jj = 1:nn
135 tcsd = csdao.index(ii,jj).y;
136 % willing to work with columns
137 [aa,bb] = size(tcsd);
138 if aa<bb
139 tcsd = tcsd.';
140 end
141 csd(ii,jj,:) = tcsd;
142 end
143 end
144
145 end
146
147
148 % call csd2tf
149 % ostruct is a struct array whose fields contain the residues and poles
150 % of estimated TFs. Since the fit is porformed on the columns of the TF
151 % matrix, each element of the array contains the residues and poles
152 % corresponding to the functions on the given column of the TF matrix.
153
154 %ostruct = utils.math.csd2tf(csd,freq,params);
155
156 ostruct = utils.math.csd2tf2(csd,freq,params);
157
158
159 % the filter for each channel is implemented by the rows of the TF matrix
160
161 switch dim
162 case 'one'
163
164 switch target
165 case 'miir'
166 % --- filter ---
167 res = ostruct.res;
168 poles = ostruct.poles;
169 % construct a struct array of miir filters vectors
170 pfilts(numel(res),1) = miir;
171 for kk=1:numel(res)
172 ft = miir(res(kk), [ 1 -poles(kk)], fs);
173 ft.setIunits(unit(tgiunit));
174 ft.setOunits(unit(tgounit));
175 pfilts(kk,1) = ft;
176 clear ft
177 end
178 filt = filterbank(pfilts,'parallel');
179
180 a = matrix(filt);
181
182 % Add history
183 a.addHistory(iobj, pl, [], [csdao(:).hist]);
184
185 case 'parfrac'
186 res = ostruct.res;
187 poles = ostruct.poles;
188
189 fbk = parfrac(res,poles,0);
190 fbk.setIunits(unit(tgiunit));
191 fbk.setOunits(unit(tgounit));
192
193 a = matrix(fbk);
194
195 % Add history
196 a.addHistory(iobj, pl, [], [csdm(:).hist]);
197 end
198
199
200 case 'multi'
201
202 switch target
203 case 'miir'
204 % init filters array
205 %fbk(nn*nn,1) = filterbank;
206 %fbk = filterbank.newarray([nn nn]);
207
208 for zz=1:nn*nn % run over system dimension
209 % --- get column filter coefficients ---
210 % each column of mres\mpoles are the coefficients of a given filter
211 clear res poles
212 res = ostruct(zz).res;
213 poles = ostruct(zz).poles;
214
215 % construct a struct array of miir filters vectors
216 %ft(numel(res),1) = miir;
217 for kk=1:numel(res)
218 ft(kk,1) = miir(res(kk), [1 -poles(kk)], fs);
219 ft(kk,1).setIunits(unit(tgiunit));
220 ft(kk,1).setOunits(unit(tgounit));
221 end
222
223 fbk(zz,1) = filterbank(ft,'parallel');
224 clear ft
225
226 end
227
228 mfbk = reshape(fbk,nn,nn);
229
230 a = matrix(mfbk);
231
232 % % overide default name
233 % if strcmpi(pl.find('name'), 'None')
234 % pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
235 % end
236 % % overide default description
237 % if isempty(pl.find('description'))
238 % pl.pset('description', csdao.description);
239 % end
240
241 % Add history
242 a.addHistory(iobj, pl, [], [csdao(:).hist]);
243
244 case 'parfrac'
245 % init filters array
246 %fbk(nn*nn,1) = parfrac;
247
248 for zz=1:nn*nn % run over system dimension
249 % --- get column filter coefficients ---
250 % each column of mres\mpoles are the coefficients of a given filter
251 clear res poles
252 res = ostruct(zz).res;
253 poles = ostruct(zz).poles;
254
255 fbk(zz,1) = parfrac(res,poles,0);
256 fbk(zz,1).setIunits(unit(tgiunit));
257 fbk(zz,1).setOunits(unit(tgounit));
258
259 end
260
261 mfbk = reshape(fbk,nn,nn);
262
263 a = matrix(mfbk);
264
265 % % overide default name
266 % if strcmpi(pl.find('name'), 'None')
267 % pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
268 % end
269 % % overide default description
270 % if isempty(pl.find('description'))
271 % pl.pset('description', csdao.description);
272 % end
273
274 % Add history
275 a.addHistory(iobj, pl, [], [csdao(:).hist]);
276 end
277
278 end
279
280 % Set properties from the plist
281 warning('off', utils.const.warnings.METHOD_NOT_FOUND);
282 % remove parameters we already used
283 pl_set = copy(pl,1);
284 if pl_set.isparam('csd')
285 pl_set.remove('csd');
286 end
287 if pl_set.isparam('fs')
288 pl_set.remove('fs');
289 end
290 if pl_set.isparam('targetobj')
291 pl_set.remove('targetobj');
292 end
293 if pl_set.isparam('UseSym')
294 pl_set.remove('UseSym');
295 end
296 if pl_set.isparam('iunits')
297 pl_set.remove('iunits');
298 end
299 if pl_set.isparam('ounits')
300 pl_set.remove('ounits');
301 end
302 if pl_set.isparam('MaxIter')
303 pl_set.remove('MaxIter');
304 end
305 if pl_set.isparam('MinOrder')
306 pl_set.remove('MinOrder');
307 end
308 if pl_set.isparam('MaxOrder')
309 pl_set.remove('MaxOrder');
310 end
311 if pl_set.isparam('PoleType')
312 pl_set.remove('PoleType');
313 end
314 if pl_set.isparam('Weights')
315 pl_set.remove('Weights');
316 end
317 if pl_set.isparam('FITTOL')
318 pl_set.remove('FITTOL');
319 end
320 if pl_set.isparam('MSEVARTOL')
321 pl_set.remove('MSEVARTOL');
322 end
323 if pl_set.isparam('plot')
324 pl_set.remove('plot');
325 end
326 a.setProperties(pl_set);
327 warning('on', utils.const.warnings.METHOD_NOT_FOUND);
328
329
330 end
331
332