wp=0.2*pi;
  ws=0.3*pi;
  rp=1;
  as=15;
  r=10^(-rp/20);
  a=10^(-as/20);
  n=ceil((log10((10^(rp/10)-1)/(10^(as/10)-1)))/(2*log10(wp/ws)));
  fprintf('\n***Butterworth Filter Order=%2.0f',n)
  omega=wp/((10^(rp/10)-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);
  p=p*omega;
  k=k*omega^n;
  B=real(poly(z));
  b0=k;
  b=k*B;
  a=real(poly(p));
  w1=0:50;
  w=w1*pi/50;
  H=freqs(b,a,w);
  mag=abs(H);
  db=20*log10((mag+eps)/max(mag));
  pha=angle(H);
  [ha,x,t]=impulse(b,a);
  subplot(221)
  plot(w/pi,db);
  title('幅度響應');
  ylabel('|H|');
  subplot(222)
  plot(w/pi,abs(mag));
  title('幅度(dB)');
  ylabel('分貝數');
  subplot(223)
  plot(w/pi,pha);
  title('相位響應');
  ylabel('弧度');
  xlabel('以Pi為單位頻率');
  subplot(224)
  plot(t,ha);
  title('沖激響應');
  ylabel('ha(t)');
  xlabel('時間（秒)');