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

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

一、实验目的

1. 掌握用窗函数法设计FIR数字滤波器的原理和方法
2. 掌握用等波纹最佳逼近法设计FIR数字滤波器的原理和方法
3. 掌握滤波器的快速卷积实现原理
4. 学会调用MATLAB函数设计与实现FIR滤波器

二、实验原理

设计FIR数字滤波器一般采用窗函数法和等波纹最佳逼近法，其原理分别如下：

1. 窗函数法设计FIR数字滤波器的原理：

如果所希望的滤波器的理想频率响应函数为，则其对应的单位冲激响应为：

 用窗函w(n)将截断，得到：，就作为实际设计的数字滤波器的单位冲激响应序列，其频率响应函数为：。如果要求线性相位特性，则还必须满足：。

根据上式中的正负号和长度正在上传…重新上传取消取为奇数或者偶数又将线性相位滤波器分成四类。

总体的思路就是：根据正在上传…重新上传取消选择合适的窗函数，根据过渡带选择窗函数的窗口长度。

1. 等波纹最佳逼近法设计正在上传…重新上传取消数字滤波器的原理：

用正在上传…重新上传取消表示希望逼近的幅度特性函数，要求设计线性相位数字滤波器时，必须满足线性相位约束条件。用表示实际设计的滤波器幅度特性函数。定义加权误差函数为：，式中，称为误差加权函数，用来控制不同频段（一般指通带和阻带）的逼近精度。等波纹最佳逼近基于切比雪夫逼近，在通带和阻带以的最大值最小化为准则，采用多重交换迭代算法求解滤波器系数。所以取值越大的频段，逼近精度越高，开始设计时应根据逼近精度要求确定，在多重交换迭代过程中是确知函数。

用等波纹最佳逼近法设计正在上传…重新上传取消数字滤波器的过程是：

1. 根据给定的逼近指标估算滤波器阶数正在上传…重新上传取消和误差加权函数；
2. 采用正在上传…重新上传取消算法得到滤波器单位脉冲响应。

三、实验内容

1. 认真复习第正在上传…重新上传取消章中用窗函数法和等波纹最佳逼近法设计数字滤波器的原理；
2. 调用信号产生函数正在上传…重新上传取消产生具有加性噪声的信号，并自动显示及其频谱，如书上图所示。
3. 请设计低通滤波器，从高频噪声中提取正在上传…重新上传取消中的单频抑制载波调幅信号，要求信号幅频失真小于0.1，将噪声频谱衰减。观察的频谱，确定滤波器指标参数。
4. 根据滤波器指标选择合适的窗函数，计算窗函数的长度正在上传…重新上传取消，调用函数设计一个低通滤波器。并编写程序，调用快速卷积函数实现对的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性曲图和时域波形图。
5. 重复正在上传…重新上传取消，滤波器指标不变，但改变等波纹最佳逼近法设计滤波器，调用函数和设计数字滤波器。比较两种设计方法设计的滤波器阶数。

提示：

①正在上传…重新上传取消函数和的功能及其调用格式请查阅本书节和节；

②采样频率正在上传…重新上传取消，采样周期；

③根据图正在上传…重新上传取消和实验要求，可选择滤波器指标参数：通带截止频率，阻带截止频率，换算成数字频率，通带截止频率，通带最大衰减为，阻带截止频率，阻带最小衰减为。

四、程序源码与运行结果

1. 两种方法的原理分析见第二部分（实验原理）。
2. 步骤分析：根据实验数提供的代码正在上传…重新上传取消产生具有加性噪声的信号，并绘制出其时域波形和频谱，下面详细阐述结果：

正在上传…重新上传取消的波形和频谱如下：

正在上传…重新上传取消

图4.1 正在上传…重新上传取消的波形和频谱图

分析：正在上传…重新上传取消信号是由频率分别为和的载波正弦信号和单频正弦波调制信号乘积得到的信号加上高频噪声得到的，可以看到，此时的时域波形已经很混乱了，很难总重提取中的有用信号部分，但是其频谱中可以看出，由于加的是高频噪声，因此频谱中噪声频率都比较高，采用低通滤波器就可以比较准确地将中的有用信号滤出来。

同时，分析函数正在上传…重新上传取消可以发现，其设置的采样频率,载波信号和调制信号频率分别为：，,高频噪声在频率小于时几乎为零，得到的中的单频抑制载波调幅信号的频谱中所含的频率分量分别为和。

程序代码：见第七部分（附录）

1. 步骤分析：首先分析正在上传…重新上传取消信号的频谱，其有用信号（单频抑制载波调幅信号）频谱分布主要在，高频噪声的频谱主要分布在大于的范围，因此根据实验要求，选择滤波器的指标参数为：通带截止频率，阻带截止频率，换算成数字频率，通带截止频率，通带最大衰减为，阻带截止频率，阻带最小衰减为。
2. 步骤分析：根据步骤正在上传…重新上传取消里面设计的滤波器，由于阻带最小衰减为，查表得，此时可以选用布莱克曼窗函数，且知道布莱克曼函数的过渡带宽度的精确值为，因此需要满足过渡带，推得，取窗函数长度（此程序中是利用函数求出的布莱克曼窗函数的长度）。再调用函数设计一个符合上述标准的低通滤波器，调用快速卷积函数实现对的滤波，下面详细阐述结果：

1. 窗函数法设计的低通滤波器的波形和幅频特性曲线如下：

正在上传…重新上传取消

图4.2 窗函数法设计的低通滤波器的波形和幅频特性曲线图

1. 窗函数法设计的滤波器输出信号的时域波形和幅频特性曲线如下：

正在上传…重新上传取消

图4.3 窗函数法设计的滤波器输出信号的时域波形和幅频特性曲线图

分析：由滤波之后得到的波形图可以看出，该滤波器对信号正在上传…重新上传取消进行滤波之后，是能够比较准确地分离出信号中的单频抑制载波调幅信号，其分离效果还可以，同时可以得到，此时的滤波器长度为，与前面理论计算得到的一样，说明仿真符合理论分析。

程序代码【程序清单】：（其中调用的正在上传…重新上传取消函数见第七部分（附录））

1. %窗函数法设计滤波器  
2. clear;  
3. N=1000;%窗函数长度  
4. xt=xtg;%产生长度为N的x(t)信号  
5. %====滤波器指标=====  
6. fp=120;  
7. fs=150;  
8. Rp=0.1;  
9. As=60;  
10. Fs=1000;%采样频率  
11. %====设计滤波器=====  
12. wc=(fp+fs)/Fs;%理想低通滤波器截止频率（关于pi归一化）  
13. B=2*pi*(fs-fp)/Fs;%过渡带宽指标  
14. Nb=ceil(11*pi/B);%布莱克窗函数的长度N  
15. hn=fir1(Nb-1,wc,blackman(Nb));  
16. Hw=abs(fft(hn,1024));%求设计的滤波器频率特性  
17. ywt=fftfilt(hn,xt,N);%调用函数fftfilt对xt滤波  
18. %====绘图=====  
19. figure;  
20. subplot(2,1,1);  
21. n=0:length(hn)-1;  
22. stem(n,hn,'.');  
23. xlabel('n');  
24. ylabel('h(n)');  
25. title('窗函数法设计的低通滤波器h(n)波形');  
26. axis([0,190,-0.1,0.3]);  
27. subplot(2,1,2);  
28. HdB=20*log10(Hw);  
29. f=0:length(HdB)-1;  
30. plot(f,HdB);  
31. xlabel('f/Hz');  
32. ylabel('幅度/dB');  
33. title('窗函数法设计的低通滤波器的幅频特性曲线');  
34. axis([0,512,-120,20]);  
35. figure;  
36. subplot(2,1,1);  
37. t=0:length(ywt)-1;  
38. t=t/Fs;  
39. plot(t,ywt);  
40. xlabel('t/s');  
41. ylabel('y_w(t)');  
42. title('滤波器输出信号时域波形曲线');  
43. axis([0,0.5,-1,1]);  
44. subplot(2,1,2);
45. ywt=fft(ywt,1024);  
46. YwdB=20*log10(abs(ywt));  
47. f=0:length(YwdB)-1;  
48. plot(f,YwdB);  
49. xlabel('f/Hz');  
50. ylabel('幅度/dB');  
51. title('滤波器输出信号的幅频特性曲线');  
52. axis([0,512,-130,20]);  

1. 步骤分析：根据步骤4正在上传…重新上传取消里面设计的滤波器指标，采用等波纹最佳逼近法设计滤波器，调用函数和设计符合上述标准的低通滤波器，调用快速卷积函数实现对的滤波，下面详细阐述结果：

1. 等波纹最佳逼近法设计的低通滤波器的波形和幅频特性曲线如下：

正在上传…重新上传取消

图4.4 等波纹最佳逼近法设计的低通滤波器的波形和幅频特性曲线图

1. 等波纹最佳逼近法设计的滤波器输出信号的时域波形和幅频特性曲线如下：

正在上传…重新上传取消

图4.5 等波纹最佳逼近法设计的滤波器输出信号的时域波形和幅频特性曲线图

分析：由滤波之后得到的波形图可以看出，该滤波器对信号正在上传…重新上传取消进行滤波之后，是能够比较准确地分离出信号中的单频抑制载波调幅信号，其分离效果还可以，同时可以得到，此时的滤波器长度为。

程序代码【程序清单】：（其中调用的正在上传…重新上传取消函数见第七部分（附录））

1. %等波纹最佳逼近法设计滤波器  
2. clear;  
3. N=1000;%窗函数长度  
4. xt=xtg;%产生长度为N的x(t)信号  
5. %====滤波器指标=====  
6. fp=120;  
7. fs=150;  
8. Rp=0.1;  
9. As=60;  
10. Fs=1000;%采样频率  
11. %====设计滤波器=====  
12. %----确定remezord函数所需参数f,m,dev-----  
13. fb=[fp,fs];  
14. m=[1,0];  
15. dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];  
16. [Nc,fo,mo,W]=remezord(fb,m,dev,Fs);%确定remez函数所需参数  
17. hn=remez(Nc,fo,mo,W);%调用remez函数进行设计  
18. Hw=abs(fft(hn,1024));%求设计的滤波器频率特性  
19. yet=fftfilt(hn,xt,N);%调用函数fftfilt对xt滤波  
20. %====绘图=====  
21. figure;  
22. subplot(2,1,1);  
23. n=0:length(hn)-1;  
24. stem(n,hn,'.');  
25. xlabel('n');  
26. ylabel('h(n)');  
27. title('等波纹最佳逼近法设计的低通滤波器h(n)波形');  
28. axis([0,100,-0.1,0.3]);  
29. subplot(2,1,2);  
30. HdB=20*log10(Hw);  
31. f=0:length(HdB)-1;  
32. plot(f,HdB);  
33. xlabel('f/Hz');  
34. ylabel('幅度/dB');  
35. title('等波纹最佳逼近法设计的低通滤波器的幅频特性曲线');  
36. axis([0,512,-80,20]);  
37. figure;  
38. subplot(2,1,1);  
39. t=0:length(yet)-1;  
40. t=t/Fs;  
41. plot(t,yet);  
42. xlabel('t/s');  
43. ylabel('y_w(t)');  
44. title('滤波器输出信号时域波形曲线');  
45. axis([0,0.5,-1,1]);  
46. subplot(2,1,2);  
47. ywt=fft(ywt,1024); 
48. YwdB=20*log10(abs(yet));  
49. f=0:length(YwdB)-1;  
50. plot(f,YwdB);  
51. xlabel('f/Hz');  
52. ylabel('幅度/dB');  
53. title('滤波器输出信号的幅频特性曲线');  
54. axis([0,512,-130,20]);  

分析比较两种方法：两种方法设计的滤波器都能够比较有效地从含有噪声的信号中提取出有用信号，但是当标准相同时，等波纹逼近法设计的滤波器的阶数明显要比窗函数设计的滤波器的阶数低，这样等波纹逼近法设计的滤波器实现的运算量以及时延都会相对小一些，可以提高设计的滤波器的性价比，窗函数设计法的优点是：总是能使通带和阻带波纹幅度相等，原理比较简单。

五、实验主程序框图

正在上传…重新上传取消

图5.1 实验主程序框图

六、实验总结

1．结论

本次实验主要是练习了通过分析信号以及给出的指标参数，分别用窗函数法和等波纹最佳逼近法设计正在上传…重新上传取消滤波器，回顾了两种方法滤波器的原理和方法，同时掌握了滤波器快速卷积的实现原理，学会调用函数设计与实现滤波器。

通过本次实验，我对于用窗函数法和等波纹最佳逼近法设计正在上传…重新上传取消滤波器的原理和方法有了更深的理解、掌握，同时学会了如何利用设计滤波器，同时对于两种设计方法有了更深入、更直观的理解，分析了两种方法的优缺点，也重温了两种方法的设计步骤。

从实验结果可以看出，相比较于窗函数设计法（和频率采样法），等波纹最佳逼近法设计的正在上传…重新上传取消滤波器的阶数更低，性价比更好，能更好地减少资源浪费。

2.思考题

1. 问题：如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何利用窗函数法设计线性相位低通滤波器？请写出设计步骤。

答：步骤如下：

1. 根据阻带最小衰减选择合适的窗函数类型，根据过渡带指标（过渡带宽度通过通带截止频率和阻带截止频率求得）估计窗口长度正在上传…重新上传取消。
2. 构造希望逼近的频率响应函数正在上传…重新上传取消，即。
3. 计算正在上传…重新上传取消。如果给出的待求滤波器的频响函数为，那么单位脉冲响应用下式求出：。
4. 加窗得到设计结果：正在上传…重新上传取消。

1. 问题：如果要求用窗函数设计带通滤波器，且给定通带上、下截止频率为正在上传…重新上传取消和，阻带上、下截止频率为和，试求理想带通滤波器的截止频率和

答：正在上传…重新上传取消，

1. 问题：解释为什么对同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低。

答：主要有以下几个原因：

1. 用窗函数法设计的滤波器，如果在阻带截止频率附近刚好满足，则离开阻带截止频率越远，阻带衰减富裕量越大，即存在资源浪费；
2. 几种常用的典型窗函数的通带最大衰减和阻带最小衰减固定，且差别较大，又不能分别控制。所以设计的滤波器的通带最大衰减和阻带最小衰减通常都存在较大富裕。如本实验所选的布莱克曼窗函数，其阻带最小衰减为正在上传…重新上传取消,而指标仅为。
3. 用等波纹最佳逼近法设计的滤波器，其通带和阻带均为等波纹特性，且通带最大衰减和阻带最小衰减可以分别控制，所以其指标均匀分布，没有资源浪费，所以期阶数低得多。

七、附录

1. 信号产生函数正在上传…重新上传取消清单（代码）
2. function xt=xtg  
3. %实验五信号x(t)产生函数，并显示信号的时域波形和幅频特性曲线  
4. %xt=xtg产生一个长度为N，有加性高斯噪声的单频调幅信号xt，N=1000  
5. %采样频率Fs=1000Hz  
6. %载波频率fc=Fs/10=100Hz，调制正弦波频率f0=fc/10=10Hz  
7. N=1000;%信号长度  
8. Fs=1000;%采样频率  
9. T=1/Fs;%采样周期  
10. Tp=N*T;%采样总时长  
11. t=0:T:(N-1)*T;  
12. fc=Fs/10;%载波频率fc=Fs/10  
13. f0=fc/10;%单频调制信号频率为f0=fc/10  
14. mt=cos(2*pi*f0*t);%产生单频正弦波调制信号mt，频率为f0  
15. ct=cos(2*pi*fc*t);%产生载波正弦波信号ct，频率为fc  
16. xt=mt.*ct;%相乘，产生单频调幅信号xt  
17. nt=2*rand(1,N)-1;%产生随机噪声nt  
18. %====设计高通滤波器hn，用于滤除噪声nt中的低频成分，生成高通噪声=====  
19. %----滤波指标-----  
20. fp=150;  
21. fs=200;  
22. Rp=0.1;  
23. As=70;  
24. %----计算remezord函数所需参数f,m,dev-----  
25. fb=[fp,fs];  
26. m=[0,1];  
27. dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];  
28. [n,fo,mo,W]=remezord(fb,m,dev,Fs);%确定remez函数所需参数  
29. hn=remez(n,fo,mo,W);%调用remez函数进行设计，用于滤除噪声nt中的低频成分  
30. yt=filter(hn,1,10*nt);%滤除随机噪声中低频成分，生成高通噪声yt  
31. %====以下为绘图部分=====  
32. xt=xt+yt;%噪声加信号  
33. fst=fft(xt,N);  
34. k=0:N-1;  
35. f=k/Tp;  
36. subplot(3,1,1);  
37. plot(t,xt);  
38. grid;  
39. xlabel('t/s');  
40. ylabel('x(t)');  
41. axis([0,Tp/5,min(xt),max(xt)]);  
42. title('(a)信号加噪声波形');  
43. subplot(3,1,2);  
44. plot(f,abs(fst)/max(abs(fst)));%归一化  
45. grid;  
46. title('(b)信号加噪声频谱');  
47. axis([0,Fs/2,0,1.2]);  
48. xlabel('f/Hz');  
49. ylabel('幅度');