企网站的互联网,沈阳高端网站制作公司哪家好,辽宁城乡和住房建设部网站,seo的搜索排名影响因素有实验目的 加深对窗函数法设计 FIR 数字滤波器的基本原理的理解。 学习用 Matlab 语言的窗函数法编写设计 FIR 数字滤波器的程序。 了解 Matlab 语言有关窗函数法设计 FIR 数字滤波器的常用函数用法。 掌握 FIR 滤波器的快速卷积实现原理。 不同滤波器的设计方法具有不同的优…实验目的 加深对窗函数法设计 FIR 数字滤波器的基本原理的理解。 学习用 Matlab 语言的窗函数法编写设计 FIR 数字滤波器的程序。 了解 Matlab 语言有关窗函数法设计 FIR 数字滤波器的常用函数用法。 掌握 FIR 滤波器的快速卷积实现原理。 不同滤波器的设计方法具有不同的优缺点因此要全面、客观看待可能面对或出现的问题。
实验原理
用窗函数法设计FIR 数字滤波器
FIR数字滤波器的系统函数为 H ( z ) ∑ n 0 N − 1 h ( n ) z − n H(z)\sum_{n0}^{N-1} h(n) z^{-n} H(z)n0∑N−1h(n)z−n
这个公式也可以看成是离散 LSI 系统的系统函数 H ( z ) Y ( z ) X ( z ) b ( z ) a ( z ) ∑ m 0 M b m z − m 1 ∑ k 1 N a k z − k b 0 b 1 z − 1 b 2 z − 2 ⋯ b m z − m 1 a 1 z − 1 a 2 z − 2 ⋯ a k z − k H(z)\frac{Y(z)}{X(z)}\frac{b(z)}{a(z)}\frac{\sum_{m0}^M b_m z^{-m}}{1\sum_{k1}^N a_k z^{-k}}\frac{b_0b_1 z^{-1}b_2 z^{-2}\cdotsb_m z^{-m}}{1a_1 z^{-1}a_2 z^{-2}\cdotsa_k z^{-k}} H(z)X(z)Y(z)a(z)b(z)1∑k1Nakz−k∑m0Mbmz−m1a1z−1a2z−2⋯akz−kb0b1z−1b2z−2⋯bmz−m
分母 a 0 a_0 a0 为 1, 其余 a k a_k ak 全都为 0时的一个特例。由于极点全部集中在零点, 稳定和线性相位特性是 FIR滤波器的突出优点, 因此在实际中广泛使用。
FIR 滤波器的设计任务是选择有限长度的 h ( n ) h(n) h(n), 使传输函数 H ( e j ω ) H\left(e^{j\omega}\right) H(ejω)满足技术要求。主要设计方法有窗函数法、频率采样法和切比雪夫等波纹逼近法等。本实验主要介绍窗函数法。
用窗函数法设计 FIR 数字滤波器的基本步骤如下:
(1) 根据过渡带和阻带衰减指标选择窗函数的类型, 估算滤波器的阶数 N N N。
(2) 由数字滤波器的理想频率响应 H ( e j ω ) H\left(e^{j\omega}\right) H(ejω)求出其单位冲激响应 h d ( n ) h_{d}(n) hd(n)。 可用自定义函数 ideal_lp \text{ideal\_lp} 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 其它选频滤波器可以由低通频响特性合成。如一个通带在 ω c 1 ∼ ω c 2 \omega_{c1} \sim \omega_{c2} ωc1∼ωc2 之间的带通滤波器在给定 N N N值的条件下可以用下列程序实现:
Hd ideal_lp(wc2, N) - ideal_lp(wc1, N);(3) 计算数字滤波器的单位冲激响应 h ( n ) w ( n ) h d ( n ) h(n) w(n)hd(n) h(n)w(n)hd(n)。
(4) 检查设计的滤波器是否满足技术指标。
如果设计的滤波器不满足技术指标则需要重新选择或调整窗函数的类型估算滤波器的阶数 N N N。再重复前面的四个步骤直到满足指标。
常用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、切比雪夫窗、布莱克曼窗、凯塞窗等Matlab均有相应的函数可以调用。另外MATLAB 信号处理工具箱还提供了 firl’函数可以用于窗函数法设计 FIR 滤波器。
由于第一类线性相位滤波器 (类型 I)能进行低通、高通、带通、带阻滤波器的设计因此本实验所有滤波器均采用第一类线性相位滤波器。
各种窗函数特性的比较
窗函数旁瓣峰值/dB近似过渡带宽精确过渡带宽阻带最小衰减/dB矩形窗-13 4 π / N 4\pi/N 4π/N 1.8 π / N 1.8\pi/N 1.8π/N21三角形窗-25 8 π / N 8\pi/N 8π/N 6.1 π / N 6.1\pi/N 6.1π/N25汉宁窗-31 8 π / N 8\pi/N 8π/N 6.2 π / N 6.2\pi/N 6.2π/N44哈明窗-41 8 π / N 8\pi/N 8π/N 6.6 π / N 6.6\pi/N 6.6π/N53布莱克曼窗-57 12 π / N 12\pi/N 12π/N 11 π / N 11\pi/N 11π/N74凯塞窗-57无 10 π / N 10\pi/N 10π/N80
信号的整数倍零值内插
设 x ( n ) x(n) x(n) 是连续信号 x a ( t ) x_a(t) xa(t) 的采样序列其采样频率为 f 1 1 T 1 ( Hz ) , T 1 f_1 \frac{1}{T_1} (\text{Hz}), T_1 f1T11(Hz),T1是采样间隔。如果将其采样频率降低到原来的 1 D \frac{1}{D} D1 D D D 为大于 1的整数称为抽取因子最简单的方法是对 x ( n ) x(n) x(n) 每 D − 1 D-1 D−1 个点抽取 1点组成一个新的序列 y ( n ) y(n) y(n)。由于 y ( n ) y(n) y(n) 的采样间隔 T 2 D T 1 T_2 DT_1 T2DT1除非抽取后仍能满足采样定理否则会引起频谱混叠现象。信号抽取前后的频谱关系见教材第8 章的 8.2节。为了避免抽取后的频率混叠在抽取前先采用一个抗混叠低通滤波器对信号滤波把信号的频带限制在某个频率以下。
抗混叠滤波器的系统函数为 H ( e j ω ) { 1 ∣ ω ∣ π D 0 π D ≤ ∣ ω ∣ ≤ π H(e^{j\omega}) \begin{cases} 1 |\omega| \frac{\pi}{D} \\ 0 \frac{\pi}{D} \leq |\omega| \leq \pi \end{cases} H(ejω){10∣ω∣DπDπ≤∣ω∣≤π
例题
例2-2 选择合适的窗函数设计一个FIR 数字低通滤波器
要求通带截止频率为ωp0.3 π Rp0.05dB阻带截止频率为ωs0.45 π As50dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。分析根据设计指标要求并查表2-1选择哈明窗。程序清单如下
addpath(.\2024-20251《信号处理实验》资料请拷贝到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 数字低通滤波器
要求通带截止频率为ωp0.3 π Rp0.05dB阻带截止频率为ωs0.45 π As50dB。
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 数字低通滤波器 要求ωp0.24 π Rp0.1dBωs0.3 π As60dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。
close;clear;
As60;
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-20251《信号处理实验》资料请拷贝到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进行 I2的整数倍 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秒确保音频播 放完。
I2; % 实现 I2的整数倍 0值内插
for i1:length(xn);
yn1(I*i-1)xn(i);
yn1(I*i) 0;
end
% sound(yn1,I*fs); %采样频率变大了为 I*fs②取原音频某段信号如 n8000~8199 。画出该段信号模拟域幅度谱横坐标为 f Hz 画出该段信号 I2 内插后的模拟域幅度谱画出该段信号内插后再经过镜像滤波后的模拟域幅度谱。
addpath(.\2024-20251《信号处理实验》资料请拷贝到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;
bb 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;N2048;
Xn1/fs*fft(xn(8000:8199),N); % 从 xn中取 200点做谱分析 N可取 2048
plot((0:N/2-1)*fs/N,abs(Xn(1:N/2)));% 模拟域幅度谱Yn11/(I*fs)*fft(yn1(16000:16399),N); % 内插后200点长变成了400点长
plot((0:N/2-1)*I*fs/N,abs(Yn1(1:N/2)));yn2filter(b,1,yn1); % 对yn1进行滤波b为所设计的镜像滤波器
Yn21/(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))); % 检验最小阻带衰减
As43;
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))); % 检验最小阻带衰减
As50;
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))); % 检验最小阻带衰减
As45;
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;
