FIR数字滤波器的设计及软件实现

（1）认真复习第七章中用窗函数法和等波纹最佳逼近法设计 FIR 数字滤波 器的原理； （2）调用信号产生函数 xtg 产生具有加性噪声的信号 xt，并自动显示 xt 及其频谱； （3）请设计低通滤波器， 从高频噪声中提取 xt 中的单频调幅信号，要求信 号幅频失真小于 0.1dB，将噪声频谱衰减 60dB。先观察 xt 的频谱，确定滤波器 指标参数； （4）根据滤波器指标选择合适的窗函数， 计算窗函数的长度 N， 调用 MATLAB 函数 fir1 设计一个 FIR 低通滤波器。并编写程序，调用 MATLAB 快速卷积函数 fftfilt 实现对 xt 的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号 的幅频特性图和时域波形图； （5）重复（3），滤波器指标不变，但改用等波纹最佳逼近法，调用 MATLAB 函数 remezord 和 remez 设计 FIR 数字滤波器。并比较两种设计方法设计的滤波 器阶数。
0
滤 波 后 的 y 2(t)的 波 形
y 2(t) y 2(t)
0.5
-0.5
0
-0.5 -1 -1
0
0
0.05
0
0.05
0.25 0.3 0.35 0.4 t/s 0.1 0.25 0.3 0.4 等 波 纹 0.15 逼等波纹输出滤波（图 近 法0.2 低通滤 波器的 损 耗0.35 函数曲 线 5 ）
(b)信 号 加 噪 声 的 频 谱 1
幅度
0.5
0
0
50
100
150
200
250 f/Hz
300
350
400
450
500
加噪声频谱（图 2）
窗函数输出波形：
1
1 0.5
y1(t) y1(t)
滤 波 后 的 y(t)的 波 形
滤 波 后 的 y(t)的 波 形
0.5 0 0 -0.5 -0.5 -1 -1
0.1
0.15
0.2
0.45
0.5
t/s 等波纹逼近法低通滤波器的损耗函数曲线
0.45
0.5
度 （ dB） 幅 度 （幅 dB ）
-20 0 -40
-20
-60
-40
-80
-60 0 -80 0
0.1
0.2
0.3
0.4
0.5 ω/π
0.6
0.7
0.8
0.9
1
0.1
0.2
Hale Waihona Puke 0.30.40.5 ω/π
0.6
x(t)
0 -5 -10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
(b)信 号 加 噪 声 的 频 谱 滤波输出波形（图 1） 1
幅度
2
0.5
x(
0 -5 -10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
FIR 数字滤波器的设计及软件实现
一、实验目的
（1）掌握用窗函数法设计 FIR 数字滤波器的原理和方法； （2）掌握用等波纹最佳逼近法设计 FIR 数字滤波器的原理和方法； （3）掌握 FIR 滤波器的快速卷积实现原理； （4）学会调用 MATLAB 函数设计与实现 FIR 滤波器。
二、实验内容及步骤




五、实验结果
1．采样指标：采样频率 Fs=1000Hz，采样周期 T=1/Fs； 2．滤波器指标参数：通带截止频率 fp=120Hz，阻带截至频率 fs=150Hz，换 算成数字频率，通带截止频率，通带最大衰为 0.1dB，阻带截至频率 ，阻带最 小衰为 60dB。 原信号的输出：
(a)信 号 加 噪 声 波 形 10 5
0.45
0.45
0.5
0.5
0
-50
-50
-100
-100 0 0
0.1
0.1
0.2
0.2
0.3
0.3
0.4
0.4
0.5 ω/π 0.5 ω/π
0.6
0.6
0.7
0.7
0.8
0.8
0.9
0.9
1
1
窗函数输出损耗（图 4）
等波纹最佳逼近法输出波形：
3
1 0.5
滤 波 后 的 y 2(t)的 波 形
1
三、实验程序框图
Fs=1000，T=1/Fs
xt=xtg 产生信号 xt，并显示 xt 及其频谱
用窗函数法或等波纹最佳逼近法设计 FIR 滤波器 hn
对信号 xt 滤波：yt=fftfilt(hn,xt)
1.计算并绘图显示滤波器损耗函数 2.绘图显示滤波器输出信号 yt 及其频谱
结束
1
四、思考题
1) 用窗函数法设计线性相位低通滤波器的设计步骤教材中有详细的介绍. 答：（1）根据对阻带衰减及过度带的指标要求，选择窗函数的类型 w(n)和 估计窗口长度 N； （2）构造希望逼近的频率响应函数 Hd(ejw); （3）计算 hd(n)； （4）加窗好 h(n)=hd(n)w(n)； （5）求实际滤波器的频响； （6）判断是否达到要求，若未达到要求，另选窗函数，改变窗长度，重复 (3)(4)(6)； （7）若达到要求由 h(n)画出实现网络结构图； 2) 希望逼近的理想带通滤波器的截止频率 分别为： 答： wcl wsl w pl ， wcu wsu w pu 2 2 3) 解释为什么对同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低？ ①用窗函数法设计的滤波器，如果在阻带截止频率附近刚好满足，则离开阻 带截止频率越远，阻带衰减富裕量越大，即存在资源浪费； ② 几种常用的典型窗函数的通带最大衰减和阻带最小衰减固定，且差别较 大， 又不能分别控制。所以设计的滤波器的通带最大衰减和阻带最小衰减通常都 存在较大富裕。如本实验所选的 blackman 窗函数，其阻带最小衰减为 74dB,而 指标仅为 60dB； ③ 用等波纹最佳逼近法设计的滤波器，其通带和阻带均为等波纹特性，且 通带最大衰减和阻带最小衰减可以分别控制，所以其指标均匀分布，没有资源浪 费，所以期阶数低得多
4
附录：
% FIR 数字滤波器设计及软件实现 clear all; %==调用 xtg 产生信号 xt, xt 长度 N=1000,并显示 xt 及其频谱,========= N=1000;xt=xtg;Fs=1000; fp=120; fs=150;Rp=0.1;As=60; % 输入给定指标 % (1) 用窗函数法设计滤波器 wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于 pi 归一化） B=2*pi*(fs-fp)/Fs; %过渡带宽度指标 Nb=ceil(11*pi/B); %blackman 窗的长度 N hn=fir1(Nb-1,wc,blackman(Nb)); Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性 ywt=fftfilt(hn,xt,N); %调用函数 fftfilt 对 xt 滤波 %以下为用窗函数法设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形） figure(2); t=0:0.001:0.999;%绘制滤波后的信号时域波形图 subplot(2,1,1); plot(t,ywt); grid; xlabel('t/s');ylabel('y1(t)'); title('滤波后的 y(t)的波形'); axis([0 0.5 -1 1]); subplot(2,1,2); [h w]=freqz(hn);%绘制低通滤波器的损耗函数曲线 plot(w/pi,20*log10(abs(h))); grid; xlabel('ω/π');ylabel('幅度（dB）'); title('窗函数法低通滤波器的损耗函数曲线'); axis([0 1 -120 5]); %(2)用等波纹最佳逼近法设计滤波器 fb=[fp,fs];m=[1,0];%确定 remezord 函数所需参数 f,m,dev dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; [Ne,fo,mo,W]=remezord(fb,m,dev,Fs);%确定 remez 函数所需参数 hn=remez(Ne,fo,mo,W);%调用 remez 函数进行设计 y2t=fftfilt(hn,xt,N);%调用函数 fftfilt 对 xt 滤波 figure(3); t=0:0.001:0.999;%绘制滤波后的信号时域波形图 subplot(2,1,1);plot(t,y2t); grid; xlabel('t/s');ylabel('y_2(t)'); title('滤波后的 y_2(t)的波形'); axis([0 0.5 -1 1]);
6
5
subplot(2,1,2); [h w]=freqz(hn);%绘制低通滤波器的损耗函数曲线 plot(w/pi,20*log10(abs(h))); grid; xlabel('ω/π');ylabel('幅度（dB）'); title('等波纹逼近法低通滤波器的损耗函数曲线'); axis([0 1 -80 5]); %子程序 %信号产生程序 function xt=xtg N=1000;Fs=1000;T=1/Fs;Tp=N*T; t=0:T:(N-1)*T; fc=Fs/10;f0=fc/10; mt=cos(2*pi*f0*t); ct=cos(2*pi*fc*t); xt=mt.*ct; nt=2*rand(1,N)-1; fp=150;fs=200;Rp=0.1;As=70; fb=[fp,fs];m=[0,1]; dev=[10^(-As/20),(10^Rp/20-1)/10^(Rp/20)+1]; [n,fo,mo,W]=remezord(fb,m,dev,Fs); hn=remez(n,fo,mo,W); yt=filter(hn,1,10*nt); xt=xt+yt; fst=fft(xt,N); k=0:N-1;f=k/Tp; figure(1) subplot(2,1,1);plot(t,xt); grid; title('t/s');ylabel('x(t)'); axis([0,Tp/5,min(xt),max(xt)]); title('(a)信号加噪声波形') subplot(2,1,2);plot(f,abs(fst)/max(abs(fst))); grid;title('(b)信号加噪声的频谱') axis([0,Fs/2,0,1.2]); xlabel('f/Hz');ylabel('幅度')