--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/cpf.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,399 @@
+% CPF finds the partial fraction expansion of the ratio of two polynomials A(s)/B(s).
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: CPF finds the residues, poles and direct terms of the
+% partial fraction expansion of the ratio of two polynomials A(s)/B(s).
+% This function assumes that the input continous filter is written in the
+% rational form or in poles, zeros and gain factorization:
+%
+%         A(s)    a(1)s^m + a(2)s^{m-1} + ... + a(m+1)
+%   H(s)= ---- = --------------------------------------
+%         B(s)    b(1)s^n + b(2)s^{n-1} + ... + b(n+1)
+%
+% or
+%
+%         A(s)     (s-z1)...(s-zn)
+%   H(s)= ---- = g ----------------
+%         B(s)     (s-p1)...(s-pn)
+%
+%
+%   It inputs a plist containing the coefficients vectors and the
+%   repeated-root tolerance.
+%   Eg:
+%       A = [a(1), a(2), ..., a(m+1)]
+%       B = [b(1), b(2), ..., b(m+1)]
+%
+% or
+%
+%       Z = [z(1), z(2), ..., z(m)]
+%       P = [p(1), p(2), ..., p(m)]
+%       G = g (Gain is a scalar)
+%
+%
+%   If there are no multiple roots,
+%
+%      A(s)       R(1)       R(2)             R(n)
+%      ----  =  -------- + -------- + ... + -------- + K(s)
+%      B(s)     s - P(1)   s - P(2)         s - P(n)
+%
+%   The number of poles is n = length(B)-1 = length(R) = length(P).
+%   The direct term coefficient vector is empty if length(A) < length(B),
+%   otherwise length(K) = length(A)-length(B)+1.
+%   K(s) is returned in the form:
+%
+%       K(s) = k(1)*s^(m-n) + ... + k(m-n)*s + k(m-n+1)
+%
+%   so that the output vector of direct terms is:
+%       K = [k(1), ..., k(m-n), k(m-n+1)]
+%
+%   If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the
+%   expansion includes terms of the form
+%                R(j)        R(j+1)                R(j+m-1)
+%              -------- + ------------   + ... + ------------
+%              s - P(j)   (s - P(j))^2           (s - P(j))^m
+%
+%   The function is also capable to convert a partial fraction expanded
+%   function to its rational form by setting the 'PARFRACT' input option.
+%   In this case the output is composed by a plist containing the vactors
+%   of numerator and denominator polynomial coefficients.
+%
+%
+%
+% CALL:     varargout = cpf(varargin)
+%
+%
+%
+% INPUTS:
+%           Input options are:
+%           'INOPT' define the input function type
+%             'RAT' input the continous function in rational form.
+%             then you have to input the vector of coefficients:
+%               'NUM' is the vector with numerator coefficients.
+%               'DEN' is the vector of denominator coefficienets.
+%             'PZ' input the continuous function in poles and
+%             zeros form. Then you have to input the vectors with poles
+%             and zeros:
+%               'POLES' the vector with poles
+%               'ZEROS' the vector with zeros
+%               'GAIN' the value of the gain
+%             'PF' input the coefficients of a partial fraction
+%             expansion of the transfer function. When this option is
+%             setted the function performs the conversion from partial
+%             fraction to rational transfer function. You have to input the
+%             vectors containing the residues, poles and direct terms:
+%               'RES' the vector with residues
+%               'POLES' the vector with poles
+%               'DTERMS' the vector with direct terms
+%           'MODE' Is the used mode for the calculation of the roots of a
+%           polynomial. It is an useful option only with rational functions
+%           at the input. Admitted values are:
+%             'SYM' uses symbolic roots calculation (you need symbolic
+%             math toolbox to use this option)
+%             'DBL' uses the standard numerical matlab style roots
+%             claculation (double precision)
+%           'RRTOL' the repeated-root tolerance default value is
+%           1e-15. If two roots differs less than rrtolerance value, they
+%           are reported as multiple roots
+%
+%
+%
+% OUTPUTS:
+%
+%           When 'INOPT' is set to 'RAT' or 'PZ', outputs
+%           are:
+%           RES vector of residues coefficients
+%           POLES vector of poles coefficients
+%           DTERMS vector of direct terms coefficients
+%           PMul vector of poles multiplicity
+%
+%           When 'INOPT' is setted  to 'PF', outputs are:
+%           NUM the vector with the numerator polynomial
+%           coefficients
+%           DEN the vector with the denominator polynomial
+%           coefficents
+%
+%
+%
+% NOTE:
+%         - 'SYM' option for 'MODE' requires the Symblic Math Toolbox. It
+%         is used only for rational function input
+%
+%
+%
+% EXAMPLES:
+%           - Input a function in rational form and output the partial
+%           fraction expansion
+%           [Res, Poles, DTerms, PMul] = cpf('INOPT', 'RAT',
+%           'NUM', [], 'DEN', [], 'MODE','SYM',  'RRTOL', 1e-15)
+%           - Input a function in poles and zeros and output the partial
+%           fraction expansion
+%           [Res, Poles, DTerms, PMul] = cpf('INOPT', 'PZ',
+%           'POLES', [], 'ZEROS', [], 'GAIN', #,  'RRTOL', 1e-15)
+%           - Input a function in partial fractions and output the rational
+%           expression
+%           [Num, Den] = cpf('INOPT', 'PF', 'POLES', [], 'RES',
+%           [], 'DTERMS', [],  'RRTOL', 1e-15)
+%
+%
+%
+% REFERENCES:
+%         [1] Alan V. Oppenheim, Allan S. Willsky and Ian T. Young, Signals
+%         and Systems, Prentice-Hall Signal Processing Series, Prentice
+%         Hall (June 1982), ISBN-10: 0138097313. Pages 767 - 776.
+%
+%
+%
+% VERSION:     $Id: cpf.m,v 1.4 2008/11/19 11:37:51 luigi Exp $
+%
+%
+% HISTORY:     16-04-2008 L Ferraioli
+%                 Creation
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = cpf(varargin)
+  
+  %% VERSION
+  
+  VERSION  = '$Id: cpf.m,v 1.4 2008/11/19 11:37:51 luigi Exp $';
+  
+  %% Extracting parameters
+  
+  % default parameters
+  inopt = 'RAT';
+  mode = 'DBL';
+  tol = 1e-15;
+  
+  % Finding input parameters
+  if ~isempty(varargin)
+    for j=1:length(varargin)
+      if strcmp(varargin{j},'INOPT')
+        inopt = varargin{j+1};
+      end
+      if strcmp(varargin{j},'MODE')
+        mode = varargin{j+1};
+      end
+      if strcmp(varargin{j},'RRTOL')
+        tol = varargin{j+1};
+      end
+    end
+  end
+  
+  % Switching between input options and setup inputs for proper calculation
+  switch inopt;
+    case 'RAT'
+      % etracting numerator and denominator vectors
+      for jj=1:length(varargin)
+        if strcmp(varargin{jj},'NUM')
+          u = varargin{jj+1};
+        end
+        if strcmp(varargin{jj},'DEN')
+          v = varargin{jj+1};
+        end
+      end
+      
+      % For the conversion we need the denominator factored in poles and the
+      % numerator in polynomial form
+      switch mode
+        case 'DBL'
+          % adopt numerical calculation
+          poles_vect = roots(v);
+        case 'SYM'
+          % adopt symbolic calculation
+%           syms s
+%           % Construct the symbolic polynomial
+%           numel = length(v);
+%           PowerVector = [];
+%           for ii=1:numel
+%             PowerVector = [PowerVector v(ii)*s^(numel-ii)];
+%           end
+%           PowerMatrix = diag(PowerVector);
+%           Polyv = trace(PowerMatrix);
+%           % Solve the polynomial in order to find the roots
+%           Sp = solve(Polyv,s);
+%           % output of the poles vector in Matlab double format
+%           numpoles = length(Sp);
+%           poles_vect = zeros(1,numpoles);
+%           for jj=1:numpoles
+%             poles_vect(jj) = double(Sp(jj));
+%           end
+
+        n = length(v);
+        cN = -1.*v(2:end)./v(1);
+        A = sym(diag(ones(1,n-2),-1));
+        A(1,:) = cN;
+        sol = eig(A);
+        poles_vect = double(sol);
+      end
+      % setting the output option to residues and poles
+      outopt = 'RP';
+      
+    case 'PZ'
+      % Extracting zeros, poles and gain from inputs
+      for jj=1:length(varargin)
+        if strcmp(varargin{jj},'ZEROS')
+          zeros_vect = varargin{jj+1};
+        end
+        if strcmp(varargin{jj},'POLES')
+          poles_vect = varargin{jj+1};
+        end
+        if strcmp(varargin{jj},'GAIN')
+          gain = varargin{jj+1};
+        end
+      end
+      
+      u = poly(zeros_vect).*gain;
+      v = poly(poles_vect);
+      if ((~isempty(v)) && (v(1)~=1))
+        u = u ./ v(1); v = v ./ v(1);   % Normalize.
+      end
+      % setting the output option to residues and poles
+      outopt = 'RP';
+      
+    case 'PF'
+      % Calculate numerator and denominator of a transfer function expanded
+      % in partial fractions
+      % etracting residues, poles and direct terms
+      for jj=1:length(varargin)
+        if strcmp(varargin{jj},'RES')
+          u = varargin{jj+1};
+        end
+        if strcmp(varargin{jj},'POLES')
+          v = varargin{jj+1};
+        end
+        if strcmp(varargin{jj},'DTERMS')
+          k = varargin{jj+1};
+        end
+      end
+      % setting the output option to Transfer Function
+      outopt = 'TF';
+  end
+  
+  
+  
+  %% Partial fractions expansion
+  
+  % Switching between output cases
+  % Note: rational input and poles zeros are equivalent from this point on,
+  % so PF expansion is calculated in the same way
+  switch outopt
+    
+    case 'RP'
+      % Direct terms calculation
+      if length(u) >= length(v)
+        [dterms,new_u]=deconv(u,v);
+      else
+        dterms = 0;
+        new_u = u;
+      end
+      
+      % identification of multiple poles
+      poles_vect = sort(poles_vect); % sort the poles in ascending order
+      mul = mpoles(poles_vect,tol,0); % find the multiplicity
+      
+      mmul = mul;
+      for kk=1:length(mmul)
+        if mmul(kk)>1
+          for hh=1:mmul(kk)
+            mmul(kk-hh+1)=mmul(kk);
+          end
+        end
+      end
+      
+      % finding the residues
+      resids = zeros(length(poles_vect),1);
+      
+      for ii=1:length(poles_vect)
+        
+        den = v;
+        p = [1 -poles_vect(ii)];
+        for hh=1:mmul(ii)
+          den = deconv(den,p);
+        end
+        
+        dnum = new_u;
+        dden = den;
+        
+        c = 1;
+        if mmul(ii)>mul(ii)
+          c = prod(1:(mmul(ii)-mul(ii)));
+          
+          for jj=1:(mmul(ii)-mul(ii))
+            [dnum,dden] = polyder(dnum,dden);
+          end
+          
+        end
+        
+        resids(ii)=(polyval(dnum,poles_vect(ii))./polyval(dden,poles_vect(ii)))./c;
+      end
+      
+      % Converting from partial fractions to rational function
+    case 'TF'
+      % This code is directly taken from matlab 'residue' function
+      [mults,i]=mpoles(v,tol,0);
+      p=v(i); r=u(i);
+      n = length(p);
+      q = [p(:).' ; mults(:).'];   % Poles and multiplicities.
+      v = poly(p); u = zeros(1,n,class(u));
+      for indx = 1:n
+        ptemp = q(1,:);
+        i = indx;
+        for j = 1:q(2,indx), ptemp(i) = nan; i = i-1; end
+        ptemp = ptemp(find(~isnan(ptemp))); temp = poly(ptemp);
+        j = length(temp);
+        if j < n, temp = [zeros(1,n-j) temp]; end
+        u = u + (r(indx) .* temp);
+      end
+      if ~isempty(k)
+        if any(k ~= 0)
+          u = [zeros(1,length(k)) u];
+          k = k(:).';
+          temp = conv(k,v);
+          u = u + temp;
+        end
+      end
+      num = u; den = v;    % Rename.
+  end
+  
+  %% Output data
+  
+  switch outopt
+    case 'RP'
+      if nargout == 1
+        varargout{1} = [resids poles_vect dterms mul];
+      elseif nargout == 2
+        varargout{1} = resids;
+        varargout{2} = poles_vect;
+      elseif nargout == 3
+        varargout{1} = resids;
+        varargout{2} = poles_vect;
+        varargout{3} = dterms;
+      elseif nargout == 4
+        varargout{1} = resids;
+        varargout{2} = poles_vect;
+        varargout{3} = dterms;
+        varargout{4} = mul;
+      else
+        error('Unespected number of outputs! Set 1, 2, 3 or 4')
+      end
+      %     plout = plist('RESIDUES', resids, 'POLES', poles_vect, 'PMul', mul, 'DIRECT_TERMS', dterms);
+      % %     plout = combine(plout, pl);
+    case 'TF'
+      if nargout == 1
+        varargout{1} = [num den];
+      elseif nargout == 2
+        varargout{1} = num;
+        varargout{2} = den;
+      else
+        error('Unespected number of outputs! Set 1 or 2')
+      end
+      %     plout = plist('NUMERATOR', num, 'DENOMINATOR', den);
+  end
+  
+end
+% END
+
+
+
+