biline.m

xnsy算法的matlab实现

function [z1,p1,k1]=biline(z,p,k,method)
% 双线性变换：[z1,p1,k1]=biline(z,p,k,method)
%   如果method为'z2w'；则从z域变换到w域
%             为'w2z'；则从w域变换到z域
% 一般形式为：1. [z1,p1,k1]=biline(z,p,k,'z2w')
%             2. [num,den]=biline(z,p,k,'w2z')
%             3. [z1,p1,k1]=biline(num,den,'z2w')
%             4. [num1,den1]=biline(num,den,'z2w')




if (nargin==3)
   method=k;
   [z,p,k]=tf2zp(z,p);
   z=z';
   p=p';
   
end;

if (nargin>4|nargin<3)
   disp('Parament input error!');
   break;
end;

if method(1)=='z'  % 'z->w'; z=(1+w)/(1-w)
   z_n=length(z);
   p_n=length(p);

   z_is_1=z_n-nnz(z+1);
   p_is_1=p_n-nnz(p+1);
   
   z=nonzeros(z+1)'-1;
   p=nonzeros(p+1)'-1;

   z1=(z-1)./(z+1);
   p1=(p-1)./(p+1);
   k1=(-1)^(p_n-z_n)*k*prod(z+1)*2^z_is_1/(prod(p+1)*2^p_is_1);

   z1=[z1,ones(1,p_n-z_n)];

elseif method(1)=='w'  % 'w->z'; w=(z-1)/(z+1)
   z_n=length(z);
   p_n=length(p);

   z_is_1=z_n-nnz(z-1);
   p_is_1=p_n-nnz(p-1);

   z=nonzeros(z-1)'+1;
   p=nonzeros(p-1)'+1;

   z1=(1+z)./(1-z);
   p1=(1+p)./(1-p);
   k1=k*prod(1-z)*(-2)^z_is_1/prod(1-p)*(-2)^p_is_1;
 
   z1=[z1,-1*ones(1,p_n-z_n)];
end;


if (nargout==2)
   [z1,p1]=zp2tf(z1',p1',k1);
end;