clear;close all;
pi=3.1416
fs=2*pi*1.5*10000;
fst1= 2*pi*0.8*1000; 
fp1=2*pi*2.3*1000;
fp2=2*pi*4.55*1000; %給定的阻帶的截止頻率
fstr2=2*pi*6.05*1000;%給定的通帶的截止頻率
wls=2*pi*fst1/fs; %低端阻帶截止頻率
wlp=2*pi*fp1/fs; %低端通帶截止頻率
wup=2*pi*fp2/fs;%高端通帶截止頻率
B=wlp-wls; %過渡帶寬
N=ceil(8*pi/B);
n=0:N-1; 
w_han=(hanning(N))';             %求窗函數(shù)
wp=[wlp/pi-6/N,wup/pi+6/N];
hn=fir1(N-1,wp,'bandpass',hanning(N));
n=0:N-1;
b=hn; 
figure(1)
[H,f]=freqz(b,1,512,15000); %采用 15000 Hz 的采樣頻率繪出該濾波器的幅頻和相頻響應(yīng)
subplot(2,1,1),plot(f,20*log10(abs(H)))
axis([0 7500 -100 50])
xlabel('頻率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('頻率/Hz');ylabel('相位/^o');grid on;
%impz(b,1);                    %可采用此函數(shù)給出濾波器的脈沖響應(yīng)
%zplane(b,1);                  %可采用此語句給出濾波器的零極點(diǎn)圖
%grpdelay(b,1);                %可采用此函數(shù)給出濾波器的群延遲
f1=3;f2=20;                   %檢測(cè)輸入信號(hào)含有兩種頻率成分
dt=0.02; t=0:dt:3;    %采樣間隔和檢測(cè)信號(hào)的時(shí)間序列
x=sin(2*pi*f1* t)+cos(2* pi*f2* t);               %檢測(cè)信號(hào)
%y=filter(b,1,x);                    %可采用此函數(shù)給出濾波器的輸出
y=fftfilt(b,x);                      %給出濾波器的輸出
figure(2)
subplot(2,1,1);
n=0:N-1; 
stem(n,hn,'.'); 
axis([0 40 -0.6 0.6]);
xlabel('n');ylabel('h(n)')
grid on;
subplot(2,1,2);stem(n,w_han);title('漢寧窗w(n)')
axis([0 40 0 1.2])

t=(0:100)/fs;
x=sin(2*pi*t*4000)+sin(2*pi*t*20000)+sin(2*pi*t*40000)+sin(2*pi*t*30000);
q=filter(hn,1,x);
[a,f1]=freqz(x);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
figure(3);
subplot(2,1,1);
plot(f1,abs(a));
title('輸入波形頻譜圖');
xlabel('頻率');ylabel('幅度')
subplot(2,1,2);
plot(f2,abs(b));
title('輸出波形頻譜圖');
xlabel('頻率');ylabel('幅度')