wp=0.2*pi;ws=0.3*pi;fs=50;rp=1;rs=15;
t=1/fs;
op=wp*fs;
os=ws*fs;
[N,Wn]=buttord(op,os,rp,rs,'s');
[Z,P,K]=buttap(N);
[A,B,C,D]=zp2ss(Z,P,K);
[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn);
[num1,den1]=ss2tf(AT,BT,CT,DT);
[num2,den2]=bilinear(num1,den1,fs);
[H,W]=freqz(num2,den2,50);
w=W*fs/pi;
h=20*log10(abs(H));
w=w';
fil=abs(H');
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]
[M,N]=size(xn);
xk=fft(xn,N);
[J,L]=size(fil);
if N>L fil=[fil,zeros(1,N-L)]
  elseif L<N xk=[xk,zeros(1,L-N)]
end
xk_fil=fil.*xk;
n=[0:1:N-1];
figure(1)
subplot(2,1,1);
stem(n,xn);
xlabel('n');
ylabel('xn');
title('原始信號');
k=[0:1:N-1];
subplot(2,1,2);
stem(k,xk);
xlabel('k');
ylabel('xk');
title('原始信號的DFT');
l=[0:1:N-1];
figure(2)
subplot(2,1,1);
stem(l,xk_fil);
xlabel('l');
ylabel('xk_fil');
title('濾波后的頻譜');
xn_fil=ifft(xk_fil,N);
subplot(2,1,2)
stem(n,xn_fil);
xlabel('n');
ylabel('xn_fil');
title('濾波后的信號');