function [order,wn] = cheb2ord(wp,ws,rp,rs,opt)
%CHEB2ORD Chebyshev type II filter order selection.
%	[N, Wn] = CHEB2ORD(Wp, Ws, Rp, Rs) returns the order N of the
%	lowest order digital Chebyshev type I filter that loses no
%	more than Rp dB in the passband and has at least Rs dB of
%	attenuation in the stopband.  The passband runs from 0 to
%	Wp and the stopband extends from Ws to 1.0, the Nyquist
%	frequency.  CHEB2ORD also returns Wn, the Chebyshev 
%	natural frequency to use with CHEBY2 to achieve the 
%	specifications.
%	[N, Wn] = CHEB1ORD(Wp, Ws, Rp, Rs, 's') does the computation
%	for an analog filter.
%	See also CHEBY2, BUTTORD, ELLIPORD, and CHEB1ORD.

%	Reference: Rabiner and Gold, p 241.

%	L. Shure 6-9-88
%	updated 11-18-92, T. Krauss
%	Copyright (c) 1988, 1992 by the MathWorks, Inc.

if nargin == 4
	opt = 'z';
elseif nargin == 5
	if ~strcmp(opt,'z') & ~strcmp(opt,'s')
		error('Invalid option for final argument.');
	end
end
np1 = length(wp);
ns1 = length(ws);
if (np1 ~= ns1)
	error('The frequency vectors must both be the same length.')
end
% figure out filter type
ftype = 2*(np1 - 1);
if wp(1) < ws(1)
	ftype = ftype + 1;	% low (1) or reject (3)
else
	ftype = ftype + 2;	% high (2) or pass (4)
end

% first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
if strcmp(opt,'z')	% digital
	WPA=tan(pi*wp/2);
	WSA=tan(pi*ws/2);
else  % don't have to if analog already
	WPA=wp;
	WSA=ws;
end

% next, transform to low pass prototype with passband edge of 1 and stopband
% edges determined by the following: (see Rabiner and Gold, p.258)
if ftype == 1	% low
	WA=WSA/WPA;
elseif ftype == 2	% high
	WA=WPA/WSA;
elseif ftype == 3	% stop
	WA=(WSA*(WPA(1)-WPA(2)))./(WSA.^2 - WPA(1)*WPA(2));
elseif ftype == 4	% pass
	WA=(WSA.^2 - WPA(1)*WPA(2))./(WSA*(WPA(1)-WPA(2)));
end

% find the minimum order cheby. type 2 filter to meet the more demanding spec:
WA=min(abs(WA));
order=ceil(acosh(sqrt((10^(.1*abs(rs))-1)/(10^(.1*abs(rp))-1)))/acosh(WA));
% ref: M.E. Van Valkenburg, "Analog Filter Design", p.232, eqn 8.39

% next find the frequency "new_wp" at which the response of an analog low-pass
% prototype is exactly -rp dB.  The prototype is the one for which the beginning
% of the stop-band is at frequency 1.
% (new_wp will be less than 1):
new_wp=1/cosh(1/order*acosh(sqrt((10^(.1*abs(rs)) - 1)/(10^(.1*abs(rp)) - 1))));

% Now convert the stop band frequency back from lowpass prototype 
% to the original analog filter.
% Here we use the mapping which maps the frequency "new_wp" to the original WP, 
% to map the (+/-)1 frequency to WN (WN will be between WP and WS):
if ftype == 1	% low
	WN=WPA/new_wp;
elseif ftype == 2	% high
	WN=WPA*new_wp;
elseif ftype == 3	% stop
	WN=(WPA(1)-WPA(2))*new_wp/2 + ...
		sqrt( (WPA(2)-WPA(1))^2*new_wp^2/4 + WPA(1)*WPA(2));
	WN(2)=WPA(1)*WPA(2)/WN(1);
elseif ftype == 4	% pass
	WN=(WPA(1)-WPA(2))/(2*new_wp) + ...
		sqrt( (WPA(2)-WPA(1))^2/(4*new_wp^2) + WPA(1)*WPA(2));
	WN(2)=WPA(1)*WPA(2)/WN(1);
%  WA=(WP.^2 - WN(1)*WN(2))./(WP*(WN(2)-WN(1))) <--- to check, this should be
%		-/+ new_wp
end

% finally, transform frequencies from analog to digital if necessary:
if strcmp(opt,'z')	% digital
	wn=(2/pi)*atan(WN);  % bilinear transform
else
	wn=WN;
end
