home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / CHEB2ORD.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  99 lines
  1. function [order,wn] = cheb2ord(wp,ws,rp,rs,opt)
  2. %CHEB2ORD Chebyshev type II filter order selection.
  3. %    [N, Wn] = CHEB2ORD(Wp, Ws, Rp, Rs) returns the order N of the
  4. %    lowest order digital Chebyshev type I filter that loses no
  5. %    more than Rp dB in the passband and has at least Rs dB of
  6. %    attenuation in the stopband.  The passband runs from 0 to
  7. %    Wp and the stopband extends from Ws to 1.0, the Nyquist
  8. %    frequency.  CHEB2ORD also returns Wn, the Chebyshev 
  9. %    natural frequency to use with CHEBY2 to achieve the 
  10. %    specifications.
  11. %    [N, Wn] = CHEB1ORD(Wp, Ws, Rp, Rs, 's') does the computation
  12. %    for an analog filter.
  13. %    See also CHEBY2, BUTTORD, ELLIPORD, and CHEB1ORD.
  14.  
  15. %    Reference: Rabiner and Gold, p 241.
  16.  
  17. %    L. Shure 6-9-88
  18. %    updated 11-18-92, T. Krauss
  19. %    Copyright (c) 1988, 1992 by the MathWorks, Inc.
  20.  
  21. if nargin == 4
  22.     opt = 'z';
  23. elseif nargin == 5
  24.     if ~strcmp(opt,'z') & ~strcmp(opt,'s')
  25.         error('Invalid option for final argument.');
  26.     end
  27. end
  28. np1 = length(wp);
  29. ns1 = length(ws);
  30. if (np1 ~= ns1)
  31.     error('The frequency vectors must both be the same length.')
  32. end
  33. % figure out filter type
  34. ftype = 2*(np1 - 1);
  35. if wp(1) < ws(1)
  36.     ftype = ftype + 1;    % low (1) or reject (3)
  37. else
  38.     ftype = ftype + 2;    % high (2) or pass (4)
  39. end
  40.  
  41. % first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
  42. if strcmp(opt,'z')    % digital
  43.     WPA=tan(pi*wp/2);
  44.     WSA=tan(pi*ws/2);
  45. else  % don't have to if analog already
  46.     WPA=wp;
  47.     WSA=ws;
  48. end
  49.  
  50. % next, transform to low pass prototype with passband edge of 1 and stopband
  51. % edges determined by the following: (see Rabiner and Gold, p.258)
  52. if ftype == 1    % low
  53.     WA=WSA/WPA;
  54. elseif ftype == 2    % high
  55.     WA=WPA/WSA;
  56. elseif ftype == 3    % stop
  57.     WA=(WSA*(WPA(1)-WPA(2)))./(WSA.^2 - WPA(1)*WPA(2));
  58. elseif ftype == 4    % pass
  59.     WA=(WSA.^2 - WPA(1)*WPA(2))./(WSA*(WPA(1)-WPA(2)));
  60. end
  61.  
  62. % find the minimum order cheby. type 2 filter to meet the more demanding spec:
  63. WA=min(abs(WA));
  64. order=ceil(acosh(sqrt((10^(.1*abs(rs))-1)/(10^(.1*abs(rp))-1)))/acosh(WA));
  65. % ref: M.E. Van Valkenburg, "Analog Filter Design", p.232, eqn 8.39
  66.  
  67. % next find the frequency "new_wp" at which the response of an analog low-pass
  68. % prototype is exactly -rp dB.  The prototype is the one for which the beginning
  69. % of the stop-band is at frequency 1.
  70. % (new_wp will be less than 1):
  71. new_wp=1/cosh(1/order*acosh(sqrt((10^(.1*abs(rs)) - 1)/(10^(.1*abs(rp)) - 1))));
  72.  
  73. % Now convert the stop band frequency back from lowpass prototype 
  74. % to the original analog filter.
  75. % Here we use the mapping which maps the frequency "new_wp" to the original WP, 
  76. % to map the (+/-)1 frequency to WN (WN will be between WP and WS):
  77. if ftype == 1    % low
  78.     WN=WPA/new_wp;
  79. elseif ftype == 2    % high
  80.     WN=WPA*new_wp;
  81. elseif ftype == 3    % stop
  82.     WN=(WPA(1)-WPA(2))*new_wp/2 + ...
  83.         sqrt( (WPA(2)-WPA(1))^2*new_wp^2/4 + WPA(1)*WPA(2));
  84.     WN(2)=WPA(1)*WPA(2)/WN(1);
  85. elseif ftype == 4    % pass
  86.     WN=(WPA(1)-WPA(2))/(2*new_wp) + ...
  87.         sqrt( (WPA(2)-WPA(1))^2/(4*new_wp^2) + WPA(1)*WPA(2));
  88.     WN(2)=WPA(1)*WPA(2)/WN(1);
  89. %  WA=(WP.^2 - WN(1)*WN(2))./(WP*(WN(2)-WN(1))) <--- to check, this should be
  90. %        -/+ new_wp
  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.