%-----------------------------------------------------------------
%     *** 2.161 Signal Processing - Continuous and Discrete  ***
%                      Tutorial MATLAB Function
%
%   lpcheby2 -  Chebyshev Type 2 continuous lowpass filter design.
%
%   Source: Class Handout: Introduction to Continuous Time Filter Design
%
%   Usage:  [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs) or
%           [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs,'zpk') returns the zeros,
%                        poles, gain and order for a Chebyshev Type 1
%                        continuous lowpass filter.
%           [b,a,n]   = lpcheby2(wc, Rc, ws, Rs,'tf') returns the transfer
%                        function numerator and denominator coefficients,
%                        and order for a Chebyshev Type 1 continuous lowpass
%                        filter.
%           where wc  - passband cut-off frequency (rad/s)
%                 Rc  - attenuation at passband cut-off (dB) (Rc > 0)
%                 ws  - stopband frequency (rad/s)
%                 Rs  - attenuation at stopband (dB) (Rs > 0).
%
%   Author:   D. Rowell
%   Revision: 2.0 9-23-2007
%
%---------------------------------------------------------------
%
function [varargout] = lpcheby2(wc,Rc, ws, Rs,varargin)
%
% Compute the required order n
%
  epsilon = sqrt(10^(Rc/10)-1);
  lambda  = sqrt(10^(Rs/10)-1);
  n       = ceil(acosh(lambda/epsilon)/acosh(ws/wc));
  wprod   = wc*ws;
  epsilonhat= 1/(epsilon*cos(n*acos(ws/wc)));
%
% Create the finite (imaginary) zeros
% Handle even and odd orders separately to eliminate the zero at infinity
% for n odd
if (rem(n,2))
	m = n - 1;
	z = cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))';
else
	m = n;
	z = cos((1:2:2*n-1)*pi/(2*n))';
end
z = (z - flipud(z))./2;
% Transform back to the s-plane
z = i*ws./z;
% Organize zeros in complex pairs:
j = [1:m/2; m:-1:m/2+1];
z = z(j(:));
%
% Determine the poles - choosing first only complex conjugates in the
% left-half plane
%
alpha = asinh(1/epsilonhat)/n;
p = [];
for j=1:floor(n/2)
   gamma = (2*j - 1)*pi/(2*n);
   newp = -wc*(sinh(alpha)*sin(gamma)+ i*cosh(alpha)*cos(gamma));
   p = [p; newp; conj(newp)];
end
% If n is odd - add a real pole
  if rem(n,2)==1   % n is odd
     gamma = (2*(floor(n/2)+1) - 1)*pi/(2*n);
     p = [p; -wc*(sinh(alpha)*sin(gamma))];
  end
  p=wprod./p;
%
% compute the gain
%
k = real(prod(-p)/prod(-z));
%
% Return the appropriate system representation
%
  if (nargin == 5 )
      if (varargin{1}(1:2) == 'tf')
         sys = zpk(z,p,k);
         [b,a] = tfdata(sys);
         varargout(1) = b;
         varargout(2) = a;
         varargout(3) = {n};
      elseif (varargin{1}(1:3) == 'zpk')
         varargout(1)={z};
         varargout(2)={p};
         varargout(3)={k};
         varargout(4)={n};
      end
  else
      varargout(1)={z};
      varargout(2)={p};
      varargout(3)={k};
      varargout(4)={n};
  end

