--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/zDomainFit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,719 @@
+% zDomainFit performs a fitting loop to identify model order and
+% parameters.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: zDomainFit fit a partial fraction model to frequency
+%     response data using the function utils.math.vdfit.
+%
+%     The function performs a fitting loop to automatically identify model
+%     order and parameters in z-domain. Output is a z-domain model expanded
+%     in partial fractions:
+%
+%             z*r1            z*rN
+%     f(s) = ------- + ... + -------
+%            z - p1          z - pN
+%
+%     The identification loop stop when the stop condition is reached.
+% 
+%     Output poles and residues are those with minimum Mean Square Error.
+% 
+%     The stop criterion is based on three different approaches, that can be
+%     chosen by setting the value of the CONDTYPE parameter:
+%
+%     1) Mean Squared Error and variation 
+%     (CONDTYPE = 'MSE')
+%     Check if the normalized mean squared error is lower than the value specified in
+%     the parameter FITTOL and if the relative variation of the mean squared error 
+%     is lower than the value specified in the parameter MSEVARTOL.
+%     E.g. FITTOL = 1e-3, MSEVARTOL = 1e-2 will search for a fit with
+%     normalized magnitude error lower than 1e-3 and MSE relative
+%     variation lower than 1e-2.
+%
+%     2) Residuals Log difference and mean squared error variation
+%     (CONDTYPE = 'RLD')
+%     Log Residuals difference
+%     Check if the minimum of the logarithmic difference between data and
+%     residuals is larger than a specified value. 
+%     E.g. if the tolerance value is set to 2 (FITTOL = 2), the function 
+%     ensures that the difference between data and residuals is at lest 2 
+%     orders of magnitude lower than data itselves.
+%     Mean Squared Error Variation
+%     Check if the relative variation of the mean squared error is lower than
+%     MSEVARTOL.
+%
+%     3) Residuals spectral flatness and mean squared error variation
+%     (CONDTYPE = 'RSF')
+%     Residuals Spectral Flatness
+%     In case of a fit on noisy data, the residuals from a good fit are
+%     expected to be as much as possible similar to a white noise. This
+%     property can be used to test the accuracy of a fit procedure. In
+%     particular it can be tested that the spectral flatness coefficient of
+%     the residuals is larger than a certain quantity sf such that 0<sf<1.
+%     E.g. if the tolerance value is set to 0.2 (FITTOL = 0.2), the function 
+%     ensures that the spectral flatness coefficient of the residuals is 
+%     larger than 0.2.
+%     Root Mean Squared Error
+%     Check if the relative variation of the mean squared error is lower than
+%     MSEVARTOL.
+%
+%     Both in the first, second and third approaches the fitting loop ends 
+%     when the two stopping conditions are satisfied.
+%
+%     The function can also perform a single loop without taking care of
+%     the stop conditions. This happens when the 'AUTOSEARCH' parameter is
+%     set to 'off'.
+%
+%     If you provide more than one AO as input, they will be fitted
+%     together with a common set of poles.
+%
+% CALL:         mod = zDomainFit(a, pl)
+%
+% INPUTS:      a  - input AOs to fit to. If you provide more than one AO as
+%                   input, they will be fitted together with a common set
+%                   of poles. Only frequency domain (fsdata) data can be
+%                   fitted. Each non-fsdata object will be ignored. Input
+%                   objects must have the same number of elements.
+%              pl - parameter list (see below)
+%
+% OUTPUTS:
+%               mod - matrix object containing filterbanks of
+%                     parallel miir filters for each input AO.
+%                     Useful fit information are stored in the objects
+%                     procinfo:
+%                     FIT_RESP  - model frequency response.
+%                     FIT_RESIDUALS - analysis object containing the fit
+%                     residuals.
+%                     FIT_MSE - analysis object containing the mean squared
+%                     error progression during the fitting loop.
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'zDomainFit')">Parameters Description</a>
+%
+% Note: all the input objects are assumed to caontain the same X
+% (frequencies) values
+%
+%
+% EXAMPLES:
+%
+% 1) Fit to a frequency-series using the 'MSE' conditioning criterion for
+% fit accuracy
+%
+%   % Create a frequency-series AO
+%   pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000);
+%   a = ao(pl_data);
+%
+%   % Fitting parameter list
+%   pl_fit = plist('FS',[],...
+%   'AutoSearch','on',...
+%   'StartPoles',[],...
+%   'StartPolesOpt','clog',...
+%   'maxiter',5,...
+%   'minorder',2,...
+%   'maxorder',20,...
+%   'weights',[],...
+%   'weightparam','abs',...
+%   'CONDTYPE','MSE',...
+%   'FITTOL',1e-3,... % check if MSE is lower than 1e-3
+%   'MSEVARTOL',1e-2,...
+%   'Plot','off',...
+%   'ForceStability','off',...
+%   'CheckProgress','off');
+%
+%   % Do fit
+%   b = zDomainFit(a, pl_fit);
+%
+% 2) Fit to a frequency-series using the 'RLD' conditioning criterion for
+% fit accuracy
+%
+%   % Create a frequency-series AO
+%   pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000);
+%   a = ao(pl_data);
+%
+%   % Fitting parameter list
+%   pl_fit = plist('FS',[],...
+%   'AutoSearch','on',...
+%   'StartPoles',[],...
+%   'StartPolesOpt','clog',...
+%   'maxiter',5,...
+%   'minorder',2,...
+%   'maxorder',20,...
+%   'weights',[],...
+%   'weightparam','abs',...
+%   'CONDTYPE','RLD',...
+%   'FITTOL',2,... % check if log10(abs(data))-log10(abs(fit_residuals)) > 2
+%   'MSEVARTOL',1e-2,...
+%   'Plot','off',...
+%   'ForceStability','off',...
+%   'CheckProgress','off');
+%
+%   % Do fit
+%   b = zDomainFit(a, pl_fit);
+%
+% 3) Fit to a frequency-series using the 'RSF' conditioning criterion for
+% fit accuracy
+%
+%   % Create a frequency-series AO
+%   pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000);
+%   a = ao(pl_data);
+%
+%   % Fitting parameter list
+%   pl_fit = plist('FS',[],...
+%   'AutoSearch','on',...
+%   'StartPoles',[],...
+%   'StartPolesOpt','clog',...
+%   'maxiter',5,...
+%   'minorder',2,...
+%   'maxorder',20,...
+%   'weights',[],...
+%   'weightparam','abs',...
+%   'CONDTYPE','RSF',...
+%   'FITTOL',0.7,... % check if residuals spectral flatness is larger than 0.7
+%   'MSEVARTOL',1e-2,...
+%   'Plot','off',...
+%   'ForceStability','off',...
+%   'CheckProgress','off');
+%
+%   % Do fit
+%   b = zDomainFit(a, pl_fit);
+%
+% VERSION:     $Id: zDomainFit.m,v 1.39 2011/08/15 09:46:44 hewitson Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = zDomainFit(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
+  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+
+  %%% Decide on a deep copy or a modify
+  bs = copy(as, nargout);
+  inhists = [as.hist];
+
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+
+  if nargout == 0
+    error('### zDomainFit cannot be used as a modifier. Please give an output variable.');
+  end
+
+  %%%%% Extract necessary parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  extpoles = find(pl, 'StartPoles'); % Check if external poles are providied
+  spolesopt = 0;
+  if isempty(extpoles) % if no external poles set them internally
+    splopt = find(pl, 'StartPolesOpt');
+    switch lower(splopt)
+      case 'real'
+        spolesopt = 1;
+      case 'clog'
+        spolesopt = 2;
+      case 'clin'
+        spolesopt = 3;
+      otherwise
+        error('### Unknown value for parameter ''StartPolesOpt''');
+    end
+  end
+
+  maxiter = find(pl, 'maxiter'); % set the maximum number of iterations
+  minorder = find(pl, 'minorder'); % set the minimum function order
+  maxorder = find(pl, 'maxorder');% set the maximum function order
+
+  extweights = find(pl, 'weights'); % check if external weights are provided
+  if isa(extweights, 'ao')
+    extweights = extweights.y;
+  end
+  weightparam = 0;
+  if isempty(extweights) % set internally the weights on the basis of the input options
+    wtparam = find(pl, 'weightparam');
+    switch lower(wtparam)
+      case 'ones'
+        weightparam = 1;
+      case 'abs'
+        weightparam = 2;
+      case 'sqrt'
+        weightparam = 3;
+      otherwise
+        error('### Unknown value for parameter ''weightparam''');
+    end
+  end
+
+  % decide to plot or not
+  plt = find(pl, 'plot');
+  switch lower(plt)
+    case 'on'
+      showplot = 1;
+    case 'off'
+      showplot = 0;
+    otherwise
+      error('### Unknown value for parameter ''plot''');
+  end
+
+  % Make a decision between Fit conditioning type
+  condtype = find(pl, 'CONDTYPE');
+  condtype = upper(condtype);
+  switch condtype
+    case 'MSE'
+      ctp = 'chivar'; % use normalized mean squared error value and relative variation
+      lrscond = find(pl, 'FITTOL');
+      % give an error for strange values of lrscond
+      if lrscond<0
+        error('!!! Negative values for FITTOL are not allowed !!!')
+      end
+      % handling data
+      lrscond = -1*log10(lrscond);
+      % give a warning for strange values of lrscond
+      if lrscond<0
+        warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
+      end
+    case 'RLD'
+      ctp = 'lrsmse'; % use residuals log difference and MSE relative variation
+      lrscond = find(pl, 'FITTOL');
+      % give a warning for strange values of lrscond
+      if lrscond<0
+        error('!!! Negative values for FITTOL are not allowed !!!')
+      end
+      if lrscond<1
+        warning('You are searching for a frequency by frequency residuals log difference of %s', num2str(lrscond))
+      end
+    case 'RSF'
+      ctp = 'rftmse'; % use residuals spectral flatness and MSE relative variation
+      lrscond = find(pl, 'FITTOL');
+      % give a warning for strange values of lrscond
+      if lrscond<0 || lrscond>1
+        error('!!! Values <0 or >1 for FITTOL are not allowed when CONDTYPE is RSF !!!')
+      end
+    otherwise
+      error('### Unknown value for parameter ''CONDTYPE''');
+  end
+
+  % Tolerance for the MSE relative variation
+  msevar = find(pl, 'MSEVARTOL');
+  % handling data
+  msevar = -1*log10(msevar);
+  % give a warning for strange values of msevar
+  if msevar<0
+    warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
+  end
+
+  % decide to stabilize or not the model
+  stab = find(pl, 'ForceStability');
+  switch lower(stab)
+    case 'on'
+      stabfit = 1;
+    case 'off'
+      stabfit = 0;
+    otherwise
+      error('### Unknown value for parameter ''ForceStability''');
+  end
+
+  % decide to disp or not the fitting progress in matlab command window
+  prg = find(pl, 'CheckProgress');
+  switch lower(prg)
+    case 'on'
+      spy = 1;
+    case 'off'
+      spy = 0;
+  end
+
+  % decide to perform or not a full automatic model search
+  autos = find(pl, 'AutoSearch');
+  switch lower(autos)
+    case 'on'
+      fullauto = 1;
+    case 'off'
+      fullauto = 0;
+    otherwise
+      error('### Unknown value for parameter ''AutoSearch''');
+  end
+
+  %%%%% End Extract necessary parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+  %%%%% Fitting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  %%% Fit parameters
+  params = struct('spolesopt',spolesopt,...
+    'extpoles', extpoles,...
+    'Nmaxiter',maxiter,...
+    'minorder',minorder,...
+    'maxorder',maxorder,...
+    'weightparam',weightparam,...
+    'extweights', extweights,...
+    'plot',showplot,...
+    'ctp',ctp,...
+    'lrscond',lrscond,...
+    'msevar',msevar,...
+    'stabfit',stabfit,...
+    'dterm',0,...
+    'spy',spy,...
+    'fullauto',fullauto);
+
+  %%% extracting elements from AOs
+
+  % Finding the index of the first fsdata
+  for gg = 1:numel(bs)
+    if isa(bs(gg).data, 'fsdata')
+      prm = gg;
+      break
+    end
+  end
+
+  fs = find(pl, 'FS');
+  if isempty(fs) && isnan(bs(prm).data.fs)
+    fs = max(bs(prm).data.x)*2;
+  elseif isempty(fs) && ~isnan(bs(prm).data.fs)
+    fs = bs(prm).data.fs;
+  end
+
+  y = zeros(length(bs(prm).data.getY),numel(bs)); % initialize input vector
+  k = numel(bs(prm).data.getY); % setting a comparison constant
+  idx = true(numel(bs),1); % initialize the control index
+  for jj=1:numel(bs)
+    % checking that AOs are fsdata and skipping non fsdata objects
+    if ~isa(bs(jj).data, 'fsdata')
+      % skipping data if non fsdata
+      warning('!!! %s expects ao/fsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
+      idx(jj) = false; % set the corresponding value of the control index to false
+    else
+      % preparing data for fit
+      yt = bs(jj).data.getY;
+      if numel(yt)~=k
+        error('Input AOs must have the same number of elements')
+      end
+      if size(yt,2)>1 % wish to work with columns
+        y(:,jj) = yt.';
+      else
+        y(:,jj) = yt;
+      end
+      if fs < max(bs(jj).data.x)*2
+        warning('!!! %s cannot fit data when max(f) > fs/2 ', mfilename);
+      end
+    end
+  end
+  % reshaping y to contain only Y from fsdata
+  y = y(:,idx);
+
+  %%% extracting frequencies
+  % Note: all the objects are assumed to caontain the same X (frequencies) values
+  f = bs(prm).data.getX;
+
+  %%% Fitting loop
+  [res,poles,dterm,mresp,rdl,mse] = utils.math.autodfit(y,f,fs,params);
+
+  %%%%% End Fitting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+  %%%%% Building output AOs with model responses, model parameters are
+  % added in the procinfo as parallel bank of miir objects %%%%%%%%%%%%%%%
+
+  %   [a,b] = size(mresp);
+  for kk = 1:numel(bs)
+    if idx(kk) % build outputs
+
+      % Constructing a vector of miir objects
+      %       pfilts = [];
+      for hh=1:length(res(:,kk))
+        mod(hh,kk) = miir(res(hh,kk), [ 1 -poles(hh)], fs);
+        %         mod(hh,kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+        mod(hh,kk).setName(sprintf('fit(%s)', ao_invars{kk}));
+        %         pfilts = [pfilts ft];
+      end
+
+      bmod(kk) = filterbank(plist('filters',mod(:,kk),'type','parallel'));
+      
+      bmod(kk).setName(sprintf('fit(%s)', ao_invars{kk}));
+      
+      % bmod(kk).addHistory(getInfo('None'), pl, ao_invars(kk), inhist);
+
+
+
+      % Output also, model response, residuals and mse in the procinfo
+      rsp = mresp(:,kk);
+      bs(kk).data.setY(rsp);
+      bs(kk).setFs(fs);
+      % clear errors
+      bs(kk).clearErrors;
+
+      % Set output AO name
+      bs(kk).setName(sprintf('fit_resp(%s)', ao_invars{kk}));
+      % Add history
+      %         bs(kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+
+      res_ao = copy(bs(kk),1);
+      trdl = rdl(:,kk);
+      res_ao.data.setY(trdl);
+      res_ao.setFs(fs);
+
+      % Set output AO name
+      res_ao.setName(sprintf('fit_residuals(%s)', ao_invars{kk}));
+      % Add history
+      %         res_ao(kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+
+      d = cdata();
+      tmse = mse(:,kk);
+      d.setY(tmse);
+      mse_ao = ao(d);
+
+      % Set output AO name
+      mse_ao.setName(sprintf('fit_mse(%s)', ao_invars{kk}));
+
+      procpl = plist('fit_resp',bs(kk),...
+        'fit_residuals',res_ao,...
+        'fit_mse',mse_ao);
+
+      bmod(kk).setProcinfo(procpl);
+
+    else % in case of non fsdata input
+      bmod(kk) = filterbank(mirr());
+
+    end
+
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  %%%%%% Set outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Set output
+  if nargout == 1
+    if numel(bs)==1
+      bmod.setName(sprintf('fit(%s)', ao_invars{:}));
+      bmod.addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+      varargout{1} = bmod;
+    else
+      mmod = matrix(bmod);
+      mmod.setName(sprintf('fit(%s)', ao_invars{:}));
+      mmod.addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+      varargout{1} = mmod;
+    end
+  else
+    % multiple output is not supported
+    error('### Multiple output is not supported ###')
+  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, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: zDomainFit.m,v 1.39 2011/08/15 09:46:44 hewitson Exp $', sets, pl);
+  ii.setModifier(false);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;  
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function pl = buildplist()
+  pl = plist();
+
+  % AutoSearch
+  p = param({'AutoSearch', ['''on'': Parform a full automatic search for the<br>'...
+    'transfer function order. The fitting<br>'...
+    'procedure will stop when stop conditions<br>'...
+    'defined are satisfied.<br>'...
+    '''off'': Perform a fitting loop as long as the<br>'...
+    'number of iteration reach ''maxiter''. The order<br>'...
+    'of the fitting function will be that<br>'...
+    'specified in ''MINORDER''.']}, ...
+    {1, {'on', 'off'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % StartPoles
+  p = param({'StartPoles', ['A vector of starting poles. Providing a fixed<br>'...
+    'set of starting poles fixes the function<br>'...
+    'order. If it is left empty starting poles are<br>'...
+    'internally assigned.']}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+
+  % StartPolesOpt
+  p = param({'StartPolesOpt', ['Define the characteristics of internally<br>'...
+    'assigned starting poles. Admitted values<br>'...
+    'are:<ul>'...
+    '<li>''real'' linear-spaced real poles</li>'...
+    '<li>''clog'' log-spaced complex poles</li>'...
+    '<li>''clin'' linear-spaced complex poles</li></ul>']}, ...
+    {2, {'real', 'clog', 'clin'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % MaxIter
+  p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(50));
+  pl.append(p);
+
+  % MinOrder
+  p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
+  pl.append(p);
+
+  % MaxOrder
+  p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(20));
+  pl.append(p);
+
+  % Weights
+  p = param({'Weights', ['A vector with the desired weights. If a single<br>'...
+    'Ao is input weights must be a Nx1 vector where<br>'...
+    'N is the number of elements in the input Ao. If<br>'...
+    'M Aos are passed as input, then weights must<br>'...
+    'be a NxM matrix. If it is leaved empty weights<br>'...
+    'are internally assigned basing on the input<br>'...
+    'parameters']}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+
+  % Weightparam
+  p = param({'weightparam', ['Specify the characteristics of the internally<br>'...
+    'assigned weights. Admitted values are:<ul>'...
+    '<li>''ones'' assigns weights equal to 1 to all data.</li>'...
+    '<li>''abs'' weighs data with <tt>1./abs(y)</tt></li>'...
+    '<li>''sqrt'' weighs data with <tt>1./sqrt(abs(y))</tt></li>']}, ...
+    {2, {'ones', 'abs', 'sqrt'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % CONDTYPE
+  p = param({'CONDTYPE', ['Fit conditioning type. Admitted values are:<ul>'...
+    '<li>''MSE'' Mean Squared Error and variation</li>'...
+    '<li>''RLD'' Log residuals difference and mean squared error variation</li>'...
+    '<li>''RSF'' Residuals spectral flatness and mean squared error variation</li></ul>']}, ...
+    {1, {'MSE', 'RLD', 'RSF'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % FITTOL
+  p = param({'FITTOL', 'Fit tolerance.'}, paramValue.DOUBLE_VALUE(1e-3));
+  pl.append(p);
+
+  % MSEVARTOL
+  p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
+    'relative variation of the mean squared error is<br>'...
+    'smaller than the value specified. This<br>'...
+    'option is useful for finding the minimum of the Chi-squared.']}, ...
+    paramValue.DOUBLE_VALUE(1e-2));
+  pl.append(p);
+
+  % Plot
+  p = param({'Plot', 'Plot results of each fitting step.'}, ...
+    {2, {'on', 'off'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % ForceStability
+  p = param({'ForceStability', 'Force poles to be stable'}, ...
+    {2, {'on', 'off'}, paramValue.SINGLE});
+  pl.append(p);
+
+  % CheckProgress
+  p = param({'CheckProgress', 'Display the status of the fit iteration.'}, ...
+    {2, {'on', 'off'}, paramValue.SINGLE});
+  pl.append(p);
+
+  %   pl = plist('FS',[],...
+  %   'AutoSearch','on',...
+  %   'StartPoles',[],...
+  %   'StartPolesOpt','clog',...
+  %   'maxiter',50,...
+  %   'minorder',2,...
+  %   'maxorder',20,...
+  %   'weights',[],...
+  %   'weightparam','abs',...
+  %   'CONDTYPE','MSE',...
+  %   'FITTOL',1e-3,...
+  %   'MSEVARTOL',1e-2,...
+  %   'Plot','off',...
+  %   'ForceStability','off',...
+  %   'CheckProgress','off');
+end
+% END
+
+
+% PARAMETERS:
+%
+%           'FS'          - It is the sampling frequency. If it is left
+%                           empty sampling frequency is searched in the
+%                           input AOs or is is calculated as 2 of the
+%                           maximum frequency reported in AOs xvalues.
+%                           [Default []].
+%           'AutoSearch'  - 'on': Parform a full automatic search for the
+%                           transfer function order. The fitting
+%                           procedure will stop when stop conditions are
+%                           satisfied. [Default]
+%                           'off': Perform a fitting loop as long as the
+%                           number of iteration reach 'maxiter'. The order
+%                           of the fitting function will be that
+%                           specified in 'minorder'.
+%           'StartPoles'  - A vector of starting poles. Providing a fixed
+%                           set of starting poles fixes the function
+%                           order. If it is left empty starting poles are
+%                           internally assigned. [Default []]
+%         'StartPolesOpt' - Define the characteristics of internally
+%                           assigned starting poles. Admitted values
+%                           are:
+%                            - 'real' linspaced real poles
+%                            - 'c1' complex poles on unit circle. First
+%                              method [Default]. See help of
+%                              utils.math.startpoles for additional info
+%                            - 'c2' complex poles on unit circle. second
+%                              method. See help of utils.math.startpoles
+%                              for additional info.
+%             'maxiter'   - Maximum number of allowed iteration. [Deafult
+%                           50].
+%                           [default: -inf for each parameter];
+%             'minorder'  - Minimum model function order. [Default 2]
+%             'maxorder'  - Maximum model function order. [Default 20]
+%             'weights'   - A vector with the desired weights. If a single
+%                           Ao is input weights must be a Nx1 vector where
+%                           N is the number of elements in the input Ao. If
+%                           M Aos are passed as input, then weights must
+%                           be a NxM matrix. If it is leaved empty weights
+%                           are internally assigned basing on the input
+%                           parameters. [Default []]
+%           'weightparam' - Specify the characteristics of the internally
+%                           assigned weights. Admitted values are:
+%                           'ones' assigns weights equal to 1 to all data.
+%                           'abs' weights data with 1./abs(y) [Default]
+%                           'sqrt' weights data with 1./sqrt(abs(y))
+%             'CONDTYPE'  - Fit conditioning type. Admitted values are:
+%                             - 'MSE' Mean Squared Error and variation
+%                             [Default]
+%                             - 'RLD' Log residuals difference and mean
+%                             squared error variation
+%                             - 'RSF' Residuals spectral flatness and mean
+%                             squared error variation
+%               'FITTOL'  - Fit tolerance [Default, 1e-3]
+%           'MSEVARTOL'   - This allow to check if the relative variation
+%                           of mean squared error is lower than the value
+%                           sepcified. [Default 1e-2]
+%           'Plot'        - Plot fit result: 'on' or 'off' [default]
+%       'ForceStability'  - Force poles to be stable, values are
+%                           'on' or 'off'. [Default 'off']
+%         'CheckProgress' - Disply the status of the fit iteration.
+%                           Values are 'on and 'off'. [Default 'off']
+%