+
−% STARTPOLES defines starting poles for fitting procedures ctfit, dtfit.
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% DESCRIPTION
+
−%
+
−%     Defines the starting poles for the fitting procedure with ctfit or
+
−%     dtfit. Starting poles definition process is different for s-domain
+
−%     and z-domain.
+
−%     s-domain identification:
+
−%     Starting poles can be real or complex in conjugate couples. Real
+
−%     poles are chosen on the [-2*pi*f(1),-2*pi*f(end)] intervall. Complex poles can
+
−%     be defined with the real and imaginary parts logspaced or linespaced
+
−%     on the intervall [2\pi f(1),2\pi f(end)]. Ratio between real and
+
−%     imaginary part can be setted by the user.
+
−%     z-domain identification:
+
−%     Starting poles can be real or come in complex conjugate couples. Real
+
−%     poles are chosen on the [-1,1] intervall. Complex poles are
+
−%     chosen inside the unit circle as:
+
−%     \alfa e^{j\pi\theta}
+
−%     where \theta is linespaced inside the intervall [0,2\pi]. In this
+
−%     case two different methods are used: the first method define angles
+
−%     as \theta = linspace(0,pi,N/2+1), then take out the first element and
+
−%     construct the complex conjugate couples. If N is odd the first
+
−%     element is added as real pole. With this method the last two
+
−%     conjugates poles have a real part very similar to that of the real
+
−%     pole. This may generate problems so a second method is implemented in
+
−%     which the angle are generated as \theta = linspace(0,pi,N/2+2). Then
+
−%     the first and the last elements of the set are taken out and the
+
−%     first element is used only for N odd. The last element instead is
+
−%     never used. This allow to have well distributed poles on the unit
+
−%     circle. The amplitude parameter \alfa can be set by the user.
+
−%
+
−% CALL:
+
−%
+
−%     spoles = startpoles(order,f,params)
+
−%
+
−% INPUT:
+
−%
+
−%     order: is the function order, ie. the number of desired poles
+
−%     f: is the frequency vector in Hz
+
−%     params: is a struct with the setting parameters
+
−%
+
−%       params.type = 'CONT' --> Output a set of poles for s-domain
+
−%       identification
+
−%       params.type = 'DISC' --> Output a set of poles for z-domain
+
−%       identification
+
−%
+
−%       params.spolesopt = 1 --> generate linear spaced real poles on
+
−%       the intervall [-2*pi*f(1),-2*pi*f(end)] for s-domain and [-1,1] for z-domain.
+
−%       params.spolesopt = 2 --> in case of s-domain generates logspaced
+
−%       complex conjugates poles. In case of z-domain generates complex
+
−%       conjugates poles with \theta = linspace(0,pi,N/2+1).
+
−%       params.spolesopt = 3 --> in case of s-domain generates linespaced
+
−%       complex conjugates poles.  In case of z-domain generates complex
+
−%       conjugates poles with \theta = linspace(0,pi,N/2+2). We advice to
+
−%       make use of this option for z-domain identification.
+
−%
+
−%       params.pamp = # --> s-domain: set the ratio \alfa/\beta between
+
−%       poles real and imaginary parts. Adviced value is 0.01.
+
−%       params.pamp = # --> z-domain: set the amplitude of the poles.
+
−%       Adviced value is 0.98.
+
−%
+
−% OUTPUT:
+
−%
+
−%     spoles: is the set of starting poles
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% VERSION: $Id: startpoles.m,v 1.6 2010/04/29 09:00:00 luigi Exp $
+
−% HISTORY:     08-10-2008 L Ferraioli
+
−%                 Creation
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−
+
−function spoles = startpoles(order,f,params)
+
−  
+
−  % Default input struct
+
−  defaultparams = struct('spolesopt',1, 'type','CONT', 'pamp', 0.98);
+
−  
+
−  names = {'spolesopt', 'type', 'pamp'};
+
−  
+
−  % collecting input and default params
+
−  if ~isempty(params)
+
−    for jj=1:length(names)
+
−      if isfield(params, names(jj))
+
−        defaultparams.(names{1,jj}) = params.(names{1,jj});
+
−      end
+
−    end
+
−  end
+
−  
+
−  type = defaultparams.type;
+
−  spolesopt = defaultparams.spolesopt;
+
−  pamp = defaultparams.pamp;
+
−  
+
−  N = order;
+
−  
+
−  % switching between continuous and discrete
+
−  switch type
+
−    case 'CONT'
+
−      switch spolesopt
+
−        
+
−        case 0
+
−          disp(' Using external starting poles')
+
−        
+
−        case 1 % real starting poles
+
−          spoles = -1.*2.*pi.*linspace(f(1),f(end),N).';
+
−          
+
−        case 2 % complex logspaced starting poles
+
−          if f(1)==0
+
−            bet=2.*pi.*logspace(log10(f(2)),log10(f(end)),N/2);
+
−          else
+
−            bet=2.*pi.*logspace(log10(f(1)),log10(f(end)),N/2);
+
−          end
+
−          spoles=[];
+
−          for n=1:length(bet)
+
−            alf=-bet(n)*pamp;
+
−            spoles=[spoles;(alf-1i*bet(n));(alf+1i*bet(n))];
+
−          end
+
−          if (N-2*floor(N/2)) ~= 0
+
−%             rpl = rand(1,1);
+
−%             if rpl > 0
+
−%               rpl = -1*rpl;
+
−%             end
+
−            rpl = -1;
+
−            spoles = [rpl; spoles];
+
−          end
+
−          
+
−        case 3 % complex linspaced starting poles
+
−          bet=linspace(2*pi*f(1),2*pi*f(end),N/2);
+
−          spoles=[];
+
−          for n=1:length(bet)
+
−            alf=-bet(n)*pamp;
+
−            spoles=[spoles;(alf-1i*bet(n));(alf+1i*bet(n))];
+
−          end
+
−          if (N-2*floor(N/2)) ~= 0
+
−%             rpl = rand(1,1);
+
−%             if rpl > 0
+
−%               rpl = -1*rpl;
+
−%             end
+
−            rpl = -1;
+
−            spoles = [rpl; spoles];
+
−          end
+
−          
+
−      end
+
−    case 'DISC'
+
−      switch spolesopt
+
−        
+
−        case 0
+
−          disp(' Using external starting poles')
+
−          
+
−        case 1 % real starting poles
+
−          spoles = linspace(-0.99,0.99,N).';
+
−          
+
−        case 2 % complex starting poles
+
−          ang = linspace(0,pi,N/2+1);
+
−          spoles=[];
+
−          for nn=2:length(ang)
+
−            spoles=[spoles; exp(1i*ang(nn)); exp(-1i*ang(nn))]; % Taking complex conjugate pairs on the unit circle
+
−          end
+
−          if (N-2*floor(N/2)) ~= 0
+
−            rpl = exp(1i*ang(1));
+
−            spoles = [rpl; spoles];
+
−          end
+
−          spoles = spoles.*pamp; % shifting starting ?poles a little inside the unit circle
+
−          
+
−        case 3 % complex starting poles
+
−          ang = linspace(0,pi,N/2+2);
+
−          spoles=[];
+
−          for nn=2:length(ang)-1
+
−            spoles=[spoles; exp(1i*ang(nn)); exp(-1i*ang(nn))]; % Taking complex conjugate pairs on the unit circle
+
−          end
+
−          if (N-2*floor(N/2)) ~= 0
+
−            rpl = exp(1i*ang(1));
+
−            spoles = [rpl; spoles];
+
−          end
+
−          spoles = spoles.*pamp; % shifting starting ?poles a little inside the unit circle
+
−          
+
−      end
+
−      
+
−  end
+
−  
+
−end