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

% choose any value of the passband frequency 
% of the original lowpass filter 
wclp=0.2;

N=M/2;
% translate  wclp to wplp
OmegaC = tan(pi*wclp/2); 
OmegaP =OmegaC/cosh(1/N*acosh(1/sqrt(10^(Ap/10)-1)));
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 method to design filter
[b,a]=cheby1(N,Ap,[wp1 wp2],'stop')
fvtool(b,a);