home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / ELLIPAP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.4 KB  |  78 lines
  1. function [z,p,k] = ellipap(n, rp, rs)
  2. %ELLIPAP Elliptic analog lowpass filter prototype.
  3. %    [Z,P,K] = ELLIPAP(N,Rp,Rs) returns the zeros, poles, and gain
  4. %    of an N-th order normalized prototype elliptic analog lowpass
  5. %    filter with Rp decibels of ripple in the passband and a
  6. %    stopband Rs decibels down.
  7.  
  8. %    L. Shure 3-8-88
  9. %    (c) Copyright 1988, by The MathWorks, Inc.
  10.  
  11. %    References:
  12. %      [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
  13. %          John Wiley & Sons, 1987, chapter 7, section 7.3.7-8.
  14.  
  15. epsilon = sqrt(10^(0.1*rp)-1);    %rp, dB of passband ripple
  16. k1 = epsilon/sqrt(10^(0.1*rs)-1);    %rs, dB of stopband ripple
  17. k1p = sqrt(1-k1^2);
  18. wp = 1;    % passband edge - normalized
  19. if abs(1-k1p^2) < eps
  20.     krat = 0;
  21. else
  22.     capk1 = ellipke([k1^2,k1p^2]);
  23.     krat = n*capk1(1)/capk1(2);    % krat = K(k)/K'(k) -- need to find relevant k
  24. end
  25.  
  26. % try to find m, elliptic parameter, so that K(m)/K'(m) = krat -- K, complete
  27. % elliptic integral of first kind
  28. global ELLIP_KRAT
  29. ELLIP_KRAT = krat;
  30. m = fmins('kratio',.5);
  31.  
  32. capk = ellipke(m);
  33. ws = wp/sqrt(m);    % stopband edge (=> transition band is ws-wp in width)
  34. m1 = 1 - m;
  35. i = sqrt(-1);
  36.  
  37. % find zeros; they are purely imaginary and paired in complex conjugates
  38. j = (1-rem(n,2)):2:n-1;
  39. [ij,jj] = size(j);
  40. % s is Jacobi elliptic function sn(u)
  41. [s,c,d] = ellipj(j*capk/n,m*ones(ij,jj));
  42. is = find(abs(s) > eps);
  43. z = 1 ./(sqrt(m)*s(is));
  44. z = i*z(:);    % make column vector
  45. z = [z ; conj(z)];
  46. % order the zeros for convenience later on
  47. z = cplxpair(z);
  48.  
  49. % poles; one purely real if n is odd - the remainder are complex conjugate pairs
  50. % put 1/epsilon, mp into global variables for calculating v0
  51. global ELLIP_EPS ELLIP_MP;
  52. ELLIP_EPS = 1 / epsilon;
  53. ELLIP_MP = k1p^2;
  54. % calculate v0, a 'fundamental' parameter for the poles related to inverse sc
  55. % function . I.e. find r so sn(r)/cn(r) = 1/epsilon for the given parameter mp
  56. r = fmins('vratio', ellipke(1-m));
  57. v0 = capk*r/(n*capk1(1));
  58. [sv,cv,dv] = ellipj(v0,1-m);
  59. p = -(c.*d*sv*cv + i*s*dv)./(1-(d*sv).^2);
  60. p = p(:);   % make column vector
  61. % check to see if there is a real pole
  62. if rem(n,2)
  63.     ip = find(abs(imag(p)) < eps*norm(p));
  64.     [pm,pn] = size(p);
  65.     pp = 1:pm;
  66.     pp(ip) = [];
  67.     p = [p ; conj(p(pp))];
  68. else
  69.     p = [p; conj(p)];
  70. end
  71. p = cplxpair(p);    % order poles for later use
  72.  
  73. % gain
  74. k = real(prod(-p)/prod(-z));
  75. if (~rem(n,2))    % n is even order so patch gain
  76.     k = k/sqrt((1 + epsilon^2));
  77. end
  78.