tf2pf.m

很多matlab的源代码

function [r,p,k] = tf2pf(u,v,ty1,tol)
% TF2PF Transfer Function to Partial Fraction Expansion.
%
%       [R,P,K] = TF2PF(NUM,DEN,TY1,TOL) PF expansion of H(x) = N(x)/D(x)
%       NUM, DEN = num and den of transfer function in descending power 
%	IF TY1 = 1, DEN is assumed to contain the DENOMINATOR ROOTS 
%	NOTE: If TY1 = 1, the leading coefficient of D(x) must be 1
%       TOL = tolerance for repeated roots. [default: TOL = 0.002].
%       WARNING: TF2PF can be unstable especially for repeated poles. 
%
%	K contains direct terms of the form [..x^2 x^1 x^0].
%       R contains the residues
%	P contains the poles. (if repeared, order is [(x-p)^1; (x-p)^2; ...])
%
%	USAGE: If H(s) = (s-9)/[(2s+1)(s+3)(s-2)], use TY = 1 and find 
%	[r,p,k]=tf2pf([0.5 -4.5],[-0.5;-3;2],1) % D(s) is normalized
%
%       TF2PF (with no input arguments) invokes the following example:
%
%       % Find the PF expansion of H(s)=[2*s*s+3s+4]/[s*s+4s+3] 
%       >>[res,pol,cons] = tf2pf([2 3 4],[1 4 3])
%       >>[nn,dd] = pf2tf(res,pol,cons)		%reassemble in TF form

% Modified version of RESIDUE to accomodate user supplied roots and tolerence
% Modified by A.Ambardar, 5-03-1991. With permission from the Mathworks, Inc.


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help tf2pf,disp('Strike a key to see results of the example')
pause,[res,pol,cons]=tf2pf([2 3 4],[1 4 3]),disp('Strike a key to continue'),
pause,[nn,dd]=pf2tf(res,pol,cons),return,end

if nargin<4,tol=0.002;end,
if nargin<3,ty1=0;end
%if ty1==1,r=sort(v(:));v=poly(r);end
if ty1==1
   v=v(:);k=find(v==0);lk=length(k);
   if lk, v(k)=[];v=[v;zeros(lk,1)];end
   lu=length(u);lv=length(v);
while u(lu)==0 & v(lv)==0,u(lu)=[];v(lv)=[];lu=lu-1;lv=lv-1;end,
r=v(:);v=poly(r);
else
   lu=length(u);lv=length(v);
   %if ty1==1
      while u(lu)==0 & v(lv)==0,
         u(lu)=[];
         v(lv)=[];
%         r(1)=[]; %% Not needed??? Sep 16,2003
         lu=lu-1;
         lv=lv-1;
      end      
   %end
end
 
   
u=u(:).';v=v(:).';k=[];
f=find(u ~= 0);g=find(v ~= 0);
if length(f),u=u(f(1):length(u));end,if length(g),v=v(g(1):length(v));end
u=u./v(1);v=v./v(1);
%u,v, %%%%%%%

if length(u)>=length(v),[k,u]=deconv(u,v);end
if ty1==0,r=roots(v);r=tol*round(20*r/tol)/20;end


nr=real(r);ni=imag(r);ra=abs(nr-round(nr))<=tol;
nnr=round(nr).*ra+nr.*(1-ra);
ia=abs(ni-round(ni))<=tol;
nni=round(ni).*ia+ni.*(1-ia);
r=nnr+sqrt(-1)*nni;
[p,i]=sort(-abs(r));p=r(i);lp=length(p);
q=zeros(2,lp);q(1,1)=p(1);q(2,1)=1;
j=1;rep=0;rep1=0;rep2=0;
for i=2:lp,
av=q(1,j)./q(2,j);
if abs(av - p(i))<=tol
q(1,j)=q(1,j)+p(i);q(2,j)=q(2,j)+1;rep=1;
else,
j=j+1;q(1,j)=p(i);q(2,j)=1;
end,
end,
if rep==0,
[pp,i]=sort(-imag(r));pp=r(i);
qq=zeros(2,lp);qq(1,1)=pp(1);qq(2,1)=1;j=1;
for i=2:lp,
av=qq(1,j)./qq(2,j);
if abs(av - pp(i))<=tol
qq(1,j)=qq(1,j)+pp(i);qq(2,j)=qq(2,j)+1;rep1=1;
else,
j=j+1;qq(1,j)=pp(i);qq(2,j)=1;
end,
end
if rep1>0,p=pp;q=qq;rep=rep1;
else
[pp,i]=sort(-real(r));pp=r(i);
qq=zeros(2,lp);qq(1,1)=pp(1);qq(2,1)=1;j=1;
for i=2:lp,
av=qq(1,j)./qq(2,j);if abs(av - pp(i))<=tol
qq(1,j)=qq(1,j)+pp(i);qq(2,j)=qq(2,j)+1;rep2=1;
else,j=j+1;qq(1,j)=pp(i);qq(2,j)=1;
end,
end
if rep2>0,q=qq;p=pp;rep=rep2;end
end
end
q(1,1:j)=q(1,1:j)./q(2,1:j);desired=1;if rep & desired,indx = 0;
for i=1:j,for ii=1:q(2,i),indx=indx+1;p(indx)=q(1,i);end,end,end
c=zeros(length(p),1);
if rep,
v=poly(p);next=0;for i=1:j,pole=q(1,i);n=q(2,i);
for indx=1:n,next=next + 1;c(next)=resi2(u,v,pole,n,indx);end,end
else,for i=1:j,temp = poly(p([1:i-1, i+1:j]));
c(i) = polyval(u,p(i)) ./ polyval(temp,p(i));end,end
r=real(c).*(abs(real(c))>100*eps)+sqrt(-1)*imag(c).*(abs(imag(c))>100*eps);
p=p(:);
p=real(p).*(abs(real(p))>100*eps)+sqrt(-1)*imag(p).*(abs(imag(p))>100*eps);