s2zni.m

很多matlab的源代码

function [nn,dd,g1] = s2zni(n,d,sf,ty,fm)
% S2ZNI Convert H(s) to H(z) using numerical algorithms
%
%	[NZ,DZ,G] = S2ZNI(NS,DS,SF,TY,FM) s2z map by numerical algorithms 
%	NS, DS are num and den of H(s)  NZ, DZ are num and den of H(z) 
%	SF = the sampling frequency (HERTZ), 
%       TY is one of the following numerical algorithms (DEFAULT: TY='trap'):
%       'forw', 'cent'		    backward/forward/central difference
%       'back' or 'rect'	    backward difference/rectangular rule
%       'simp', 'tick'  	    Simpsons/Ticks/ [THESE YIELD UNSTABLE H(z)]
%       'trap' or 'tust' or 'bili'  Trapezoidal rule (Bilinear or Tustin)
%       'ams2' 'ams3' 'ams4' 'ams5 'Adams/Moulton/Schnider's order 2,3,4 and 5
%       FM = OPTIONAL freq (Hz) at which gain of H(z) and H(s) is matched.
%	FM = [FA, FD] specifies analog and digital match frequencies in HERTZ.
%	G is a gain such that the DC gain of G*H(z) matches H(s) (if nonzero).
%
%       S2ZNI (with no input arguments) invokes the following example:
%
%       % Transform H(s)=2/(s+2) using backward difference, with SF=10Hz,  
%	match the gains of H(s) and H(z) at 3Hz and compare/plot results
%       >>[nr,dr] = s2zni(2,[1 2],10,'back',3)
%       >>[h1 p1 f1]=tfplot('s',2,[1 2],[0 10]);
%       >>[h2 p2 f2]=tfplot('z',nr,dr,[0 0.5]);plot(f1/sf,h1,f2,h2)
 

% 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 s2zni,disp('Strike a key to see results of the example')
pause,[nr,dr,gg]=s2zni(2,[1 2],10,'back',3)
[h1 p1 f1]=tfplot('s',2,[1 2],[0 10]);[h2 p2 f2]=tfplot('z',nr,dr,[0 0.5]);
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
plot(f1/10,h1,f2,h2),grid,return,end

if nargin<5,fm=[];end,
if nargin<4,ty='trap';end,
if nargin<3,sf=1;end
t0=1/sf;ty=ty(1:4);
if ty=='tust'|ty=='bili',ty='trap';end
n1=n;d1=d;
c=2/t0;
if ~isempty(fm),fa=fm(1);
if length(fm)>1,fd=fm(2);else,fd=fa;end
if 2*t0*fd>=1,error('Digital match freq. exceeds 0.5fsamp'),return,end
fd1=fd/sf;

if fa>0,
fn=fd1/fa;c=2*pi*fa/tan(pi*fd1);
if fn~=1,
if ty~='trap',[n,d]=lp2af('lp',n,d,fn);t0=1;end
end,end,end


if ty=='trap'
%if ~isempty(fm),c=2/t0;else,c=2*pi*fa/tan(pi*fd1);end
na=[1 -1];da=[1 1]/c; %Trapezoidal
elseif ty=='back'|ty=='rect',na=[1 -1];da=t0*[1 0];
elseif ty=='simp',na=[1 0 -1];da=(t0/3)*[1 4 1]; %Simpsons
elseif ty=='tick',na=[1 0 -1];da=t0*[0.3584 1.2832 0.3584]; %Ticks(2)
elseif ty=='cent',na=[1 0 -1];da=2*t0*[1 0]; %Central(2)
elseif ty=='forw',na=[1 -1];da=t0; %forward(1)
elseif ty=='ams2',na=[1 -1 0];da=(t0/12)*[5 8 -1]; %Ord2
elseif ty=='ams3',na=[1 -1 0 0];da=(t0/24)*[9 19 -5 1]; %Ord3
elseif ty=='ams4'na=[1 -1 0 0 0];da=(t0/720)*[251 646 -264 106 -19]; %Ord4
elseif ty=='ams5'na=[1 -1 0 0 0 0];da=(t0/1440)*[475 1427 -798 482 -173 27]; 
else,error('Unknown mapping type'),return,end
[nn,dd]=polymap(n,d,na,da);


if ty ~='trap'
if ~isempty(fm),
%if length(w0)==1,
j=sqrt(-1);
w2=j*2*pi*fa;
gc=abs(polyval(n1,w2)./polyval(d1,w2));
w1=exp(j*2*pi*fd1);
gd=abs(polyval(nn,w1)./polyval(dd,w1));
g0=gc/gd;nn=g0*nn;
end
end

gc=abs(polyval(n1,0)./polyval(d1,0));
if sum(dd)~=0,
gd=abs(sum(nn)/sum(dd));
if gd~=0,
g1=gc/gd;
end
end