home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / BILINEAR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.0 KB  |  95 lines
  1. function [zd, pd, kd, dd] = bilinear(z, p, k, fs, fp, fp1)
  2. %BILINEAR Bilinear transformation with optional frequency prewarping.
  3. %    [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer
  4. %    function specified by Z, P, and K to a z-transform discrete
  5. %    equivalent obtained from the bilinear transformation:
  6. %
  7. %       H(z) = H(s) |
  8. %                   | s = 2*Fs*(z-1)/(z+1)
  9. %
  10. %    where column vectors Z and P specify the zeros and poles, scalar
  11. %    K specifies the gain, and Fs is the sample frequency in Hz.
  12. %    [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are 
  13. %    row vectors containing numerator and denominator transfer
  14. %    function coefficients, NUM(s)/DEN(s), in descending powers of
  15. %    s, transforms to z-transform coefficients NUMd(z)/DENd(z).
  16. %    [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version.
  17. %    Each of the above three forms of BILINEAR accepts an optional
  18. %    additional input argument that specifies prewarping. For example,
  19. %    [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping before
  20. %    the bilinear transformation so that the frequency responses
  21. %    before and after mapping match exactly at frequency point Fp
  22. %    (match point Fp is specified in Hz).
  23.  
  24. %    J.N. Little 4-28-87 
  25. %    Revised 5-5-87 JNL
  26. %    Copyright (c) 1987 by the MathWorks, Inc.
  27.  
  28. %    Gene Franklin, Stanford Univ., motivated the state-space
  29. %    approach to the bilinear transformation.
  30.  
  31. [mn,nn] = size(z);
  32. [md,nd] = size(p);
  33.  
  34. if (nd == 1 & nn < 2) & nargout ~= 4    % In zero-pole-gain form
  35.     if mn > md
  36.         error('Numerator cannot be higher order than denominator.')
  37.     end
  38.     if nargin == 5        % Prewarp
  39.         fp = 2*pi*fp;
  40.         fs = fp/tan(fp/fs/2);
  41.     else
  42.         fs = 2*fs;
  43.     end
  44.     z = z(finite(z));     % Strip infinities from zeros
  45.     pd = (1+p/fs)./(1-p/fs); % Do bilinear transformation
  46.     zd = (1+z/fs)./(1-z/fs);
  47.     kd = real(k*prod(fs-z)./prod(fs-p));
  48.     zd = [zd;-ones(length(pd)-length(zd),1)];  % Add extra zeros at -1
  49.  
  50. elseif (md == 1 & mn == 1) | nargout == 4 %
  51.     if nargout == 4        % State-space case
  52.         a = z; b = p; c = k; d = fs; fs = fp;
  53.         error(abcdchk(a,b,c,d));
  54.         if nargin == 6            % Prewarp
  55.             fp = fp1;        % Decode arguments
  56.             fp = 2*pi*fp;
  57.             fs = fp/tan(fp/fs/2)/2;
  58.         end
  59.     else            % Transfer function case
  60.         if nn > nd
  61.             error('Numerator cannot be higher order than denominator.')
  62.         end
  63.         num = z; den = p;        % Decode arguments
  64.         if nargin == 4            % Prewarp
  65.             fp = fs; fs = k;    % Decode arguments
  66.             fp = 2*pi*fp;
  67.             fs = fp/tan(fp/fs/2)/2;
  68.         else
  69.             fs = k;            % Decode arguments
  70.         end
  71.         % Put num(s)/den(s) in state-space canonical form.  
  72.         [a,b,c,d] = tf2ss(num,den);
  73.     end
  74.     % Now do state-space version of bilinear transformation:
  75.     n = max(size(a));
  76.     t = 1/fs;
  77.     r = sqrt(t);
  78.     t1 = eye(n) + a*t/2;
  79.     t2 = eye(n) - a*t/2;
  80.     ad = t2\t1;
  81.     bd = t/r*(t2\b);
  82.     cd = r*c/t2;
  83.     dd = c/t2*b*t/2 + d;
  84.     if nargout == 4
  85.         zd = ad; pd = bd; kd = cd;
  86.     else
  87.         % Convert back to transfer function form:
  88.         p = poly(ad);
  89.         zd = poly(ad-bd*cd)+(dd-1)*p;
  90.         pd = p;
  91.     end
  92. else
  93.     error('First two arguments must have the same orientation.')
  94. end
  95.