home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / BUTTORD.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.1 KB  |  99 lines
  1. function [order,wn] = buttord(wp,ws,rp,rs,opt)
  2. %BUTTORD Butterworth filter order selection.
  3. %    [N, Wn] = BUTTORD(Wp, Ws, Rp, Rs) returns the order N of the
  4. %    lowest order digital Butterworth filter that loses no more than
  5. %    Rp dB in the passband and has at least Rs dB of attenuation
  6. %    in the stopband.  The passband runs from 0 to Wp and the
  7. %    stopband extends from Ws to 1.0, the Nyquist frequency.
  8. %    BUTTORD also returns Wn, the Butterworth natural frequency (or, 
  9. %    the "3 dB frequency") to use with BUTTER to achieve the specifications. 
  10. %    [N, Wn] = BUTTORD(Wp, Ws, Rp, Rs, 's') does the computation
  11. %    for an analog filter.
  12. %
  13. %    When Rp is chosen as 3 dB, the Wn in BUTTER is equal to Wp
  14. %    in BUTTORD.  
  15. %
  16. %    See also BUTTER, CHEB1ORD, CHEB2ORD, and ELLIPORD.
  17. % 
  18.  
  19. %    Reference: Rabiner and Gold, p 241.
  20.  
  21. %    L. Shure 6-9-88
  22. %    updated 11-13-92, T. Krauss
  23. %    Copyright (c) 1988, 1992 by the MathWorks, Inc.
  24.  
  25. if nargin == 4
  26.     opt = 'z';
  27. elseif nargin == 5
  28.     if ~strcmp(opt,'z') & ~strcmp(opt,'s')
  29.         error('Invalid option for final argument.');
  30.     end
  31. end
  32. np1 = length(wp);
  33. ns1 = length(ws);
  34. if (np1 ~= ns1)
  35.     error('The frequency vectors must both be the same length.')
  36. end
  37. ftype = 2*(np1 - 1);
  38. if wp(1) < ws(1)
  39.     ftype = ftype + 1;    % low (1) or reject (3)
  40. else
  41.     ftype = ftype + 2;    % high (2) or pass (4)
  42. end
  43.  
  44. % first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
  45. if strcmp(opt,'z')    % digital
  46.     WP=tan(pi*wp/2);
  47.     WS=tan(pi*ws/2);
  48. else  % don't have to if analog already
  49.     WP=wp;
  50.     WS=ws;
  51. end
  52. %note - on old systems that are NOT case sensitive, this will still work OK
  53.  
  54. % next, transform to low pass prototype with passband edge of 1 and stopband
  55. % edges determined by the following: (see Rabiner and Gold, p.258)
  56. if ftype == 1    % low
  57.     WA=WS/WP;
  58. elseif ftype == 2    % high
  59.     WA=WP/WS;
  60. elseif ftype == 3    % stop
  61.     WA=(WS*(WP(1)-WP(2)))./(WS.^2 - WP(1)*WP(2));
  62. elseif ftype == 4    % pass
  63.     WA=(WS.^2 - WP(1)*WP(2))./(WS*(WP(1)-WP(2)));
  64. end
  65.  
  66. % find the minimum order b'worth filter to meet the more demanding spec:
  67. WA=min(abs(WA));
  68. order = ceil( log10( (10 .^ (0.1*abs(rs)) - 1)./ ...
  69.     (10 .^ (0.1*abs(rp)) - 1) ) / (2*log10(WA)) );
  70.  
  71. % next find the butterworth natural frequency W0 (or, the "3dB frequency")
  72. % to give exactly rs dB at WA.  W0 will be between 1 and WA:
  73. W0 = WA / ( (10^(.1*abs(rs)) - 1)^(1/(2*order)) );
  74.  
  75. % now convert this frequency back from lowpass prototype 
  76. % to the original analog filter:
  77. if ftype == 1    % low
  78.     WN=W0*WP;
  79. elseif ftype == 2    % high
  80.     WN=WP/W0;
  81. elseif ftype == 3    % stop
  82.     WN(1) = ( (WP(2)-WP(1)) + sqrt((WP(2)-WP(1))^2 + ...
  83.         4*W0.^2*WP(1)*WP(2)))./(2*W0);
  84.     WN(2) = ( (WP(2)-WP(1)) - sqrt((WP(2)-WP(1))^2 + ...
  85.         4*W0.^2*WP(1)*WP(2)))./(2*W0);
  86.     WN=sort(abs(WN)); 
  87. elseif ftype == 4    % pass
  88.     W0=[-W0 W0];  % need both left and right 3dB frequencies
  89.     WN= -W0*(WP(2)-WP(1))/2 + sqrt( W0.^2/4*(WP(2)-WP(1))^2 + WP(1)*WP(2) );
  90.     WN=sort(abs(WN)); 
  91. end
  92.  
  93. % finally, transform frequencies from analog to digital if necessary:
  94. if strcmp(opt,'z')    % digital
  95.     wn=(2/pi)*atan(WN);  % bilinear transform
  96. else
  97.     wn=WN;
  98. end
  99.