fir.m

基于汉宁窗的fir带通滤波器 用matlab仿真

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))';             %求窗函数
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 的采样频率绘出该滤波器的幅频和相频响应
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);                    %可采用此函数给出滤波器的脉冲响应
%zplane(b,1);                  %可采用此语句给出滤波器的零极点图
%grpdelay(b,1);                %可采用此函数给出滤波器的群延迟
f1=3;f2=20;                   %检测输入信号含有两种频率成分
dt=0.02; t=0:dt:3;    %采样间隔和检测信号的时间序列
x=sin(2*pi*f1* t)+cos(2* pi*f2* t);               %检测信号
%y=filter(b,1,x);                    %可采用此函数给出滤波器的输出
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('幅度')