Mercurial > hg > ltpda

comparison m-toolbox/classes/@pest/tdChi2.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 % tdChi2 computes the chi-square for a parameter estimate.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: tdChi2 computes the chi-square in time domain for an input
5 % pest. The system measured outputs, inputs and models must be
6 % contained in the plist. Also whitening filters for each
7 % output may be taken into account.
8 %
9 % CALL: obj = tdChi2(objs,pl);
10 %
11 % INPUTS: obj - must be a single pest
12 %
13 % <a href="matlab:utils.helper.displayMethodInfo('pest', 'tdChi2')">Parameters Description</a>
14 %
15 % VERSION: $Id: tdChi2.m,v 1.5 2011/04/08 08:56:25 hewitson Exp $
16 %
17 % HISTORY: 08-02-2011 G. Congedo
18 %
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21 function varargout = tdChi2(varargin)
22
23 %%% Check if this is a call for parameters
24 if utils.helper.isinfocall(varargin{:})
25 varargout{1} = getInfo(varargin{3});
26 return
27 end
28
29 import utils.const.*
30 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
31
32 % Collect input variable names
33 in_names = cell(size(varargin));
34 for ii = 1:nargin,in_names{ii} = inputname(ii);end
35
36 % Collect all AOs and plists
37 [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
38 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
39
40
41 % combine plists
42 pl = parse(pl, getDefaultPlist());
43
44 % Extract necessary parameters
45 inputs = pl.find('inputs');
46 outputs = pl.find('outputs');
47 models = pl.find('models');
48 WF = pl.find('WhiteningFilters');
49 Ncut = pl.find('Ncut');
50 Npad = pl.find('Npad');
51
52 if nargout == 0
53 error('### tdChi2 cannot be used as a modifier. Please give an output variable.');
54 end
55
56 if ~all(isa(pests, 'pest'))
57 error('### tdChi2 must be only applied to pest objects.');
58 end
59
60 % Determine the class
61 outClass = class(outputs);
62
63 % Check ouputs
64 if isempty(outputs)
65 error('### Please give the outputs');
66 end
67 switch outClass
68 case 'ao'
69 N = numel(outputs);
70 if N>1
71 error('### Please give the outputs in a MATRIX');
72 end
73 case 'matrix'
74 N = numel(outputs);
75 if N>1
76 error('### Please give the outputs in a COLLECTION of MATRIXs');
77 end
78 if outputs.ncols>1
79 outputs = outputs.';
80 end
81 case 'collection'
82 N = numel(outputs.objs);
83 for ii=1:N
84 if outputs.objs{ii}.ncols>1
85 outputs.objs{ii} = outputs.objs{ii}.';
86 end
87 end
88 otherwise
89 error('### Unknown class for outputs');
90 end
91
92 % Check inputs
93 if isempty(inputs)
94 error('### Please give the inputs');
95 end
96 if ~strcmp(class(inputs),outClass)
97 error('### Please give inputs as the same class of outputs');
98 end
99 switch outClass
100 case 'ao'
101 if numel(inputs)>1
102 error('### Please give the inputs in a MATRIX');
103 end
104 case 'matrix'
105 if numel(inputs)>1
106 error('### Please give the inputs in a COLLECTION of MATRIXs');
107 end
108 if inputs.ncols>1
109 inputs = inputs.';
110 end
111 case 'collection'
112 for ii=1:numel(inputs.objs)
113 if inputs.objs{ii}.ncols>1
114 inputs.objs{ii} = inputs.objs{ii}.';
115 end
116 end
117 end
118
119 % Check models
120 if isempty(models)
121 error('### Please give the transfer function models');
122 end
123 switch outClass
124 case 'ao'
125 if ~strcmp(class(models),'smodel')
126 error('### Please, give the transfer function in a SMODEL.');
127 end
128 if numel(models)>1
129 error('### The size of the transfer function SMODEL does not match with the number of inputs/outputs.');
130 end
131 case {'matrix','collection'}
132 if ~strcmp(class(models),'matrix')
133 error('### Please, give the transfer function in a MATRIX of SMODELs.');
134 end
135 for ii=1:N
136 if strcmp(outClass,'matrix')
137 checkSz = models.nrows~=outputs.nrows || models.ncols~=inputs.nrows;
138 elseif strcmp(outClass,'collection')
139 checkSz = models.nrows~=outputs.objs{ii}.nrows || models.ncols~=inputs.objs{ii}.nrows;
140 end
141 if checkSz
142 error('### The size of the transfer function MATRIX does not match with the number of inputs/outputs.');
143 end
144 end
145 end
146
147 % Check whitening filters
148 whiten = ~isempty(WF);
149 if whiten
150 switch outClass
151 case 'ao'
152 if numel(WF)>1 || ~any(strcmp(class(WF),{'miir','fiir','filterbank'}))
153 error('### Please give the whitening filters in a FIIR, MIIR or FILTERBANK');
154 end
155 case {'matrix','collection'}
156 if ~strcmp(class(WF),'matrix')
157 error('### Please give the whitening filters in a MATRIX of FIIRs, MIIRs or FILTERBANKs');
158 end
159 for ii=1:N
160 if strcmp(outClass,'matrix')
161 checkSz = WF.nrows~=outputs.nrows && WF.ncols~=outputs.nrows;
162 elseif strcmp(outClass,'collection')
163 checkSz = WF.nrows~=outputs.objs{ii}.nrows && WF.ncols~=outputs.objs{ii}.nrows;
164 end
165 if checkSz
166 error('### The size of the whitening filters MATRIX does not match with the number of outputs.');
167 end
168 end
169 end
170 end
171
172
173 % Actual computation
174
175 chi2 = zeros(N,1);
176 Ndata = zeros(N,1);
177
178 % Subs unwanted params
179 if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
180 for kk=1:numel(models.objs)
181 models.objs(kk).setParams(pests.names,pests.y);
182 models.objs(kk).subs(setdiff(models.objs(kk).params,pests.names));
183 end
184 else
185 models.setParams(pests.names,pests.y);
186 models.subs(setdiff(models.params,pests.names));
187 end
188
189 for ii=1:N
190
191 % Time-domain template
192 if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
193 template = fftfilt(inputs.objs{ii},models,plist('Npad',Npad));
194 else
195 template = fftfilt(inputs,models,plist('Npad',Npad));
196 end
197
198 % Residues
199 if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
200 res = template-outputs.objs{ii};
201 else
202 res = template-outputs;
203 end
204
205 % Whiten
206 if whiten
207 res = filter(res,WF);
208 end
209
210 % Split-out transients
211 if ~isempty(Ncut) || Ncut~=0
212 res = split(res,plist('samples',[Ncut+1,Inf]));
213 end
214
215 % Compute chi2 & dof
216 if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
217 chi2(ii) = sum(sum((res.objs.y).^2));
218 Ndata(ii) = max(size(res.objs.y))*min(size(res.objs.y));
219 else
220 chi2(ii) = sum((res.y).^2);
221 Ndata(ii) = numel(res.y);
222 end
223
224 end
225
226 % Compute total chi2 & dof
227 chi2 = sum(chi2);
228 dof = sum(Ndata)-numel(pests.y);
229
230 % Output pest
231 out = copy(pests,1);
232 out = out.setChi2(chi2/dof);
233 out = out.setDof(dof);
234 out = out.setName(['tdChi2(' pests.name ')']);
235
236 out.addHistory(getInfo('None'), pl, pest_invars(:), [pests(:).hist]);
237
238 % Set outputs
239 if nargout > 0
240 varargout{1} = out;
241 end
242
243 end
244
245 %--------------------------------------------------------------------------
246 % Get Info Object
247 %--------------------------------------------------------------------------
248 function ii = getInfo(varargin)
249 if nargin == 1 && strcmpi(varargin{1}, 'None')
250 sets = {};
251 pl = [];
252 else
253 sets = {'Default'};
254 pl = getDefaultPlist;
255 end
256 % Build info object
257 ii = minfo(mfilename, 'pest', 'ltpda', utils.const.categories.helper, '$Id: tdChi2.m,v 1.5 2011/04/08 08:56:25 hewitson Exp $', sets, pl);
258 end
259
260 %--------------------------------------------------------------------------
261 % Get Default Plist
262 %--------------------------------------------------------------------------
263 function plout = getDefaultPlist()
264 persistent pl;
265 if exist('pl', 'var')==0 || isempty(pl)
266 pl = buildplist();
267 end
268 plout = pl;
269 end
270
271 function plo = buildplist()
272 plo = plist();
273
274 % Outputs
275 p = param({'Outputs', 'The system outputs. Must be an AO, a MATRIX or a COLLECTION of MATRIXs, one per each experiment.'}, paramValue.EMPTY_DOUBLE);
276 plo.append(p);
277
278 % Inputs
279 p = param({'Inputs', 'The system inputs. Must be an AO, a MATRIX or a COLLECTION of MATRIXs, one per each experiment.'}, paramValue.EMPTY_DOUBLE);
280 plo.append(p);
281
282 % Models
283 p = param({'Models', 'The system transfer function SMODELs. Must be a SMODEL or a MATRIX of SMODELs.'}, paramValue.EMPTY_DOUBLE);
284 plo.append(p);
285
286 % Whitening filters
287 p = param({'WhiteningFilters', 'The output whitening filters. Must be a MIIR, FIIR, FILTERBANK or a MATRIX.'}, paramValue.EMPTY_DOUBLE);
288 plo.append(p);
289
290 % Ncut
291 p = param({'Ncut', 'The number of points to cut out initial whitening filter transients.'}, paramValue.EMPTY_DOUBLE);
292 plo.append(p);
293
294 % Npad
295 p = param({'Npad', 'The number of points to zero-pad the input for ifft. If left empty, a data length is assumed.'}, paramValue.EMPTY_DOUBLE);
296 plo.append(p);
297
298 end