+
−% AUTODFIT perform a fitting loop to identify model order and parameters.
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% DESCRIPTION:
+
−%
+
−%     Perform a fitting loop to automatically identify model order and
+
−%     parameters in z-domain. Model identification is performed by 'dtfit'
+
−%     function. Output is a z-domain model expanded in partial fractions:
+
−%
+
−%             z*r1            z*rN
+
−%     f(s) = ------- + ... + ------- + d
+
−%            z - p1          z - pN
+
−%
+
−%     The function attempt to perform first the identification of a model
+
−%     with d = 0, then if the operation do not succeed, it try the
+
−%     identification with d different from zero.
+
−%     Identification loop stop when the stop condition is reached. Six
+
−%     stop criteria are available:
+
−% 
+
−%     Mean Squared Error
+
−%     Check if the normalized mean squared error is lower than the value
+
−%     specified in lrscond:
+
−%     mse < 10^(-1*lrsvarcond)
+
−%
+
−%     Mean Squared Error and variation
+
−%     Check if the normalized mean squared error is lower than the value specified in
+
−%     lrscond and that the relative variation of the mean squared error is lower
+
−%     than the value provided.
+
−%     Checking algorithm is:
+
−%     mse < 10^(-1*lrsvarcond)
+
−%     Dmse < 10^(-1*msevar)
+
−%
+
−%     Log Residuals difference
+
−%     Check if the minimum of the logarithmic difference between data and
+
−%     residuals is larger than a specified value. ie. if the conditioning
+
−%     value is 2, the function ensures that the difference between data and
+
−%     residuals is at lest 2 order of magnitude lower than data itsleves.
+
−%     Checking algorithm is:
+
−%     lsr = log10(abs(y))-log10(abs(rdl));
+
−%     min(lsr) > lrscond;
+
−%
+
−%     Log Residuals difference and Mean Squared Error Variation
+
−%     Check if the log difference between data and residuals is in
+
−%     larger than the value indicated in lsrcond and that the relative variation of
+
−%     the mean squared error is lower than 10^(-1*msevar).
+
−%     Checking algorithm is:
+
−%     lsr = log10(abs(y))-log10(abs(rdl));
+
−%     (lsr > lrscond) && (mse < 10^(-1*lrsvarcond));
+
−% 
+
−%     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 qiantity sf such that 0<sf<1.
+
−% 
+
−%     Residuals Spectral Flatness and root mean squared error variation
+
−%     Check if the spectral flatness coefficient of the residuals is
+
−%     larger than a certain qiantity sf such that 0<sf<1 and that the
+
−%     relative variation of the mean squared error is lower than
+
−%     10^(-1*msevar).
+
−% 
+
−%     Once the loop iteration stops the parameters giving best Mean Square
+
−%     Error are output.
+
−%
+
−% CALL:
+
−%
+
−%     [res,poles,dterm,mresp,rdl,mse] = autodfit(y,f,fs,params)
+
−%
+
−% INPUT:
+
−%
+
−%     - y are the data to be fitted. They represent the frequency response
+
−%     of a certain process.
+
−%     - f is the frequency vector in Hz
+
−%     - fs is the sampling frequency in Hz
+
−%     - params is a struct containing identification parameters
+
−%
+
−%       params.spolesopt = 0 --> use external starting poles
+
−%       params.spolesopt = 1 --> use real starting poles
+
−%       params.spolesopt = 2 --> generates complex conjugates poles of the
+
−%       type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1).
+
−%       params.spolesopt = 3 --> generates complex conjugates poles of the
+
−%       type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2).
+
−%       Default option.
+
−% 
+
−%       params.extpoles = [] --> a vector with the starting poles.
+
−%       Providing a fixed set of starting poles fixes the function order so
+
−%       params.minorder and params.maxorder will be internally set to the
+
−%       poles vecto length.
+
−%
+
−%       params.fullauto = 0 --> Perform a fitting loop as far as the number
+
−%       of iteration reach Nmaxiter. The order of the fitting function will
+
−%       be that specified in params.minorder. If params.dterm is setted to
+
−%       1 the function will fit only with direct term.
+
−%       params.fullauto = 1 --> Parform a full automatic search for the
+
−%       transfer function order. The fitting procedure will stop when the
+
−%       stopping condition defined in params.ctp is satisfied. Default
+
−%       value.
+
−%
+
−%       params.Nmaxiter = # --> Number of maximum iteration per model order
+
−%       parformed. Default is 50.
+
−%
+
−%       params.minorder = # --> Minimum model trial order. Default is 2.
+
−%
+
−%       params.maxorder = # --> Maximum model trial order. Default is 25.
+
−%
+
−%       params.weightparam = 0 --> use external weights
+
−%       params.weightparam = 1 --> fit with equal weights (one) for each
+
−%       data point.
+
−%       params.weightparam = 2 --> weight fit with the inverse of absolute
+
−%       value of data. Default value.
+
−%       params.weightparam = 3 --> weight fit with the square root of the
+
−%       inverse of absolute value of data.
+
−%       params.weightparam = 4 --> weight fit with inverse of the square
+
−%       mean spread
+
−% 
+
−%       params.extweights = [] --> A vector of externally provided weights.
+
−%       It has to be of the same size of input data.
+
−%
+
−%       params.plot = 0 --> no plot during fit iteration
+
−%       params.plot = 1 --> plot results at each fitting steps. default
+
−%       value.
+
−%
+
−%       params.ctp = 'chival' --> check if the value of the Mean Squared
+
−%       Error is lower than 10^(-1*lsrcond).
+
−%       params.ctp = 'chivar' --> check if the value of the Mean Squared
+
−%       Error is lower than 10^(-1*lsrcond) and if the relative variation of mean
+
−%       squared error is lower than 10^(-1*msevar).
+
−%       params.ctp = 'lrs' --> check if the log difference between data and
+
−%       residuals is point by point larger than the value indicated in
+
−%       lsrcond. This mean that residuals are lsrcond order of magnitudes
+
−%       lower than data.
+
−%       params.ctp = 'lrsmse' --> check if the log difference between data
+
−%       and residuals is larger than the value indicated in lsrcond and if
+
−%       the relative variation of root mean squared error is lower than
+
−%       10^(-1*msevar).
+
−%       params.ctp = 'rft' --> check that the residuals spectral flatness
+
−%       coefficient is lerger than the value provided in lsrcond. In this
+
−%       case lsrcond must be such that 0<lsrcond<1.
+
−%       params.ctp = 'rftmse' --> check that the residuals spectral flatness
+
−%       coefficient is lerger than the value provided in lsrcond and if
+
−%       the relative variation of root mean squared error is lower than
+
−%       10^(-1*msevar). In this case lsrcond must be such that
+
−%       0<lsrcond<1.
+
−%
+
−%       params.lrscond = # --> set conditioning value for mean squared
+
−%       error (params.ctp = 'chival' and params.ctp = 'chivar') or set
+
−%       conditioning value for point to point log residuals difference
+
−%       (params.ctp = 'lsr' and params.ctp = 'lrsmse') or set conditioning
+
−%       value for residuals spectral flateness (params.ctp = 'rft' and
+
−%       params.ctp = 'rftmse'). In the last case params.lrscond must be set
+
−%       to 0<lrscond<1. Default is 2. See help for stopfit.m for further
+
−%       remarks.
+
−%
+
−%       params.msevar = # --> set conditioning value for mean squared
+
−%       error variation. This allow to check that the variation of
+
−%       mean squared error is lower than 10^(-1*msevar).Default is 7. See
+
−%       help for stopfit.m for further remarks.
+
−%
+
−%       params.stabfit = 0 --> Fit without forcing poles stability. Default
+
−%       value.
+
−%       params.stabfit = 1 --> Fit forcing poles stability
+
−%
+
−%       params.dterm = 0 --> Try to fit without direct term
+
−%       params.dterm = 1 --> Try to fit with and without direct term
+
−%
+
−%       params.spy = 0 --> Do not display the iteration progression
+
−%       params.spy = 1 --> Display the iteration progression
+
−%
+
−% OUTPUT:
+
−%
+
−%     - res is the vector with model residues r
+
−%     - poles is the vector with model poles p
+
−%     - dterm is the model direct term d
+
−%     - mresp is the model frequency response calculated at the input
+
−%     frequencies
+
−%     - rdl are the residuals between data and model, at the input
+
−%     frequencies
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% VERSION: $Id: autodfit.m,v 1.25 2010/01/27 17:56:11 luigi Exp $
+
−%
+
−% HISTORY:     08-10-2008 L Ferraioli
+
−%                 Creation
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−
+
−function [res,poles,dterm,mresp,rdl,mse] = autodfit(y,f,fs,params)
+
−
+
−  utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
+
−
+
−  % Default input struct
+
−  defaultparams = struct('spolesopt',1, 'Nmaxiter',30, 'minorder',2,...
+
−    'maxorder',25, 'weightparam',1, 'plot',1,...
+
−    'ctp','chival','lrscond',2,'msevar',2,...
+
−    'stabfit',0,'dterm',0,'spy',0,'fullauto',1,'extweights', [],...
+
−    'extpoles', []);
+
−
+
−  names = {'spolesopt','Nmaxiter','minorder',...
+
−    'maxorder','weightparam','plot',...
+
−    'ctp','lrscond','msevar',...
+
−    'stabfit','dterm','spy','fullauto','extweights','extpoles'};
+
−
+
−  % collecting input and default params
+
−  if ~isempty(params)
+
−    for jj=1:length(names)
+
−      if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
+
−        defaultparams.(names{1,jj}) = params.(names{1,jj});
+
−      end
+
−    end
+
−  end
+
−
+
−  % collecting input params
+
−  spolesopt = defaultparams.spolesopt;
+
−  Nmaxiter = defaultparams.Nmaxiter;
+
−  minorder = defaultparams.minorder;
+
−  maxorder = defaultparams.maxorder;
+
−  weightparam = defaultparams.weightparam;
+
−  check = defaultparams.plot;
+
−  stabfit = defaultparams.stabfit;
+
−  ctp = defaultparams.ctp;
+
−  lrscond = defaultparams.lrscond;
+
−  msevar = defaultparams.msevar;
+
−  idt = defaultparams.dterm;
+
−  spy = defaultparams.spy;
+
−  autosearch = defaultparams.fullauto;
+
−  extweights = defaultparams.extweights;
+
−  extpoles = defaultparams.extpoles;
+
−
+
−  if check == 1
+
−    fitin.plot = 1;
+
−    fitin.ploth = figure; % opening new figure window
+
−  else
+
−    fitin.plot = 0;
+
−  end
+
−
+
−  if stabfit % fit with stable poles only
+
−    fitin.stable = 1;
+
−  else % fit without restrictions
+
−    fitin.stable = 0;
+
−  end
+
−
+
−  % Colum vector are preferred
+
−  [a,b] = size(y);
+
−  if a < b % shifting to column
+
−    y = y.';
+
−  end
+
−  [Nx,Ny] = size(y);
+
−
+
−  [a,b] = size(f);
+
−  if a < b % shifting to column
+
−    f = f.';
+
−  end
+
−  
+
−  % in case of externally provided poles
+
−  if ~isempty(extpoles)
+
−    spolesopt = 0;
+
−  end
+
−  if spolesopt == 0 % in case of external poles
+
−    % Colum vector are preferred
+
−    [a,b] = size(extpoles);
+
−    if a < b % shifting to column
+
−      extpoles = extpoles.';
+
−    end
+
−    [Npls,b] = size(extpoles);
+
−    minorder = Npls;
+
−    maxorder = Npls;
+
−  end
+
−  
+
−  if weightparam == 0 % in case of external weights
+
−    % Colum vector are preferred
+
−    [a,b] = size(extweights);
+
−    if a < b % shifting to column
+
−      extweights = extweights.';
+
−    end
+
−  end
+
−
+
−  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−  % Importing package
+
−  import utils.math.*
+
−
+
−  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−  % Fitting
+
−
+
−  % decide to fit with or without direct term according to the input
+
−  % options
+
−  if autosearch
+
−    if idt % full auto identification
+
−      dterm_off = 1;
+
−      dterm_on = 1;
+
−    else % auto ident without dterm
+
−      dterm_off = 1;
+
−      dterm_on = 0;
+
−    end
+
−  else
+
−    if idt % fit only with dterm
+
−      dterm_off = 0;
+
−      dterm_on = 1;
+
−    else % fit without dterm
+
−      dterm_off = 1;
+
−      dterm_on = 0;
+
−    end
+
−  end
+
−
+
−  ext = zeros(Ny,1);
+
−  
+
−  % starting banana mse
+
−  cmse = inf;
+
−  bmse = inf;
+
−
+
−  if dterm_off
+
−    utils.helper.msg(utils.const.msg.PROC1, ' Try fitting without direct term ')
+
−    fitin.dterm = 0;
+
−    fitin.fs = fs;
+
−
+
−    % Weighting coefficients
+
−    if weightparam == 0
+
−      % using external weigths
+
−      utils.helper.msg(utils.const.msg.PROC1, ' Using external weights... ')
+
−      weight = extweights;
+
−    else
+
−      weight = utils.math.wfun(y,weightparam);
+
−    end
+
−
+
−    if autosearch
+
−      order_vect = minorder:maxorder;
+
−    else
+
−      order_vect = minorder:minorder;
+
−    end
+
−    
+
−
+
−    for N = order_vect
+
−
+
−      if spy
+
−        utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
+
−      end
+
−
+
−      % Starting poles
+
−      if spolesopt == 0 % in case of external poles
+
−        utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
+
−        spoles = extpoles;
+
−      else % internally calculated starting poles
+
−        pparams = struct('spolesopt',spolesopt, 'type','DISC', 'pamp', 0.98);
+
−        spoles = utils.math.startpoles(N,f,pparams);
+
−      end
+
−
+
−      % Fitting
+
−      M = 3*N;
+
−      if M > Nmaxiter
+
−        M = Nmaxiter;
+
−      elseif not(autosearch)
+
−        M = Nmaxiter;
+
−      end
+
−      
+
−      clear mlr
+
−      
+
−
+
−      for hh = 1:M
+
−        [res,spoles,dterm,mresp,rdl,mse] = utils.math.vdfit(y,f,spoles,weight,fitin); % Fitting
+
−        
+
−        % decide to store the best result based on mse
+
−        %fprintf('iteration = %d, order = %d \n',hh,N)
+
−
+
−        if norm(mse)<cmse
+
−          %fprintf('nice job \n')
+
−          bres = res;
+
−          bpoles = spoles;
+
−          bdterm = dterm;
+
−          bmresp = mresp;
+
−          brdl = rdl;
+
−          bmse = mse;
+
−          cmse = norm(mse);
+
−        end
+
−
+
−        if spy
+
−          utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
+
−        end
+
−
+
−        %         ext = zeros(Ny,1);
+
−        if autosearch
+
−          for kk = 1:Ny
+
−            % Stop condition checking
+
−            mlr(hh,kk) = mse(:,kk);
+
−            % decide between stop conditioning
+
−            if strcmpi(ctp,'lrs')
+
−              yd = y(:,kk); % input data
+
−            elseif strcmpi(ctp,'lrsmse')
+
−              yd = y(:,kk); % input data
+
−            elseif strcmpi(ctp,'rft')
+
−              yd = mresp(:,kk); % model response
+
−            elseif strcmpi(ctp,'rftmse')
+
−              yd = mresp(:,kk); % model response
+
−            elseif strcmpi(ctp,'chival')
+
−              yd = y(:,kk); % model response              
+
−            elseif strcmpi(ctp,'chivar')
+
−              yd = y(:,kk); % model response    
+
−            else
+
−              error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
+
−            end
+
−            [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
+
−            ext(kk,1) = next;
+
−          end
+
−        else
+
−          for kk = 1:Ny
+
−            % storing mse progression
+
−            mlr(hh,kk) = mse(:,kk);
+
−          end
+
−        end
+
−
+
−        if all(ext)
+
−          utils.helper.msg(utils.const.msg.PROC1, msg)
+
−          break
+
−        end
+
−
+
−      end
+
−      if all(ext)
+
−        break
+
−      end
+
−
+
−    end
+
−  end
+
−
+
−  if dterm_on
+
−    if ~all(ext) % fit with direct term only if the fit without does not give acceptable results (in full auto mode)
+
−      utils.helper.msg(utils.const.msg.PROC1, ' Try fitting with direct term ')
+
−      fitin.dterm = 1;
+
−
+
−      if autosearch
+
−        order_vect = minorder:maxorder;
+
−      else
+
−        order_vect = minorder:minorder;
+
−      end
+
−
+
−      for N = order_vect
+
−
+
−        if spy
+
−          utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
+
−        end
+
−
+
−        % Starting poles
+
−        if spolesopt == 0 % in case of external poles
+
−          utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
+
−          spoles = extpoles;
+
−        else % internally calculated starting poles
+
−          pparams = struct('spolesopt',spolesopt, 'type','DISC', 'pamp', 0.98);
+
−          spoles = utils.math.startpoles(N,f,pparams);
+
−        end
+
−
+
−        % Fitting
+
−        M = 3*N;
+
−        if M > Nmaxiter
+
−          M = Nmaxiter;
+
−        elseif not(autosearch)
+
−          M = Nmaxiter;
+
−        end
+
−        
+
−        clear mlr
+
−        
+
−        for hh = 1:M
+
−          [res,spoles,dterm,mresp,rdl,mse] = utils.math.vdfit(y,f,spoles,weight,fitin); % Fitting
+
−          
+
−          % decide to store the best result based on mse
+
−          if norm(mse)<cmse
+
−            bres = res;
+
−            bpoles = spoles;
+
−            bdterm = dterm;
+
−            bmresp = mresp;
+
−            brdl = rdl;
+
−            bmse = mse;
+
−            cmse = norm(mse);
+
−          end
+
−
+
−          if spy
+
−            utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
+
−          end
+
−
+
−          ext = zeros(Ny,1);
+
−          if autosearch
+
−            for kk = 1:Ny
+
−              % Stop condition checking
+
−              mlr(hh,kk) = mse(:,kk);
+
−              % decide between stop conditioning
+
−              if strcmpi(ctp,'lrs')
+
−                yd = y(:,kk); % input data
+
−              elseif strcmpi(ctp,'lrsmse')
+
−                yd = y(:,kk); % input data
+
−              elseif strcmpi(ctp,'rft')
+
−                yd = mresp(:,kk); % model response
+
−              elseif strcmpi(ctp,'rftmse')
+
−                yd = mresp(:,kk); % model response
+
−              elseif strcmpi(ctp,'chival')
+
−                yd = y(:,kk); % model response
+
−              elseif strcmpi(ctp,'chivar')
+
−                yd = y(:,kk); % model response
+
−              else
+
−                error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
+
−              end
+
−              [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
+
−              ext(kk,1) = next;
+
−            end
+
−          else
+
−            for kk = 1:Ny
+
−              % storing mse progression
+
−              mlr(hh,kk) = mse(:,kk);
+
−            end
+
−          end
+
−
+
−          if all(ext)
+
−            utils.helper.msg(utils.const.msg.PROC1, msg)
+
−            break
+
−          end
+
−
+
−        end
+
−        if all(ext)
+
−          break
+
−        end
+
−
+
−      end
+
−
+
−    end
+
−  end
+
−
+
−  poles = bpoles;
+
−  clear mse
+
−  mse = mlr(:,:);
+
−  
+
−  res = bres;
+
−  dterm = bdterm;
+
−  mresp = bmresp;
+
−  rdl = brdl;
+
−  mse = bmse;
+
−
+
−  if all(ext) == 0
+
−    utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed without reaching the prescribed accuracy. Try changing Nmaxiter or maxorder or accuracy requirements ')
+
−  else
+
−    utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed successfully ')
+
−  end