pade.m

数字通信第四版原书的例程

function [a,b,c,d]=pade(T,n)
% PADE  Returns an nth order Pade approximation to a time delay, T:
%	                    2           3
%	num    2 -Ts + (-Ts) /2! + (-Ts) /3! + ....
%	--- =  --------------------------- 
%	                   2          3            
%	den    2 +Ts + (Ts) /2! + (Ts) /3!+ .....
%
%	[NUM,DEN]=PADE(T,n)  returns the nth order pade approximation in
%	transfer function form where NUM and DEN contain the polynomial 
%	coefficients in descending powers of s.
%	
%	[A,B,C,D]=PADE(T,n)  returns the nth order pade approximation as
%	a continuous-time state-space system.
%
%	When invoked without lefthand arguments,
%		PADE(T,n)
%	plots the step response of an nth order Pade approximation on
%	the screen.
%
%	See also: C2DT.

%	Andrew C.W. Grace 8-13-89
%	Copyright (c) 1986-93 by the MathWorks, Inc.

% Next command sets the type of pade approximation to be used.
% Set ptype=0 for a smoother step response and low pass 
% filter characteristics, otherwise set ptype=1. 
ptype = 1;

% Polynomial series approximation of exp(-sT) (constant first)
series=1;
for i=1:2*n;
        series(i+1)=-series(i)*T/i;
end

% DENOMINATOR COEFFICIENTS
% First form a matrix to index appropriate elements of 'series'
C=zeros(n,n+1);
for i=1:n
	C(i,:)=i+ptype:n+i+ptype;
end
C(:)=series(C(:));

% Find b, a null vector of C,  so that C*b'=0 
[Q,R,E]=qr(C');
if min(size(R))>1 R=sum(R'); end
[dum,ind]=min(abs(R));
% Denominator solution
b=Q(:,ind).';


% NUMERATOR COEFFICIENTS
% Indexing matrix
C=zeros(n+1,n+1);
for i=1:n+1
	C(i,:)=-i+2:n+2-i;
end
% Address elements of series only for +ve elements of C
C(:)=series(C(:)+n*(C(:)<1))'.*(C(:)>0);
% Numerator solution
a=b*C.'; 

if ~ptype, a(1)=[]; end

% Graphical Output if no left hand arguments.
% Step Response and Bode Plots
if nargout==0
	clg
	hold off
	subplot(211)
	t1=[0:T/40:2*T];
	y1=step(a,b,t1);
	plot(t1,y1)
	hold off
	str=['st';'nd';'rd';'th']; str=str((n<4).*n+4*(n>3),:);
	title(['Step Response of ',num2str(n),'''',str,' order Pade Approx'])
	xlabel('Time (secs)')
	ylabel('Amplitude')
	[mag,phase,w]=bode(a,b);
	subplot(223)
        semilogx(w,20*log10(mag))
        xlabel('Frequency (rad/sec)')
        ylabel('Gain dB')
        subplot(224)
        semilogx(w,phase)
        xlabel('Frequency (rad/sec)')
        ylabel('Phase deg')
	subplot(212)
	title('Bode Plot')
end
if nargout==4
	[a,b,c,d]=tf2ss(a,b);
elseif nargout==3
	[z,p,k]=tf2zp(a,b);
end