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