%-----------------------------------------------------------------
%     *** 2.161 Signal Processing - Continuous and Discrete  ***
%                      Tutorial MATLAB Function
%
%   lpbutter -  Butterworth continuous lowpass filter design based on
%       derivation in the 2.161 class notes.
%
%   Source: Class Handout: Introduction to Continuous Time Filter Design
%
%   Usage:  [z,p,k,n] = lpbutter(wc, Rc, ws, Rs) or
%           [z,p,k,n] = lpbutter(wc, Rc, ws, Rs,'zpk') returns the zeros,
%                        poles, and gain, and order for a Butterworth
%                        continuous lowpass filter.
%           [b,a,n]   = lpbutter(wc, Rc, ws, Rs,'tf') returns the transfer
%                        function numerator and denominator coefficients,
%                        and order for a Butterworth 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] = lpbutter(wc,Rc,ws,Rs,varargin)
%
% Compute the required order n
%
  epsilon = sqrt(10^(Rc/10)-1);
  lambda  = sqrt(10^(Rs/10)-1);
  n       = ceil(log(lambda/epsilon)/log(ws/wc));
%
% There are no zeros in the lpf - create an empty zeros vector:
%
  z = [];
%
% Create the complex conjugate poles on a circle
%
  p = wc*epsilon^(-1/n)*exp(i*(pi*(1:2:n-1)/(2*n) + pi/2));
  p = [p; conj(p)];
  p = p(:);
% If n is odd, add a real pole
  if rem(n,2)==1   % n is odd
    p = [p; -wc*epsilon^(-1/n)];
  end
%
% For unity low freq gain, the k is the product of the negative of the poles.
%
  k = real(prod(-p));
%
% 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