Import.
line source
+ − % 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
+ −
+ −
+ −
+ −