当前位置：首页 > news >正文

实验七用 MATLAB 设计 FIR 数字滤波器

news 文章来源：https://blog.csdn.net/weixin_73384981/article/details/144196777 2025/5/1 11:14:29

实验目的

加深对窗函数法设计 FIR 数字滤波器的基本原理的理解。
学习用 Matlab 语言的窗函数法编写设计 FIR 数字滤波器的程序。
了解 Matlab 语言有关窗函数法设计 FIR 数字滤波器的常用函数用法。
掌握 FIR 滤波器的快速卷积实现原理。
不同滤波器的设计方法具有不同的优缺点，因此要全面、客观看待可能面对或出现的问题。

实验原理

用窗函数法设计FIR 数字滤波器

FIR数字滤波器的系统函数为

$H(z)=\sum_{n=0}^{N-1} h(n) z^{-n}$

这个公式也可以看成是离散 LSI 系统的系统函数

$H(z)=\frac{Y(z)}{X(z)}=\frac{b(z)}{a(z)}=\frac{\sum_{m=0}^M b_m z^{-m}}{1+\sum_{k=1}^N a_k z^{-k}}=\frac{b_0+b_1 z^{-1}+b_2 z^{-2}+\cdots+b_m z^{-m}}{1+a_1 z^{-1}+a_2 z^{-2}+\cdots+a_k z^{-k}}$

分母 $a_0$ 为 1, 其余 $a_k$ 全都为 0时的一个特例。由于极点全部集中在零点, 稳定和线性相位特性是 FIR滤波器的突出优点, 因此在实际中广泛使用。

FIR 滤波器的设计任务是选择有限长度的 $h (n)$ , 使传输函数 $H\left(e^{j\omega}\right)$ 满足技术要求。主要设计方法有窗函数法、频率采样法和切比雪夫等波纹逼近法等。本实验主要介绍窗函数法。

用窗函数法设计 FIR 数字滤波器的基本步骤如下:

(1) 根据过渡带和阻带衰减指标选择窗函数的类型, 估算滤波器的阶数 $N$ 。

(2) 由数字滤波器的理想频率响应 $H\left(e^{j\omega}\right)$ 求出其单位冲激响应 $h_{d}(n)$ 。可用自定义函数 $ideal_lp \text{ideal\_lp}$ 实现理想数字低通滤波器频率响应的求解。程序清单如下:

function hd = ideal_lp(wc, N) % 点 0 到 N-1 之间的理想脉冲响应% wc = 截止频率 (弧度)% N = 理想滤波器的长度tao = (N-1)/2;n = [0:(N-1)];m = n - tao + eps; % 加一个小数以避免 0 作除数hd = sin(wc * m) ./ (pi * m);
end

其它选频滤波器可以由低通频响特性合成。如一个通带在 $\omega_{c1} \sim \omega_{c2}$ 之间的带通滤波器在给定 $N$ 值的条件下，可以用下列程序实现:

Hd = ideal_lp(wc2, N) - ideal_lp(wc1, N);

(3) 计算数字滤波器的单位冲激响应 $h (n) = w (n) h d (n)$ 。

(4) 检查设计的滤波器是否满足技术指标。

如果设计的滤波器不满足技术指标，则需要重新选择或调整窗函数的类型，估算滤波器的阶数 $N$ 。再重复前面的四个步骤，直到满足指标。

常用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、切比雪夫窗、布莱克曼窗、凯塞窗等，Matlab均有相应的函数可以调用。另外，MATLAB 信号处理工具箱还提供了 'firl’函数，可以用于窗函数法设计 FIR 滤波器。

由于第一类线性相位滤波器 (类型 I)能进行低通、高通、带通、带阻滤波器的设计，因此，本实验所有滤波器均采用第一类线性相位滤波器。

各种窗函数特性的比较

窗函数	旁瓣峰值/dB	近似过渡带宽	精确过渡带宽	阻带最小衰减/dB
矩形窗	-13	$4\pi/N$	$1.8\pi/N$	21
三角形窗	-25	$8\pi/N$	$6.1\pi/N$	25
汉宁窗	-31	$8\pi/N$	$6.2\pi/N$	44
哈明窗	-41	$8\pi/N$	$6.6\pi/N$	53
布莱克曼窗	-57	$12\pi/N$	$11\pi/N$	74
凯塞窗	-57	无	$10\pi/N$	80

信号的整数倍零值内插

设 $x (n)$ 是连续信号 $x_a(t)$ 的采样序列，其采样频率为 $f_1 = \frac{1}{T_1} (\text{Hz}), T_1$ 是采样间隔。如果将其采样频率降低到原来的 $\frac{1}{D}$ （ $D$ 为大于 1的整数，称为抽取因子），最简单的方法是对 $x (n)$ 每 $D - 1$ 个点抽取 1点，组成一个新的序列 $y (n)$ 。由于 $y (n)$ 的采样间隔 $T_2 = DT_1$ ，除非抽取后仍能满足采样定理，否则会引起频谱混叠现象。信号抽取前后的频谱关系见教材第8 章的 8.2节。为了避免抽取后的频率混叠，在抽取前先采用一个抗混叠低通滤波器对信号滤波，把信号的频带限制在某个频率以下。

抗混叠滤波器的系统函数为：

$H(e^{j\omega}) = \begin{cases} 1 & |\omega| < \frac{\pi}{D} \\ 0 & \frac{\pi}{D} \leq |\omega| \leq \pi \end{cases}$

例题

例2-2 选择合适的窗函数设计一个FIR 数字低通滤波器

要求：通带截止频率为ωp=0.3 π ，Rp=0.05dB；阻带截止频率为ωs=0.45 π ，As=50dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。分析：根据设计指标要求，并查表2-1，选择哈明窗。程序清单如下：

addpath('.\2024-2025（1）《信号处理实验》资料（请拷贝到u盘中，每次试验带到实验室）\实验中用到的一些函数和音频图像文件');
wp = 0.3*pi;
ws = 0.45*pi;
deltaw = ws - wp;
N0 = ceil(6.6*pi/deltaw);N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = (hamming(N))';
wc = (ws + wp) / 2;hd = ideal_lp(wc, N);
b = hd .* windows;[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));pha = angle(H);
n = 0:N-1;
dw = 2 * pi / 1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减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, -80, 10]);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp/pi, ws/pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-50, -20, -3, 0]);
grid on;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, wp/pi, ws/pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]);
grid on;

例2-3 用Matlab 信号处理工具箱的fir1 函数设计一个FIR 数字低通滤波器

要求：通带截止频率为ωp=0.3 π ，Rp=0.05dB；阻带截止频率为ωs=0.45 π ，As=50dB。

close;clear;wp = 0.3*pi;
ws = 0.45*pi;
deltaw = ws - wp;
N0 = ceil(6.6*pi/deltaw);
N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = hamming(N); % 使用哈明窗，此句可省略
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc, windows); % 用 fir1 函数求系统函数系数，windows 可省略[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);n = 0:N-1;
dw = 2*pi/1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减draw(n,b,N,windows,w,db,wp,ws,pha,As);

实验内容

1、阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。

2、选择合适的窗函数设计FIR 数字低通滤波器

要求：ωp=0.24 π ，Rp=0.1dB；ωs=0.3 π ，As=60dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

close;clear;
As=60;
wp = 0.24*pi;
ws = 0.3*pi;
deltaw = ws - wp;
N0 = ceil(11*pi/deltaw);
N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = blackman(N); % 使用布莱克曼窗
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc, windows); % 用 fir1 函数求系统函数系数，windows 可省略[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);n = 0:N-1;
dw = 2*pi/1000;% Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
% As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减draw(n,b,N,windows,w,db,wp,ws,pha,As);

3、调用信号产生函数xtg 产生具有加性噪声的信号x(t)，并显示信号及其频谱。

addpath('.\2024-2025（1）《信号处理实验》资料（请拷贝到u盘中，每次试验带到实验室）\实验中用到的一些函数和音频图像文件');
[xt,tf]=xtg(1000);

4、采用实验内容步骤2 中设计的FIR 数字低通滤波器

调用Matlab 快速卷积函数fftfilt 实现对x(t)的滤波，从高频噪声中提取x(t)中的单频调幅信号。绘图显示滤波器的频率响应特性曲线、滤波器输出信号的幅频特性图和时域波形图。

yt = fftfilt(b, xt);Hyk = abs(fft(yt));
figure;subplot(2, 1, 1);
plot(tf, yt);
axis([0, 1, -2, 2]);subplot(2, 1, 2);
stem(Hyk);axis([80, 120, min(Hyk), max(Hyk)]);

5、选做题，读取音频信号 motherland.wav，得到 xn

1 对 xn进行 I=2的整数倍 0值内插，得到音频信号 yn1

2 设计一个镜像低通滤波器（可在实验内容 2的代码上进行修改）

3 对 yn1进行滤波，得到音频信号 yn2。

①音频播放：依次播放原音频信号 xn、 yn1和 yn2，体验整数倍 0值内插后的音质。

close;clear;
[xn,fs]=audioread('motherland.wav');% 读取音频信号
% sound(xn,fs); % 播放音频信号，
% pause(length(xn)/fs); % 暂停执行程序 length(xn)/fs秒，确保音频播 放完。
I=2; % 实现 I=2的整数倍 0值内插
for i=1:length(xn);
yn1(I*i-1)=xn(i);
yn1(I*i)= 0;
end
% sound(yn1,I*fs); %采样频率变大了，为 I*fs

②取原音频某段信号，如 n=8000~8199 。画出该段信号模拟域幅度谱（横坐标为 f Hz ）；画出该段信号 I=2 内插后的模拟域幅度谱；画出该段信号内插后再经过镜像滤波后的模拟域幅度谱。

addpath('.\2024-2025（1）《信号处理实验》资料（请拷贝到u盘中，每次试验带到实验室）\实验中用到的一些函数和音频图像文件');
wp = 0.48*pi;
ws = 0.52*pi;
deltaw = ws - wp;
N0 = ceil(6.6*pi/deltaw);
N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = (hamming(N))';
wc = (ws + wp) / 2;hd = ideal_lp(wc, N);
b = hd .* windows;
b

b = 1x165    0.0000    0.0003   -0.0000   -0.0003    0.0000    0.0004   -0.0000   -0.0004    0.0000    0.0005   -0.0000   -0.0005    0.0000    0.0006    0.0000   -0.0007    0.0000    0.0009   -0.0000   -0.0010    0.0000    0.0012   -0.0000   -0.0013   -0.0000    0.0015   -0.0000   -0.0018    0.0000    0.0020    0.0000   -0.0023   -0.0000    0.0026   -0.0000   -0.0029    0.0000    0.0033   -0.0000   -0.0037   -0.0000    0.0042   -0.0000   -0.0047    0.0000    0.0052   -0.0000   -0.0059   -0.0000    0.0065

[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));pha = angle(H);
n = 0:N-1;
dw = 2 * pi / 1000;


N=2048;
Xn=1/fs*fft(xn(8000:8199),N); % 从 xn中取 200点做谱分析， N可取 2048
plot((0:N/2-1)*fs/N,abs(Xn(1:N/2)));% 模拟域幅度谱

Yn1=1/(I*fs)*fft(yn1(16000:16399),N); % 内插后，200点长变成了400点长
plot((0:N/2-1)*I*fs/N,abs(Yn1(1:N/2)));

yn2=filter(b,1,yn1); % 对yn1进行滤波，b为所设计的镜像滤波器
Yn2=1/(I*fs)*fft(yn2(16000:16399),N); % 内插后，200点长变成了400点长
plot((0:N/2-1)*I*fs/N,abs(Yn2(1:N/2)));

思考题

根据下面指标要求设计四种不同类型的FIR 线性相位数字滤波器，要求画出h(n)，幅频特性曲线、幅频衰减特性曲线，相频特性曲线，标注相关信息，包括横坐标，纵坐标的单位，曲线名称。

close;clear;wp = 0.2*pi;
ws = 0.3*pi;
deltaw = ws - wp;
N0 = ceil(6.1*pi/deltaw);
N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = triang(N); % 使用哈明窗，此句可省略
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc, windows); % 用 fir1 函数求系统函数系数，windows 可省略
[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);n = 0:N-1;
dw = 2*pi/1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减draw2(n,b,N,w,db,wp,ws,mag,pha,As);

close;clear;wp = 0.6*pi;
ws = 0.4*pi;
deltaw = wp - ws;
N0 = ceil(6.2*pi/deltaw);
N = N0 + mod(N0 + 1, 2); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = hanning(N); % 使用哈明窗，此句可省略
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc,'high', windows); % 用 fir1 函数求系统函数系数，windows 可省略[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);n = 0:N-1;
dw = 2*pi/1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
% As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减
As=43;
draw2(n,b,N,w,db,wp,ws,mag,pha,As);

close;clear;wp = [0.2,0.6]*pi;
ws = [0.15,0.65]*pi;
deltaw = wp - ws;
N0 = ceil(6.6*pi./deltaw);
N = max(N0 + mod(N0 + 1, 2)); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = hamming(N); % 使用哈明窗，此句可省略
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc,'bandpass'); % 用 fir1 函数求系统函数系数，windows 可省略[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);
nextfigure
n = 0:N-1;
dw = 2*pi/1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
% As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减
As=50;
draw3(n,b,N,w,db,wp,ws,mag,pha,As);

close;clear;wp = [0.2,0.6]*pi;
ws = [0.15,0.65]*pi;
deltaw = wp - ws;
N0 = ceil(6.6*pi./deltaw);
N = max(N0 + mod(N0 + 1, 2)); % 为实现 FIR 类型 1 偶对称滤波器，应确保 N 为奇数
windows = hamming(N); % 使用哈明窗，此句可省略
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N-1, wc,'stop'); % 用 fir1 函数求系统函数系数，windows 可省略[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';mag = abs(H);
db = 20*log10((mag + eps) / max(mag));
pha = angle(H);n = 0:N-1;
dw = 2*pi/1000;Rp = -(min(db(1:wp/dw + 1))); % 检验通带波动
% As = -round(max(db(ws/dw + 1:501))); % 检验最小阻带衰减
As=45;
draw3(n,b,N,w,db,wp,ws,mag,pha,As);

绘图函数

用到的绘图函数，可以以.m文件保存，放在同一路径下。或addpath加入函数的绝对/相对路径。

% function draw(n,b,N,windows,w,db,wp,ws,pha,As)
% 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, -80, 10]);
% title('幅度频率响应');
% xlabel('频率（单位：\pi）');
% ylabel('H(e^{j\omega})');
% set(gca, 'XTickMode', 'manual', 'XTick', [0, wp/pi, ws/pi, 1]);
% set(gca, 'YTickMode', 'manual', 'YTick', [-As, -20, -3, 0]);
% grid on;
% 
% 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, wp/pi, ws/pi, 1]);
% set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]);
% grid on;

% function draw2(n,b,N,w,db,wp,ws,mag,pha,As)
% 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);
% plot(w/pi, mag); 
% xlabel('频率（单位：\pi）');
% ylabel('|H|'); 
% title('幅频特性曲线'); 
% axis([0, 1, 0, 1.1]);
% set(gca, 'XTickMode', 'manual', 'XTick', [0, min([wp/pi,ws/pi]), max([wp/pi,ws/pi]), 1]);
% set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.707, 1]);
% grid on;
% 
% subplot(2, 2, 3);
% plot(w/pi, db);
% axis([0, 1, -80, 10]);
% title('幅频衰减特性曲线');
% xlabel('频率（单位：\pi）');
% ylabel('H(e^{j\omega})');
% set(gca, 'XTickMode', 'manual', 'XTick', [0, min([wp/pi,ws/pi]), max([wp/pi,ws/pi]), 1]);
% set(gca, 'YTickMode', 'manual', 'YTick', [min([-As,-20]), max([-As,-20]), -3, 0]);
% grid on;
% 
% 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, min([wp/pi,ws/pi]), max([wp/pi,ws/pi]), 1]);
% set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]);
% grid on;

% function draw3(n,b,N,w,db,wp,ws,mag,pha,As)
% unique_ticks = unique([0, wp/pi, ws/pi, 1]);
% 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);
% plot(w/pi, mag); 
% xlabel('频率（单位：\pi）');
% ylabel('|H|'); 
% title('幅频特性曲线'); 
% axis([0, 1, 0, 1.1]);
% set(gca, 'XTickMode', 'manual', 'XTick', unique_ticks);
% set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.707, 1]);
% grid on;
% 
% subplot(2, 2, 3);
% plot(w/pi, db);
% axis([0, 1, -80, 10]);
% title('幅频衰减特性曲线');
% xlabel('频率（单位：\pi）');
% ylabel('H(e^{j\omega})');
% set(gca, 'XTickMode', 'manual', 'XTick', unique_ticks);
% set(gca, 'YTickMode', 'manual', 'YTick', [min([-As,-20]), max([-As,-20]), -3, 0]);
% grid on;
% 
% subplot(2, 2, 4);
% plot(w/pi, pha);
% axis([0, 1, -4, 4]);
% title('相频特性曲线');
% xlabel('频率（单位：\pi）');
% ylabel('\phi(\omega)');
% 
% set(gca, 'XTickMode', 'manual', 'XTick', unique_ticks);
% set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]);
% grid on;