line source
+ − % LP2Z converts a continous TF in to a discrete TF.
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: LP2Z converts a continous (lapalce notation) transfer
+ − % function in a discrete (z transfor formalism) partial fractioned transfer
+ − % function.
+ − %
+ − % Input function can be in rational, poles and zeros or partial fractions
+ − % form:
+ − % The rational case assumes an input function of the type:
+ − %
+ − % a(1)s^m + a(2)s^{m-1} + ... + a(m+1)
+ − % H(s) = --------------------------------------
+ − % b(1)s^n + b(2)s^{n-1} + ... + b(n+1)
+ − %
+ − % The poles and zeros case assumes an impution function of the type:
+ − %
+ − % (s-Z(1))(s-Z(2))...(s-Z(n))
+ − % H(s) = K ---------------------------
+ − % (s-P(1))(s-P(2))...(s-P(n))
+ − %
+ − % The partail fraction case assumes an input function of the type:
+ − %
+ − % R(1) R(2) R(n)
+ − % H(s) = -------- + -------- + ... + -------- + K(s)
+ − % s - P(1) s - P(2) s - P(n)
+ − %
+ − % LP2Z function is also capable to handling input transfer functions with
+ − % multiple poles.
+ − %
+ − %
+ − %
+ − % CALL:
+ − % pfstruct = lp2z(varargin)
+ − % [res,poles,dterm] = lp2z(varargin)
+ − % [pfstruct,res,poles,dterm] = lp2z(varargin)
+ − %
+ − %
+ − %
+ − % INPUTS:
+ − % Input parameters are:
+ − % 'FS' sampling frequency
+ − % 'RRTOL' the repeated-root tolerance, default value is
+ − % 1e-15. If two roots differs less than rrtolerance value, they
+ − % are reported as multiple roots.
+ − % 'INOPT' define the input function type
+ − % 'PF' use this option if you want to input a function
+ − % expanded in partial fractions. Then you have to input the
+ − % vectors containing:
+ − % 'RES' the vector containig the residues
+ − % 'POLES' The vector containing the poles
+ − % 'DTERMS' The vector containing the direct terms
+ − % '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
+ − % 'MODE' Is the mode used for the calculation of the roots of a
+ − % polynomial. 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)
+ − %
+ − %
+ − % OUTPUTS:
+ − % Function output is a structure array with two fields: 'num' and
+ − % 'den' respectively. Each term of the partial fraction expansion
+ − % can be seen as rational transfer function in z domain. These
+ − % reduced transfer functions can be used to filter in parallel a
+ − % series of data. This representation is particularly useful when
+ − % poles with multiplicity higher than one are present. In this
+ − % case the corresponding terms of the partial fraction expansion
+ − % are really rational transfer functions whose order depends from
+ − % the pole multiplicity [1].
+ − % As a sencond option residues, poles and direct term are output
+ − % as standard vectors (this possibility works only with simple
+ − % poles).
+ − % The third possibility is to call the function with four output,
+ − % 1st is a struct, 2nd are residues, 3rd are ples and 4th is the
+ − % direct term.
+ − %
+ − % NOTE1: The function is capable to convert in partial fractions
+ − % only rational functions with numerator of order lower or equal
+ − % to the order of the denominator. When the order of the
+ − % numerator is higher than the order of the denominator (improper
+ − % rational functions) the expansion in partial fractions is not
+ − % useful and other methods are preferable.
+ − %
+ − % NOTE2: 'SYM' calculation mode requires Symbolic Math Toolbox
+ − % to work
+ − %
+ − %
+ − %
+ − % 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: lp2z.m,v 1.11 2009/02/10 12:14:29 luigi Exp $
+ − %
+ − % HISTORY: 16-04-2008 L Ferraioli
+ − % Creation
+ − %
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = lp2z(varargin)
+ −
+ −
+ − %% Finding input parameters
+ −
+ − % Default values for the parameters
+ − fs = 10;
+ − inopt = 'RAT';
+ − tol = 1e-15;
+ − mode = 'DBL';
+ −
+ − % 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},'FS')
+ − fs = 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
+ −
+ −
+ − %% Conversion to partial fractions and discretization
+ −
+ − % swithing between input options
+ − switch inopt
+ − case 'RAT'
+ − % Finding proper parameters
+ − for jj=1:length(varargin)
+ − if strcmp(varargin{jj},'NUM')
+ − num = varargin{jj+1};
+ − end
+ − if strcmp(varargin{jj},'DEN')
+ − den = varargin{jj+1};
+ − end
+ − end
+ − % conversion in partail fractions cP is the vector of residues, cP is the
+ − % vector of poles and cK is the vector of direct terms
+ − [cR, cP, cK, mul] = utils.math.cpf('INOPT', 'RAT', 'NUM', num,...
+ − 'DEN', den, 'MODE', mode, 'RRTOL', tol);
+ −
+ − case 'PZ'
+ − % Finding proper parameters
+ − for jj=1:length(varargin)
+ − if strcmp(varargin{jj},'ZEROS')
+ − zer = varargin{jj+1};
+ − end
+ − if strcmp(varargin{jj},'POLES')
+ − pol = varargin{jj+1};
+ − end
+ − if strcmp(varargin{jj},'GAIN')
+ − gain = varargin{jj+1};
+ − end
+ − end
+ − % conversion in partail fractions cP is the vector of residues, cP is the
+ − % vector of poles and cK is the vector of direct terms
+ − [cR, cP, cK, mul] = utils.math.cpf('INOPT', 'PZ', 'ZEROS', zer,...
+ − 'POLES', pol, 'GAIN', gain, 'RRTOL', tol);
+ −
+ − case 'PF'
+ − % Finding proper parameters
+ − for jj=1:length(varargin)
+ − if strcmp(varargin{jj},'RES')
+ − cR = varargin{jj+1};
+ − end
+ − if strcmp(varargin{jj},'POLES')
+ − cP = varargin{jj+1};
+ − end
+ − if strcmp(varargin{jj},'DTERMS')
+ − cK = varargin{jj+1};
+ − end
+ − end
+ − % finding poles multiplicity
+ − mul = mpoles(cP,tol,0);
+ − end
+ −
+ − % Checking for not proper functions
+ − if length(cK)>1
+ − error('Unable to directly discretize not proper continuous functions ')
+ − end
+ −
+ − % disctretization
+ − T = 1/fs;
+ −
+ − dpls = exp(cP.*T); % poles discretization
+ −
+ − % Construct a struct array in which the elements of the partial
+ − % fractions expansion are stored. Each partial fraction can be
+ − % considered as rational function with proper numerator and denominator
+ − % coefficients
+ − pfstruct = struct();
+ −
+ − indxmul = 0;
+ − for kk=1:length(mul)
+ − if mul(kk)==1
+ − pfstruct(kk).num = cR(kk).*T;
+ − pfstruct(kk).den = [1 -1*dpls(kk)];
+ − else
+ − coeffs = PolyLogCoeffs(1-mul(kk),exp(cP(kk)*T));
+ − coeffs = flipdim(coeffs,2);
+ − pfstruct(kk).num = (cR(kk)*T^mul(kk)/prod(1:mul(kk)-1)).*coeffs;
+ − pfstruct(kk).den = BinomialExpansion(mul(kk),exp(cP(kk)*T));
+ − indxmul = 1;
+ − end
+ − end
+ −
+ − % Inserting the coefficients relative to the presence of a direct term
+ − if ((length(cK)==1)&&(not(cK(1)==0)))
+ − pfstruct(length(mul)+1).num = cK.*T;
+ − pfstruct(length(mul)+1).den = [1 0];
+ − end
+ −
+ − % Output also the vector of residues, poles and direct term if non multiple
+ − % poles are found
+ − if indxmul == 0
+ − res = cR.*T;
+ − poles = dpls;
+ − dterm = cK.*T;
+ − else
+ − res = [];
+ − poles = [];
+ − dterm = [];
+ − end
+ −
+ − % Output empty fields struct in case of no input
+ − if isempty(cR)||isempty(cP)
+ − pfstruct.num = [];
+ − pfstruct.den = [];
+ − end
+ −
+ − %% Output data
+ −
+ − if indxmul == 1
+ − disp( ' Found poles with multiplicity higher than one. Results are output in struct form only. ')
+ − end
+ −
+ − if nargout == 1
+ − varargout{1} = pfstruct;
+ − elseif nargout == 3
+ − varargout{1} = res;
+ − varargout{2} = poles;
+ − varargout{3} = dterm;
+ − elseif nargout == 4
+ − varargout{1} = pfstruct;
+ − varargout{2} = res;
+ − varargout{3} = poles;
+ − varargout{4} = dterm;
+ − else
+ − error('Unespected number of outputs! Output must be 1, 3 or 4')
+ − end
+ −
+ −
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: Eulerian Number
+ − %
+ − % DESCRIPTION: Calculates the Eulerian Number
+ − %
+ − % HISTORY: 12-05-2008 L Ferraioli
+ − % Creation
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − function enum = eulerian(n,k)
+ −
+ − enum = 0;
+ − for jj = 0:k+1;
+ − enum = enum + ((-1)^jj)*(prod(1:n+1)/(prod(1:n+1-jj)*prod(1:jj)))*(k-jj+1)^n;
+ − end
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: PloyLogCoeffs
+ − %
+ − % DESCRIPTION: Computes the coefficients of a poly-logarithm expansion for
+ − % negative n values
+ − %
+ − % HISTORY: 12-05-2008 L Ferraioli
+ − % Creation
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − function coeffs = PolyLogCoeffs(n,x)
+ −
+ − if n>=0
+ − error('n must be a negative integer!')
+ − else
+ − % Now we keep the absolute value of n
+ − n = -1*n;
+ − end
+ −
+ − coeffs = [];
+ − for ii=0:n;
+ − coeffs = [coeffs eulerian(n,ii)*(x^(n-ii))];
+ − end
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: BinomialExpansion
+ − %
+ − % DESCRIPTION: Performs the binomial expamsion of the expression (1-x)^n
+ − %
+ − % HISTORY: 12-05-2008 L Ferraioli
+ − % Creation
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − function binexp = BinomialExpansion(n,x)
+ −
+ − binexp = [];
+ − for jj=0:n
+ − binexp = [binexp ((-1)^jj)*(prod(1:n)/(prod(1:n-jj)*prod(1:jj))*(x^jj))];
+ − end