home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / BUTTER.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.5 KB  |  96 lines
  1. function [num, den, z, p] = butter(n, Wn, ftype)
  2. %BUTTER    Butterworth digital filter design.
  3. %    [B,A] = BUTTER(N,Wn) designs an N'th order lowpass digital
  4. %    Butterworth filter and returns the filter coefficients in length
  5. %    N+1 vectors B and A.  The cut-off frequency Wn must be
  6. %    0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate.
  7. %
  8. %    If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an 
  9. %    order 2N bandpass filter with passband  W1 < W < W2.
  10. %    [B,A] = BUTTER(N,Wn,'high') designs a highpass filter.
  11. %    [B,A] = BUTTER(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
  12. %    
  13. %    When used with three left-hand arguments, as in
  14. %    [Z,P,K] = BUTTER(...), the zeros and poles are returned in
  15. %    length N column vectors Z and P, and the gain in scalar K. 
  16. %
  17. %    When used with four left-hand arguments, as in
  18. %    [A,B,C,D] = BUTTER(...), state-space matrices are returned.
  19. %
  20. %    See also BUTTORD, CHEBY1, CHEBY2, FREQZ and FILTER.
  21.  
  22. %    J.N. Little 1-14-87
  23. %    Revised 1-14-88 JNL, 4-29-88 LS
  24. %    (c) Copyright 1987-88, by The MathWorks, Inc.
  25.  
  26. %    References:
  27. %      [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
  28. %          John Wiley & Sons, 1987, chapter 7, section 7.3.3.
  29.  
  30. btype = 1;
  31. if nargin == 3
  32.     btype = 3;
  33. end
  34. if max(size(Wn)) == 2
  35.     btype = btype + 1;
  36. end
  37.  
  38. % step 1: get analog, pre-warped frequencies
  39. fs = 2;
  40. u = 2*fs*tan(pi*Wn/fs);
  41.  
  42. % step 2: convert to low-pass prototype estimate
  43. if btype == 1    % lowpass
  44.     Wn = u;
  45. elseif btype == 2    % bandpass
  46.     Bw = u(2) - u(1);
  47.     Wn = sqrt(u(1)*u(2));    % center frequency
  48. elseif btype == 3    % highpass
  49.     Wn = u;
  50. elseif btype == 4    % bandstop
  51.     Bw = u(2) - u(1);
  52.     Wn = sqrt(u(1)*u(2));    % center frequency
  53. end
  54.  
  55. % step 3: Get N-th order Butterworth analog lowpass prototype
  56. [z,p,k] = buttap(n);
  57.  
  58. % Transform to state-space
  59. [a,b,c,d] = zp2ss(z,p,k);
  60.  
  61. % step 4: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
  62. if btype == 1        % Lowpass
  63.     [a,b,c,d] = lp2lp(a,b,c,d,Wn);
  64.  
  65. elseif btype == 2    % Bandpass
  66.     [a,b,c,d] = lp2bp(a,b,c,d,Wn,Bw);
  67.  
  68. elseif btype == 3    % Highpass
  69.     [a,b,c,d] = lp2hp(a,b,c,d,Wn);
  70.  
  71. elseif btype == 4    % Bandstop
  72.     [a,b,c,d] = lp2bs(a,b,c,d,Wn,Bw);
  73. end
  74.  
  75. % step 5: Use Bilinear transformation to find discrete equivalent:
  76. [a,b,c,d] = bilinear(a,b,c,d,fs);
  77.  
  78. if nargout == 4
  79.     num = a;
  80.     den = b;
  81.     z = c;
  82.     p = d;
  83. else    % nargout <= 3
  84. % Transform to zero-pole-gain and polynomial forms:
  85.     if nargout == 3
  86.         [z,p,k] = ss2zp(a,b,c,d,1);
  87.         num = z;
  88.         den = p;
  89.         z = k;
  90.     else % nargout <= 2
  91.         den = poly(a);
  92.         num = poly(a-b*c)+(d-1)*den;
  93.     end
  94. end
  95.  
  96.