home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / ELLIPORD.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  72 lines
  1. function [order,wn] = ellipord(wp,ws,rp,rs,opt)
  2. %ELLIPORD Elliptic filter order selection.
  3. %    [N, Wn] = ELLIPORD(Wp, Ws, Rp, Rs) returns the order N of the
  4. %    lowest order digital elliptic 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. %    ELLIPORD also returns Wn, the elliptic natural frequency to
  9. %    use with ELLIP to achieve the specifications.
  10. %    [N, Wn] = ELLIPORD(Wp, Ws, Rp, Rs, 's') does the computation
  11. %    for an analog filter.
  12. %    See also ELLIP, BUTTORD, CHEB1ORD, and CHEB2ORD.
  13.  
  14. %    Reference: Rabiner and Gold, p 241.
  15.  
  16. %    L. Shure 6-9-88
  17. %    updated 11-18-92, T. Krauss
  18. %    Copyright (c) 1988, 1992 by the MathWorks, Inc.
  19.  
  20. if nargin == 4
  21.     opt = 'z';
  22. elseif nargin == 5
  23.     if ~strcmp(opt,'z') & ~strcmp(opt,'s')
  24.         error('Invalid option for final argument.');
  25.     end
  26. end
  27. np1 = length(wp);
  28. ns1 = length(ws);
  29. if (np1 ~= ns1)
  30.     error('The frequency vectors must both be the same length.')
  31. end
  32. ftype = 2*(np1 - 1);
  33. if wp(1) < ws(1)
  34.     ftype = ftype + 1;    % low (1) or reject (3)
  35. else
  36.     ftype = ftype + 2;    % high (2) or pass (4)
  37. end
  38.  
  39. % natural frequencies are simply the given band-edges:
  40. wn=wp;
  41.  
  42. % first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
  43. if strcmp(opt,'z')    % digital
  44.     WP=tan(pi*wp/2);
  45.     WS=tan(pi*ws/2);
  46. else  % don't have to if analog already
  47.     WP=wp;
  48.     WS=ws;
  49. end
  50.  
  51. % next, transform to low pass prototype with passband edge of 1 and stopband
  52. % edges determined by the following: (see Rabiner and Gold, p.258)
  53. if ftype == 1    % low
  54.     WA=WS/WP;
  55. elseif ftype == 2    % high
  56.     WA=WP/WS;
  57. elseif ftype == 3    % stop
  58.     WA=(WS*(WP(1)-WP(2)))./(WS.^2 - WP(1)*WP(2));
  59. elseif ftype == 4    % pass
  60.     WA=(WS.^2 - WP(1)*WP(2))./(WS*(WP(1)-WP(2)));
  61. end
  62.  
  63. % find the minimum order elliptic filter to meet the more demanding spec:
  64. WA = min(abs(WA));
  65. epsilon = sqrt(10 .^ (0.1*rp)-1);
  66. k1 = epsilon./sqrt(10 .^ (0.1*rs)-1);
  67. k = 1/WA;
  68. capk = ellipke([k.^2 1-k.^2]);
  69. capk1 = ellipke([(k1 .^2) 1-(k1 .^2)]);
  70. order = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));
  71.  
  72.