butterbir(Rp,As,F1,F2,Fs) %[N,omegac,b,G,a1,a2]

function butterbir(Rp,As,F1,F2,Fs) %[N,omegac,b,G,a1,a2]
A1=1/(10^(Rp/20));
A2=1/(10^(As/20));
w1=(2*pi*F1)/Fs;
w2=(2*pi*F2)/Fs;
%Compute Order of filter, N
N=ceil(log((((1/A2)^2)-1)/(((1/A1)^2)-1))/(2*log(tan(w2/2)/tan(w1/2))));
%Compute analog cut-off frequency, omegac
T=1/Fs;
omegac=((2/T)*tan(w1/2))/((((1/A1)^2)-1)^(1/(2*N)));
%Compute coefficients of analog filter tranfer function
if (rem(N,2)==0)
    G=(omegac^2)^(N/2);
    for k=1:(N/2)
        b(k)=2*sin((((2*k)-1)*pi)/(2*N));
        a1(k)=b(k)*omegac;
        a2(k)=omegac^2;
    end
elseif (rem(N,2)~=0)
    G=omegac*((omegac^2)^((N-1)/2));
    a2(1)=omegac;
    for k=1:((N-1)/2)
        b(k)=2*sin((((2*k)-1)*pi)/(2*N));
        a1(k)=b(k)*omegac;
        a2(k+1)=omegac^2;
    end
end
%Perform polynomial multiplication in analog domain
sn=[zeros(1,N),G];
if (rem(N,2)==0)
    for k=1:(N/2)
        sk(k,:)=[1,a1(k),a2(k)];
    end
    sd=sk(1,:);
    if ((N/2)>1)
        for k=2:(N/2)
            sd=conv(sd,sk(k,:));
        end
    else
    end
elseif (rem(N,2)~=0)
    sk(1,:)=[0,1,a2(1)];
    for k=1:((N-1)/2)
        sk(k+1,:)=[1,a1(k),a2(k+1)];
    end
    sd=sk(1,:);
    for k=1:((N-1)/2)
        sd=conv(sd,sk(k+1,:));
    end
    sd=sd(2:length(sd));
end
%Obtain digital filter coefficients
[zn,zd]=bilinear(sn,sd,Fs);
%Display values
disp(['Sampling period, T=',num2str(T),' secs']);
disp(['Order of filter, N=',num2str(N)]);
disp(['Analog cut-off frequency, omegac=',num2str(omegac),' Hz']);
disp(['Gain, G=',num2str(G)]);
disp('Filter coefficients:');
if (rem(N,2)==0)
    for k=1:(N/2)
        disp(['             b',num2str(k),'=',num2str(b(k))]);
        disp(['             a1',num2str(k),'=',num2str(a1(k))]);
        disp(['             a2',num2str(k),'=',num2str(a2(k))]);
    end
elseif (rem(N,2)~=0)
    disp(['             a20=',num2str(a2(1))]);
    for k=1:((N-1)/2)
        disp(['             b',num2str(k),'=',num2str(b(k))]);
        disp(['             a1',num2str(k),'=',num2str(a1(k))]);
        disp(['             a2',num2str(k),'=',num2str(a2(k+1))]);
    end
end
disp('Numerator polynomial coefficients in analog domain:');
disp(['         ',num2str(sn)]);
disp('Denominator polynomial coefficients in analog domain:');
disp(['         ',num2str(sd)]);
disp('Numerator polynomial coefficients in digital domain:');
disp(['         ',num2str(zn)]);
disp('Denominator polynomial coefficients in digital domain:');
disp(['         ',num2str(zd)]);