%-----------------------------------------------------------------
%     *** 2.161 Signal Processing - Continuous and Discrete  ***
%                      Tutorial MATLAB Function
%
%   lpcheby1 -  Chebyshev Type 1 continuous lowpass filter design.
%
%   Source: Class Handout: Introduction to Continuous Time Filter Design
%
%   Usage:  [z,p,k,n] = lpcheby1(wc, Rc, ws, Rs) or
%           [z,p,k,n] = lpcheby1(wc, Rc, ws, Rs,'zpk') returns the zeros,
%                        poles, gain and order for a Chebyshev Type 1
%                        continuous lowpass filter.
%           [b,a,n]   = lpcheby1(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] = lpcheby1(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));
%
% There are no zeros in the lpf - create an empty zeros vector:
%
  z = [];
%
% Determine the poles - choosing first only complex conjugates in the
% left-half plane
%
alpha = asinh(1/epsilon)/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
%
% compute the gain
%
k = real(prod(-p));
%
% If n is even adjust the gain to unity at zero frequency
%
if ~rem(n,2)	% n is even so patch k
	k = k/sqrt((1 + epsilon^2));
end
%
% 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