home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / PADE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.5 KB  |  102 lines
  1. function [a,b,c,d]=pade(T,n)
  2. % PADE  Returns an nth order Pade approximation to a time delay, T:
  3. %                        2           3
  4. %    num    2 -Ts + (-Ts) /2! + (-Ts) /3! + ....
  5. %    --- =  --------------------------- 
  6. %                       2          3            
  7. %    den    2 +Ts + (Ts) /2! + (Ts) /3!+ .....
  8. %
  9. %    [NUM,DEN]=PADE(T,n)  returns the nth order pade approximation in
  10. %    transfer function form where NUM and DEN contain the polynomial 
  11. %    coefficients in descending powers of s.
  12. %    
  13. %    [A,B,C,D]=PADE(T,n)  returns the nth order pade approximation as
  14. %    a continuous-time state-space system.
  15. %
  16. %    When invoked without lefthand arguments,
  17. %        PADE(T,n)
  18. %    plots the step response of an nth order Pade approximation on
  19. %    the screen.
  20. %
  21. %    See also: C2DT.
  22.  
  23. %    Andrew C.W. Grace 8-13-89
  24. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  25.  
  26. % Next command sets the type of pade approximation to be used.
  27. % Set ptype=0 for a smoother step response and low pass 
  28. % filter characteristics, otherwise set ptype=1. 
  29. ptype = 1;
  30.  
  31. % Polynomial series approximation of exp(-sT) (constant first)
  32. series=1;
  33. for i=1:2*n;
  34.         series(i+1)=-series(i)*T/i;
  35. end
  36.  
  37. % DENOMINATOR COEFFICIENTS
  38. % First form a matrix to index appropriate elements of 'series'
  39. C=zeros(n,n+1);
  40. for i=1:n
  41.     C(i,:)=i+ptype:n+i+ptype;
  42. end
  43. C(:)=series(C(:));
  44.  
  45. % Find b, a null vector of C,  so that C*b'=0 
  46. [Q,R,E]=qr(C');
  47. if min(size(R))>1 R=sum(R'); end
  48. [dum,ind]=min(abs(R));
  49. % Denominator solution
  50. b=Q(:,ind).';
  51.  
  52.  
  53. % NUMERATOR COEFFICIENTS
  54. % Indexing matrix
  55. C=zeros(n+1,n+1);
  56. for i=1:n+1
  57.     C(i,:)=-i+2:n+2-i;
  58. end
  59. % Address elements of series only for +ve elements of C
  60. C(:)=series(C(:)+n*(C(:)<1))'.*(C(:)>0);
  61. % Numerator solution
  62. a=b*C.'; 
  63.  
  64. if ~ptype, a(1)=[]; end
  65.  
  66. % Graphical Output if no left hand arguments.
  67. % Step Response and Bode Plots
  68. if nargout==0
  69.     clg
  70.     hold off
  71.     subplot(211)
  72.     t1=[0:T/40:2*T];
  73.     y1=step(a,b,t1);
  74.     plot(t1,y1)
  75.     hold off
  76.     str=['st';'nd';'rd';'th']; str=str((n<4).*n+4*(n>3),:);
  77.     title(['Step Response of ',num2str(n),'''',str,' order Pade Approx'])
  78.     xlabel('Time (secs)')
  79.     ylabel('Amplitude')
  80.     [mag,phase,w]=bode(a,b);
  81.     subplot(223)
  82.         semilogx(w,20*log10(mag))
  83.         xlabel('Frequency (rad/sec)')
  84.         ylabel('Gain dB')
  85.         subplot(224)
  86.         semilogx(w,phase)
  87.         xlabel('Frequency (rad/sec)')
  88.         ylabel('Phase deg')
  89.     subplot(212)
  90.     title('Bode Plot')
  91. end
  92. if nargout==4
  93.     [a,b,c,d]=tf2ss(a,b);
  94. elseif nargout==3
  95.     [z,p,k]=tf2zp(a,b);
  96. end
  97.  
  98.  
  99.  
  100.     
  101.  
  102.