clear,close all;
N=input('輸入FIR濾波器的長度N=');
wp=input('輸入FIR濾波器的通帶截止頻率（單位：弧度）wp=');
ws=input('輸入FIR濾波器的阻帶截止頻率（單位：弧度） wp=');
wc=(wp+ws)/2;               %確定理想低通濾波器的截止頻率
L=ceil(wc*N/(2*pi));  
k=0:N-1;  M=2*pi/N;
H=u(k)-u(k-L);               %構造FIR濾波器的N抽樣的頻率響應特性
phai=-pi*k*(N-1)/N;
%H(L+1)=0.5;%H(L+2)=0.1065;%H(L+2)=0.031;   %用于插入過渡采樣點
H=H.*exp(i*phai);
kc=fix(N/2);
for i=0:kc;            
    H(N-(i))=conj(H(i+2));
end 
h=ifft(H,N);                   %求DFT反變換得到單位沖激響應h(n)
subplot(221)                   %繪制FIR濾波器的幅度特性
k=0:2*pi/N:2*pi-2*pi/N;
stem(k/pi,abs(H),'k.')
axis([0,2,0,max(abs(H))]);  grid on,  xlabel('Frequency *pi')
title('The |H(k)|of the FIR lowpass filter')
subplot(222)                     %繪制單位沖激響應波形
k=0:N-1;
stem(k,h,'k.'),axis([0,N-1,1.1*min(real(h)),1.1*max(real(h))])
xlabel('Time index n'),  title('The impulse response h(n) of the FIR filter'),  grid on
n=0:2*pi/512:2*pi-2*pi/512;      %計算FIR濾波器的幅度響應
subplot(223)                     %繪制FIR濾波器的幅度響應
Hd=fft(h,512);                 
plot(n/(pi),abs(Hd)),title('The frequency response |H(ejw)|')
hold on,grid on
k=0:2*pi/N:2*pi-2*pi/N;
plot(k/pi,abs(H),'k.'),axis([0,1,0,1.2*max(abs(H))])
xlabel('Frequency *pi')
subplot(224)                   %繪制FIR濾波器衰的減幅度特性曲線
plot(n/pi,20*log10(abs(Hd)/max(abs(Hd))),'r-')
title('The frequency response in dB')
axis([0,1,-80,0])
grid on
xlabel('Frequency *pi')