Mercurial > hg > ltpda
view m-toolbox/classes/+utils/@math/cpf.m @ 33:5e7477b94d94 database-connection-manager
Add known repositories list to LTPDAPreferences
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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