function f = er_lpfc1(N,rip)
% function f = er_lpfc1(N,rip)
%   Function to create prototype Chebyshev Type I digital LPF using bilinear
%   transformation with maximally flat characteristics in the stopband.  Takes
%   filter order in N (even/odd integer) and maximum allowable passband ripple 
%   in RIP (in dB).  The poles will be evenly spaced about an ellipse in the left
%   half plane for the analog prototype.  The cutoff frequency (in rads, corresponding 
%   to filter magnitude response value of 10^(-RIP/20)) is normalized.  Function 
%   returns transfer function coefficient vectors in struct F.
%
%   EXAMPLE: >> f = er_lpfc1(6,0.1), yields 6th-order Type-I Chebyshev
%   prototype digital filter with normalized passband edge at pi/4 with
%   maximum ripple of 0.1 dB in stopband.
%
%   Ref: T.W. Parks and C.S. Burrus, Digital Filter Design, New York: John Wiley & Sons, 1987. Chapter 7.
%
%   Author:  Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB assignment, FA2004

% Filter must be greater than or equal to 2
if (N < 2)
    error('Filter order must be >= 2');
end

% Extract parameters from desired passband ripple
epsilon = sqrt(10^(0.1*rip)-1);
mu = asinh(1/epsilon)/N;

% Generate pole/zero vectors for normalized analog prototype
p1 = exp(j*(pi*(1:2:2*N-1)/(2*N) + pi/2))';
p = sinh(mu)*real(p1) + j*cosh(mu)*imag(p1);
z = []; % No zeros

% Generate gain constant
k = prod(-p);
if (rem(N,2) ~= 0)  % Check and adjust for for even order
	k = k/sqrt((1 + epsilon^2));
end

% Do zero/pole bilinear transformation, T = 2
fs = 1;
pd = (1+p/fs)./(1-p/fs); % Flip for z^(-1)
zd = (1+z/fs)./(1-z/fs);
kd = (k*prod(fs-z)./prod(fs-p));
zd = [zd;-ones(length(pd)-length(zd),1)];  % Add extra zeros at -1
numd = real(kd*poly(zd));
dend = real(poly(pd));

% Pack output parameters into struct
f = struct('tf_complete',{numd,dend});