polymap.m

很多matlab的源代码

function [nn,dd] = polymap(n,d,nz,dz,ty)
% POLYMAP Mapping of polynomials.
%
%         [By,Ay] = POLYMAP(NUMx,DENx,NUMy,DENy,ty) maps H(x)=NUMx(x)/DENx(x) 
%         to H(y)=B(y)/A(y) using the polynomial transformation 
%         x = f(y) = NUMy(y)/DENy(y)
%
%         NUMx,DENx, NUMy,DENy must be in DESCENDING order of x and y.
%         ty=0 for DFT method and ty=1 for convolution method (Def. ty=1)
%         [B,A] are returned as numerator and denominator arrays.
%
%         POLYMAP (with no input arguments) invokes the following example:
%
%         % Transform H(x)=(x+2)/(x*x+4x+3) to H(s) using x=(s*s+4)/s, 
%         >>[ny,dy] = polymap([1 2],[1 4 3],[1 0 4],[0 1 0])


% 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


%REF: Higher order s2z mapping functions: Schneider et al 
%     IEEE Proceedings, vol 79, Nov. 1991

if nargin==0,help polymap,disp('Strike a key to see results of the example')
pause,[ny,dy]=polymap([1 2],[1 4 3],[1 0 4],[0 1 0]),return,end

if nargin<5,ty=1;end
if nargin<4,dz=1;end
while n(1)==0,n(1)=[];end,
while d(1)==0,d(1)=[];end
l1=length(n);
l2=length(d);

%% Do checks
if l1<=1 & l2<=1,nn=n;dd=d;return,end

if ~any(nz) & ~any(dz),error('Invalid mapping'),return,end
if ~any(nz),nn=polyval(n,0);dd=polyval(d,0);return,end
if ~any(dz),
if l1<l2,nn=0;dd=1;
elseif l1>l2, nn=inf,dd=1;
else,nn=n(1);dd=d(1);
end
return,end

while nz(1)==0,nz(1)=[];end
while dz(1)==0,dz(1)=[];end

ln=length(nz);ld=length(dz);

t=max(ln,ld)-1;
if l1>l2,d=[zeros(1,l1-l2) d];elseif l2>l1,n=[zeros(1,l2-l1) n];end

if ty==0
j=sqrt(-1);
m=length(n)-1;l=t*m;k=0:l;z=exp(j*2*k*pi/(l+1));
lz=length(z);na=zeros(1,lz);da=na;
for i=1:ln;na=na.*z+nz(i);end,for i=1:ld;da=da.*z+dz(i);end
x=0;y=0;
for i=0:m,a=(na.^(m-i)).*(da.^i);y=y+n(i+1)*a;x=x+d(i+1)*a;end
nn=real(ifft(y));
nn=[nn(2:l+1) nn(1)];
dd=real(ifft(x));
dd=[dd(2:l+1) dd(1)];
return
end



no = length(n)-1;do = length(d)-1; nn=0; dd=0;

for k = 0:no
 nterm = 1;
   for p = 0:k-1
       nterm = conv(nterm,dz);
   end
 dterm = 1;
   for p = 0:no-k-1
       dterm = conv(dterm,nz);
   end
 b1=n(k+1)*conv(nterm,dterm);l1=length(nn);l2=length(b1);
 if l1>l2,b1=[zeros(1,l1-l2) b1];else nn=[zeros(1,l2-l1) nn];end
 nn=nn+b1;
end


for k = 0:do
 nterm = 1;
   for p = 0:k-1
       nterm = conv(nterm,dz);
   end
 dterm = 1;
   for p = 0:do-k-1
       dterm = conv(dterm,nz);
   end
 b1=d(k+1)*conv(nterm,dterm);l1=length(dd);l2=length(b1);
 if l1>l2,b1=[zeros(1,l1-l2) b1];else dd=[zeros(1,l2-l1) dd];end
 dd=dd+b1;
end