function [b,a] = cheby2_3db_bs(M,wc1,wc2,As)
% Inputs:
% N is the order of filter
% wc1 is the low 3-db cutoff frequency   
% wc2 is the up 3-db cutoff frequency   
% As is the attenuation of stopband

% choose any value of the passband frequency 
% of the original lowpass filter 

% 任意選取wclp的值
wclp=0.7;

% halve the order
N=M/2;

% translate  wclp to wplp
OmegaC = tan(pi*wclp/2); 
ep = sqrt(10^(As/10)-1);
OmegaP = OmegaC*cosh(1/N*acosh(ep));
wplp = 2*atan(OmegaP)/pi


% figure out wp1 and wp2
B=cos((wc1+wc2)*pi/2)/cos((wc2-wc1)*pi/2);
K=tan((wc2-wc1)*pi/2)*tan(wclp*pi/2)
alfa1=-2*B/(K+1);
alfa2=(1-K)/(K+1);
B=-exp(-j*wplp*pi);
delta=(B*alfa1+alfa1)^2-4*(alfa2+B)*(1+B*alfa2);
Z2=((alfa1+B*alfa1)+sqrt(delta))/(2*(1+B*alfa2));
Z1=((alfa1+B*alfa1)-sqrt(delta))/(2*(1+B*alfa2));
wp1=-angle(Z1)/pi;
wp2=-angle(Z2)/pi;
wp1=1+wp1
wp2=1-wp2




% use classic methods to design filter
[b,a]=cheby2(N,As,[wp1 wp2],'stop');
fvtool(b,a);
axis(0,1,-5,1)