--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@pest/tdChi2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,298 @@
+% tdChi2 computes the chi-square for a parameter estimate.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: tdChi2 computes the chi-square in time domain for an input
+%              pest. The system measured outputs, inputs and models must be 
+%              contained in the plist. Also whitening filters for each
+%              output may be taken into account.
+%
+% CALL:        obj = tdChi2(objs,pl);
+%              
+% INPUTS:      obj - must be a single pest
+%
+% <a href="matlab:utils.helper.displayMethodInfo('pest', 'tdChi2')">Parameters Description</a>
+%
+% VERSION:     $Id: tdChi2.m,v 1.5 2011/04/08 08:56:25 hewitson Exp $
+%
+% HISTORY:     08-02-2011 G. Congedo
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = tdChi2(varargin)
+
+  %%% Check if this is a call for parameters
+  if utils.helper.isinfocall(varargin{:})
+    varargout{1} = getInfo(varargin{3});
+    return
+  end
+  
+  import utils.const.*
+  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all AOs and plists
+  [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ 
+  
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+  
+  % Extract necessary parameters
+  inputs    = pl.find('inputs');
+  outputs   = pl.find('outputs');
+  models    = pl.find('models');
+  WF        = pl.find('WhiteningFilters');
+  Ncut      = pl.find('Ncut');
+  Npad      = pl.find('Npad');
+    
+  if nargout == 0
+    error('### tdChi2 cannot be used as a modifier. Please give an output variable.');
+  end
+  
+  if ~all(isa(pests, 'pest'))
+    error('### tdChi2 must be only applied to pest objects.');
+  end
+  
+  % Determine the class
+  outClass = class(outputs);
+  
+  % Check ouputs
+  if isempty(outputs)
+    error('### Please give the outputs');
+  end
+  switch outClass
+    case 'ao'
+      N = numel(outputs);
+      if N>1
+        error('### Please give the outputs in a MATRIX');
+      end
+    case 'matrix'
+      N = numel(outputs);
+      if N>1
+        error('### Please give the outputs in a COLLECTION of MATRIXs');
+      end
+      if outputs.ncols>1
+        outputs = outputs.';
+      end
+    case 'collection'
+      N = numel(outputs.objs);
+      for ii=1:N
+        if outputs.objs{ii}.ncols>1
+          outputs.objs{ii} = outputs.objs{ii}.';
+        end
+      end
+    otherwise
+      error('### Unknown class for outputs');
+  end
+    
+  % Check inputs
+  if isempty(inputs)
+    error('### Please give the inputs');
+  end
+  if ~strcmp(class(inputs),outClass)
+    error('### Please give inputs as the same class of outputs');
+  end
+  switch outClass
+    case 'ao'
+      if numel(inputs)>1
+        error('### Please give the inputs in a MATRIX');
+      end
+    case 'matrix'
+      if numel(inputs)>1
+        error('### Please give the inputs in a COLLECTION of MATRIXs');
+      end
+      if inputs.ncols>1
+        inputs = inputs.';
+      end
+    case 'collection'
+      for ii=1:numel(inputs.objs)
+        if inputs.objs{ii}.ncols>1
+          inputs.objs{ii} = inputs.objs{ii}.';
+        end
+      end
+  end
+  
+  % Check models
+  if isempty(models)
+    error('### Please give the transfer function models');
+  end
+  switch outClass
+    case 'ao'
+      if ~strcmp(class(models),'smodel')
+        error('### Please, give the transfer function in a SMODEL.');
+      end
+      if numel(models)>1
+        error('### The size of the transfer function SMODEL does not match with the number of inputs/outputs.');
+      end
+    case {'matrix','collection'}
+      if ~strcmp(class(models),'matrix')
+        error('### Please, give the transfer function in a MATRIX of SMODELs.');
+      end
+      for ii=1:N
+        if strcmp(outClass,'matrix')
+          checkSz = models.nrows~=outputs.nrows || models.ncols~=inputs.nrows;
+        elseif strcmp(outClass,'collection')
+          checkSz = models.nrows~=outputs.objs{ii}.nrows || models.ncols~=inputs.objs{ii}.nrows;
+        end
+        if checkSz
+          error('### The size of the transfer function MATRIX does not match with the number of inputs/outputs.');
+        end
+      end
+  end
+   
+  % Check whitening filters
+  whiten = ~isempty(WF);
+  if whiten
+    switch outClass
+    case 'ao'
+      if numel(WF)>1 || ~any(strcmp(class(WF),{'miir','fiir','filterbank'}))
+        error('### Please give the whitening filters in a FIIR, MIIR or FILTERBANK');
+      end
+    case {'matrix','collection'}
+      if ~strcmp(class(WF),'matrix')
+        error('### Please give the whitening filters in a MATRIX of FIIRs, MIIRs or FILTERBANKs');
+      end      
+      for ii=1:N
+        if strcmp(outClass,'matrix')
+          checkSz = WF.nrows~=outputs.nrows && WF.ncols~=outputs.nrows;
+        elseif strcmp(outClass,'collection')
+          checkSz = WF.nrows~=outputs.objs{ii}.nrows && WF.ncols~=outputs.objs{ii}.nrows;
+        end
+        if checkSz
+          error('### The size of the whitening filters MATRIX does not match with the number of outputs.');
+        end
+      end
+    end
+  end   
+    
+  
+  % Actual computation
+  
+  chi2 = zeros(N,1);
+  Ndata = zeros(N,1);
+  
+  % Subs unwanted params
+  if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
+    for kk=1:numel(models.objs)
+      models.objs(kk).setParams(pests.names,pests.y);
+      models.objs(kk).subs(setdiff(models.objs(kk).params,pests.names));
+    end
+  else
+    models.setParams(pests.names,pests.y);
+    models.subs(setdiff(models.params,pests.names));
+  end
+  
+  for ii=1:N
+        
+        % Time-domain template
+        if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
+          template = fftfilt(inputs.objs{ii},models,plist('Npad',Npad));
+        else
+          template = fftfilt(inputs,models,plist('Npad',Npad));
+        end
+        
+        % Residues
+        if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
+          res = template-outputs.objs{ii};
+        else
+          res = template-outputs;
+        end
+
+        % Whiten
+        if whiten
+          res = filter(res,WF);
+        end
+
+        % Split-out transients
+        if ~isempty(Ncut) || Ncut~=0
+          res = split(res,plist('samples',[Ncut+1,Inf]));
+        end
+
+        % Compute chi2 & dof
+        if strcmp(outClass,'matrix') || strcmp(outClass,'collection')
+          chi2(ii) = sum(sum((res.objs.y).^2));
+          Ndata(ii) = max(size(res.objs.y))*min(size(res.objs.y));
+        else
+          chi2(ii) = sum((res.y).^2);
+          Ndata(ii) = numel(res.y);
+        end
+    
+  end 
+   
+  % Compute total chi2 & dof
+  chi2 = sum(chi2);
+  dof = sum(Ndata)-numel(pests.y);
+  
+  % Output pest
+  out = copy(pests,1);
+  out = out.setChi2(chi2/dof);
+  out = out.setDof(dof);
+  out = out.setName(['tdChi2(' pests.name ')']);
+  
+  out.addHistory(getInfo('None'), pl, pest_invars(:), [pests(:).hist]);
+     
+  % Set outputs
+  if nargout > 0
+    varargout{1} = out;
+  end
+  
+  end
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pl   = [];
+  else
+    sets = {'Default'};
+    pl   = getDefaultPlist;
+  end
+  % Build info object
+  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);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;  
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function plo = buildplist()
+  plo = plist();
+  
+  % Outputs
+  p = param({'Outputs', 'The system outputs. Must be an AO, a MATRIX or a COLLECTION of MATRIXs, one per each experiment.'}, paramValue.EMPTY_DOUBLE);
+  plo.append(p);
+  
+  % Inputs
+  p = param({'Inputs', 'The system inputs. Must be an AO, a MATRIX or a COLLECTION of MATRIXs, one per each experiment.'}, paramValue.EMPTY_DOUBLE);
+  plo.append(p);
+   
+  % Models
+  p = param({'Models', 'The system transfer function SMODELs. Must be a SMODEL or a MATRIX of SMODELs.'}, paramValue.EMPTY_DOUBLE);
+  plo.append(p);
+  
+  % Whitening filters
+  p = param({'WhiteningFilters', 'The output whitening filters. Must be a MIIR, FIIR, FILTERBANK or a MATRIX.'}, paramValue.EMPTY_DOUBLE);
+  plo.append(p);
+  
+  % Ncut
+  p = param({'Ncut', 'The number of points to cut out initial whitening filter transients.'}, paramValue.EMPTY_DOUBLE);
+  plo.append(p);
+  
+  % Npad
+  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);
+  plo.append(p);
+  
+end