s2zinvar.m

很多matlab的源代码

function [n,d,g1]=s2zinvar(nn,dd,sf,ty,w0,tol)
% S2ZINVAR Response invariant mapping of H(s) to H(z)
%
%       [N,D,G] = S2ZINVAR(NN,DD,SF,TY,FM,TOL)
%	Response invariant mapping of H(s) to H(z)
%	NN and DD are the numerator and denominator of H(s)
%	SF is the sampling frequency in HERTZ (default: SF=1)
%	TY= 'i' (impulse), 's' (step), or 'r' (ramp) invariance.
%	FM is an optional match frequency (HERTZ) for gain of H(s) and H(z)
%       If FM = 0, the D.C. gain of H(z) and H(s) are matched.
%	FM = [FA,FD] specifies analog and digital match frequencies in HERTZ.
%       TOL specifies the tolerance for repeated roots (default: TOL= 0.002)
%       N and D return numerator and denominator of H(z).
%       G is a gain such that DC gain of G*H(z) matches H(s) (if nonzero).
%
%       S2INVAR (with no input arguments) invokes the following example:
%
%	% Transform H(s)=2/(s+2) using ramp invariance, SF=10Hz and a 
%       % a match frequency of FM = 3Hz and plot results.
%          >>[nr,dr,gr] = s2zinvar(2,[1 2],10,'r',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 s2zinvar,disp('Strike a key to see results of the example')
pause,[nr,dr,gr]=s2zinvar(2,[1 2],10,'r',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<6,tol=0.002;end,
if nargin<5,w0=[];end
if nargin<4,ty='i';end,
if nargin<3,sf=1;end,
ts=1/sf;n1=nn;d1=dd;ch=0;ty1=ty;ty=ty(1);

if ty=='i',if length(ty1)>1,ch=1;end,end

if ~isempty(w0),fa=w0(1);if length(w0)>1,fd=w0(2);else,fd=fa;end,
if 2*ts*fd>=1,error('Digital match freq. exceeds 0.5fsamp'),return,end

fd1=fd/sf;

%if ty=='i',
if fa>0,fn=fd1/fa;
if fn~=1,
[nn,dd]=lp2af('lp',nn,dd,fn);sf=1;ts=1;end
end
end
%end


if ty=='s',dd=[dd 0];end,if ty=='r',dd=[dd 0 0];end
ln=length(nn);ld=length(dd);
while nn(ln)==0 & dd(ld)==0  %Check for cancellation of s terms
nn(ln)=[];dd(ld)=[];
ln=length(nn);
ld=length(dd);
end

[a,p,k]=tf2pf(nn,dd,0,tol);
if ch==1,%if ln>=ld-1,lk=length(k);sp=-0.5*sum(real(a));
if length(a)>=length(p)-1,lk=length(k);sp=-0.5*sum(real(a));
if lk>0,k(lk)=k(lk)+sp;else,if sp~=0,k=sp;end,end,end,end

w=p;[p,i]=sort(-abs(w));p=w(i);a=a(i);lp=length(p);q=ones(lp,1);
for j=2:lp,if abs(p(j)-p(j-1))<tol,q(j)=q(j-1)+1;end,end %Multiple poles.
q=q(lp:-1:1);p=p(lp:-1:1);a=a(lp:-1:1);qm=max(q);

if qm==1,p=exp(p*ts);lp=0;
else
x=zeros(qm);x(1,1)=1;n0=qm-1;c=zeros(n0);c(1,1)=1;
for j=1:n0,for i=1:j
if i==1,c(i,j)=c(1,1);elseif j>i & j<=n0 & i>1,c(i,j)=i*c(i,j-1)+i*c(i-1,j-1);
elseif i==j,c(i,j)=i*c(i-1,j-1);else,end,end,end,
x(2:qm,2:qm)=c';fk=diag(x)';
t=ts.^(0:qm-1);pp=[];aa=[];
while lp>0,m=q(1);
lp=lp-m;r1=exp(p(1)*ts);
r=r1.^(0:m-1);y=x(1:m,1:m);
for j=1:m,y(:,j)=y(:,j)*r(j);end,
am=a(1:m);am=am(m:-1:1);
for j=1:m,y(j,:)=y(j,:)*am(j)*t(j)/fk(j);end,
s=sum(y);aa=[aa;s.'];
pp=[pp;r1*ones(m,1)];a(1:m)=[];p(1:m)=[];q(1:m)=[];
end,
a=aa;p=pp;
end

lk=length(k);
if lk>0,a=[k(lk);a];p=[0;p];k(lk)=[];end  %Modify PFExp form
[n,d]=pf2tf(a,p,k,tol); %Find H(z)/z

while n(1)==0,n(1)=[];end,
while d(1)==0,d(1)=[];end

if ty=='i',n=[n 0];end  %Find H(z)
if ty=='s',
if length(d)>2,d=deconv(d,[1 -1]);
else,n=conv(n,[1 -1]);end
end
if ty=='r',n=n/ts;
if length(d)>3,d=deconv(d,[1 -2 1]);
else,n=conv(n,[1 -2 1]);end
end
n=real(n);d=real(d);ln=length(n);ld=length(d);
while n(ln)==0 & d(ld)==0,
n(ln)=[];d(ld)=[];ln=length(n);ld=length(d);
end
n=[zeros(1,ld-ln) n];d=[zeros(1,ln-ld) d];
if ty=='i',if ts ~=1,n=n*ts;end,end

if ~isempty(w0),
%if fa==fd
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(n,w1)./polyval(d,w1));g0=gc/gd;n=g0*n;
%end,
end
gc=abs(polyval(n1,0)./polyval(d1,0));
if sum(d)~=0,
gd=abs(sum(n)/sum(d));
if gd~=0,g1=gc/gd;end
end