wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15;
T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(wr/2);

%function [b,a]=afd_butt(Omegap,Omegar,Ap,Ar);
if Omegap<=0
    error('通帶邊緣必須大于0')
end
if (Omegar<=Omegap)
    error('阻帶邊緣必須大于通帶邊緣必須')
end
N=ceil((log10((10^(0.1*Ap)-1)/(10^(0.1*Ar)-1)))/(2*log10(Omegap/Omegar)));
fprintf('\n***巴特沃斯模擬低通濾波器階數=%2.0f\n',N)
OmegaC=Omegap/((10^(0.1*Ap)-1)^(1/(2*N)))

 z=[];
  p=exp(i*(pi*(1:2:N-1)/(2*N)+pi/2));
  p=[p;conj(p)];
  p=p(:);
if rem(N,2)==1
    p=[p;-1];
end
  k=real(-p);

%function [b,a]=u_buttap(N,OmegaC);
[z,p,k]=buttap(N);
p=p*OmegaC;
k=k*OmegaC^N;
B=real(poly(z));
b0=k;
b=k*B;
a=real(poly(p));

%function[A,B,C]=sdir2cas(b,a);
Na=length(a)-1;Nb=length(b)-1;
b0=b(1);b=b/b0;
a0=a(1);a=a/a0;
C=b0/a0;
p=cplxpair(roots(a));k=floor(Na/2);
if k*2==Na
    A=zeros(k,3);
    for n=1:2:Na
        Arow=p(n:1:n+1,:);Arow=poly(Arow);
        A(fix((n+1)/2),:)=real(Arow);
    end
elseif Na==1
    A=[0 real(poly(Arow))];
else
    A=zeros(k+1,3);
   for n=1:2:2*k
        Arow=p(n:1:n+1,:);Arow=poly(Arow);
        A(fix((n+1)/2),:)=real(Arow);
    end 
      A(k+1,:)=[0 real(poly(p(Na)))];
end
z=cplxpair(roots(b));k=floor(Nb/2);
if Nb==0
    B=[0 0 poly(z)];
elseif k*2==Nb
    B=zeros(k,3);
     for n=1:2:Nb
        Brow=z(n:1:n+1,:);Brow=poly(Brow);
        B(fix((n+1)/2),:)=real(Brow);
     end
elseif Nb==1
    A=[0 real(poly(z))];
else 
    B=zeros(k+1,3);
     for n=1:2:2*k
        Brow=z(n:1:n+1,:);Brow=poly(Brow);
        B(fix((n+1)/2),:)=real(Brow);
     end 
       B(k+1,:)=[0 real(poly(z(Nb)))];
end  

%function [db,mag,pha,Omega]=freqs_m(b,a,Omega_max);
Omega_max=0:50;
Omega=Omega_max.*pi/50;
H=freqs(b,a,Omega);
mag=abs(H);
dB=20*log10((mag+eps)/max(mag));
pha=angle(H);

subplot(224);plot(Omega/pi,dB);
title('模擬濾波器幅度響應|H|(j/Omegar)');
[b,a]=bilinear(b,a,T);

%function[C,B,A]=dir2cas(b,a);
b0=b(1);b=b/b0;
a0=a(1);a=a/a0;
C=b0/a0;
M=length(a)-1;N=length(b)-1;
p=cplxpair(roots(a));k=floor(Na/2);
if N>M
    b=[b zeros(1,N-M)];
  
elseif M>N
    a=[a zeros(1,M-N)];N=M;  
else
   NM=0;
end
k=floor(N/2);B=zeros(k,3);A=zeros(k,3);
if k*2==N;
   b=[b 0];
   a=[a 0];
end
broots=cplxpair(roots(b));
aroots=cplxpair(roots(a));    

for i=1:2:2*k
     Brow=broots(i:1:i+1,:);Brow=real(poly(Brow));
     B(fix((i+1)/2),:)=Brow;
     Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));
     A(fix((i+1)/2),:)=Arow;
end

%function [db,mag,pha,w]=freqz_m(b,a);
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';
w=(w(1:1:501))';
mag=abs(H)
db=20*log10((mag+eps)/max(mag))
pha=angle(H)

subplot(221);plot(w/pi,mag);title('數字濾波器幅度響應|H|(ej/Omegar)');
subplot(222);plot(w/pi,db);title('數字濾波器幅度響應(db)');
subplot(223);plot(w/pi,pha/pi);title('數字濾波器相位響應');
delta_w=2*pi/1000
Ap=-(min(db(1:1:wp/delta_w+1)))
Ar=-round(max(db(delta_w+1:1:501)))