% Lowpass digital filter with Chebyshev-I analog prototype
%

function lowpass
% Digital Filter Specifications:
wp = 0.125*2*pi; % digital passband frequency in Hz (normalized)
ws = 0.1375*2*pi; % digital stopband frequency in Hz (normalized)
Rp = 0.5; % passband ripple in dB
As = 20; % stopband attenuation in dB
% Analog Prototype Specifications:
Fs = 1; T = 1/Fs;
OmegaP = (2/T)*tan(wp/2); % prewarp prototype passband frequency
OmegaS = (2/T)*tan(ws/2); % prewarp prototype stopband frequency
% Analog Chebyshev-1 Prototype Filter Calculation:
[c, d] = chb1(OmegaP, OmegaS, Rp, As);
% Bilinear Transformation:
[b, a] = bilinear(cs, ds, Fs);
%
[db,mag,pha,grd,w] = freqz(b,a);
plot(w*8000/2/pi,db);
xlabel('frequency (Hz)'); ylabel('decibels'); title('Magnitude in dB');

endfunction

function [b,a] = chb1(Wp, Ws, Rp, As);
  % Analog Lowpass Filter Design: Chebyshev-1
  %
  % [b,a] = chb1(Wp, Ws, Rp, As);
  % b = Numerator coefficients of Ha(s)
  % a = Denominator coefficients of Ha(s)
  % Wp = Passband edge frequency in rad/sec
  % Ws = Stopband edge frequency in rad/sec
  % Rp = Passband ripple in dB
  % As = Stopband attenuation in dB
  %
  if Wp <= 0
  error('Passband edge must be larger than 0')
  end
  if Ws <= Wp
  error('Stopband edge must be larger than Passband edge')
  end
  if (Rp <= 0) | (As < 0)
  error('PB ripple and/or SB attenuation must be larger than 0')
  end
  ep = sqrt(10^(Rp/10)-1);
  A = 10^(As/20);
  OmegaC = Wp;
  OmegaR = Ws/Wp;
  g = sqrt(A*A-1)/ep;
  N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));
  fprintf('\n*** Chebyshev-1 Filter Order = %2.0f \n',N);
  [b,a] = ap_chb1(N, Rp, OmegaC);
endfunction


function [b,a] = ap_chb1(N, Rp, Omegac);
% Chebyshev-1 Analog Lowpass Filter Prototype
%
% [b,a] = ap_chb1(N, Rp, Omegac);
% b = nemerator polynomial coefficients
% a = denominator polynomial coefficients
% N = Order of the elliptical filter
% Rp = Passband Ripple in dB
% Omegac = cutoff frequency in rad/sec
%
[z,p,k] = cheb1ap(N,Rp);
a = real(poly(p));
aNn = a(N+1);
p = p*Omegac;
a = real(poly(p));
aNu = a(N+1);
k = k*aNu/aNn;
b0 = k;
B = real(poly(z));
b = k*B;

endfunction