--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/startpoles.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,179 @@
+% 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
\ No newline at end of file