实验目的

1. 加深对窗函数法设计FIR数字滤波器基本方法的理解。
2. 学会用MATLAB语言的窗函数法设计FIR数字滤波器的程序思想。
3. 了解MATLAB有关窗函数法设计的常用子函数。

实验环境

Windows 11, MATLAB2021a

示例程序

设计FIR数字带阻滤波器

用凯塞窗设计一个长度为75的FIR数字带阻滤波器，要求：下通带截止频率为$ω_{p1}=0.2π$,$R_p=0.1dB$; 阻带低端截止频率$ω_{s1}=0.3π$,$A_s=60dB$; 阻带高端截止频率$ω_{s2}=0.7π$,$A_s=60dB$; 上通带截止频率为$ω_{p2}=0.8π$, $R_p=0.1dB$.编写程序实现该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

beta=0.1102*(60-8.7);
wp1=0.2*pi;wp2=0.8*pi;
ws1=0.3*pi;ws2=0.7*pi;
wp=[wp1,wp2];ws=[ws1,ws2];
deltaw=ws1-wp1;%过渡带带宽计算
N0=ceil(7.5*pi/deltaw);%滤波器长度
N=N0+mod(N0+1,2);%为实现FIR类型I偶对称滤波器，应确保N为奇数
windows=(kaiser(75,beta))';
wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;%将通阻带频率的平均值作为截止频率
hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N);%建立理想带阻
b=hd.*windows;
[db,mag,pha,grd,w]=freqz_m(b,1);%求解频率特性
n=0:N-1;dw=2*pi/1000;%设置频率分辨率为2*pi/1000
wp0=[1:wp1/dw+1,wp2/dw+1:501];%建立通带频率样点数组
As=-round(max(db(ws1/dw+1:ws2/dw+1)));%检验最小阻带衰减
Rp=-(min(db(wp0)));%检验带通波动

%% 绘图
subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]);
title('实际脉冲响应');xlabel('n');ylabel('h(n)');
subplot(2,2,2);stem(n,windows);
axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)');
subplot(2,2,3);plot(w/pi,db);axis([0,1,-150,10]);
title('幅度频率响应');xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})');
set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]);
set(gca,'YTickMode','manual','YTick',[-150,-40,-3,0]);grid
subplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]);
title('相位频率响应');xlabel('频率（单位：\pi）');ylabel('\phi(\omega)');
set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]);
set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid

%% 自定义函数
function hd=ideal_lp(wc,M);alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
end

function [db,mag,pha,grd,w]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';
w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
end

其中，有子程序ideal_lp.m和freqz_m.m定义对应函数：

%% 子程序‘ideal_lp.m’
function [db,mag,pha,grd,w] = freqz_m(b,a)
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
grd = grpdelay(b,a,w);

%% 子程序‘freqz_m.m’
function hd=ideal_lp(wc,M)
alpha = (M-1)/2;
n = [0:1:(M-1)];
m = n - alpha +eps;
%加一个小数以避免0 作为除数
hd = sin(wc*m)./(pi*m);