使用频率采样法设计的频率采样型FIR滤波器课程设计

1、长沙理工大学数字信号处理课程设计报告 学 院 计算机与通信工程 专 业 通信工程 班 级 学 号 学生姓名 指导教师 课程成绩 完成日期 课程设计成绩评定学 院 计算机与通信工程 专 业 通信工程 班 级 学 号 学生姓名 指导教师 完成日期 2015 年 3月 13 日 指导教师对学生在课程设计中的评价评分项目优良中及格不及格课程设计中的创造性成果学生掌握课程内容的程度课程设计完成情况课程设计动手能力文字表达学习态度规范要求课程设计论文的质量指导教师对课程设计的评定意见综合成绩 指导教师签字 2015 年 3 月 15 日课程设计任务书计算机与通信工程 系 通信工程 专业 课程名称数字信号处

2、理课程设计时间20142015学年第二学期12周学生姓名 指导老师 题 目语音信号滤波去噪使用频率采样法设计的频率采样型FIR滤波器主要内容： 用麦克风采集一段语音信号，绘制波形并观察其频谱特点，加入一个带外单频噪声，使用频率采样法设计一个满足指标的FIR滤波器，对该含噪语音信号进行滤波去噪处理，画出频率采样型滤波器结构图。比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。与不同滤波方法的同学比较各种滤波方法性能的优劣。要求：（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课

3、程设计并按要求编写课程设计报告书。应当提交的文件：（1）课程设计学年论文。（2）课程设计附件（主要是源程序）。语音信号滤波去噪使用频率采样法设计的频率采样型FIR滤波器学生姓名： 指导老师： 摘 要 此次课程设计是使用频率采样法设计一个频率采样型滤波器，对一段加噪的语音信号进行滤波去噪处理，并且对滤波前后的时域波形和频谱分析滤波性能。此次课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先利用Windows录音机录制一段语音信号，并人为加入单频噪声，然后对信号进行频谱分析以确定所加噪声频率，再设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析，并用VISIO绘制滤波

4、器的结构图。由分析结果可知，滤波后的语音信号与原始语音信号基本一致，即设计的FIR 滤波器能够去除信号中所加单频噪声，即达到了设计目的。关键词 MATLAB；滤波去噪；频率采样法；频率采样型滤波器；VISIO1 引言滤波去噪1是信号处理中一种最基本但十分重要的技术。利用滤波可以从复杂的信号中提取所需的信号，抑制不需要的信号。滤波器就是这样一种可以在时域和频域对信号进行滤波处理的系统。通常情况下，有用信号和干扰信号是在不同频段上的，于是通过对滤波器的频率特性精心设计就能达到滤波的目的。本课程设计是采用频率采样法设计频率抽样型滤波器，从而对语音信号滤波去噪。通过对比滤波前后的波形图及回放滤波前后的

5、语音信号，来判断滤波器对噪声信号确实有滤除作用。1.1 课程设计目的(1)熟悉使用MATLAB；(2)了解FIR 滤波器原理及结构；(3)利用所学数字信号处理相关知识用MATLAB 设计一个FIR滤波器；(4)提高自己动手的能力。(5)对加噪的语音信号进行滤波去噪处理，比较滤波前后的时域波形和频谱并进行分析；1.2 课程设计的要求(1)滤波器指标必须符合工程设计。(2)设计完后应检查其频率响应曲线是否满足指标。(3)处理结果和分析结论应该一致，而且应符合理论。 (4)独立完成课程设计并按要求编写课程设计报告。1.3 设计平台本课程设计仿真平台为MATLAB7.0。MATLAB 的名称源自Mat

6、rix Laboratory，1984年由美国Mathworks 公司推向市场。它是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB 将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。1993 年MathWorks 公司从加拿大滑铁卢大学购得MAPLE 软件的使用权，从而以MAPLE 为“引擎”开发了符号数学工具箱（Symbolic Math Toolbox） 2。2 设计原理用windows自带的录音机录制一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用频率采样法设计的一个满足指标的频率采

7、样型FIR 滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 FIR 滤波器的设计FIR(Finite Impulse Response)滤波器3：有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。有限长单位冲激响应（FIR）滤波器有以下特点：(1) 系统的单位冲激响应h (n)在有限个n值处不为零(2) 系统函数H(z)在|z|>0处收敛，极点全部在

8、z = 0处（因果系统）(3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。2.2频率采样型结构把一个有限长序列（长度为N点）的z变换H (z)在单位圆上作N等分抽样，就得到H (k)，其主值序列就等于h (n)的离散傅里叶变换H (k)。那里也说到用H (k)表示的H (z)的内插公式为 (2.1)这个公式就为FIR滤波器提供了另外一种结构，这种结构由两部分级联组成。 (2.2)其中级联的第一部分为梳状滤波器，其结构如下图所示： (2.3)图2-1 梳状滤波器结构图第二部分由N个谐振器组成的谐振柜。它是由N个一阶网络并联组成，而这每一

9、个一阶网络都是一个谐振器 (2.4) 其结构如下图所示：图2-2 一阶谐振器 频率抽样型结构特点：(1)它的系数H(k)直接就是滤波器在处的频率响应。因此，控制滤波器的频率响应是很直接的。(2)结构有两个主要缺点：(a)所有的相乘系数及H(k)都是复数，应将它们先化成二阶的实数，这样乘起来较复杂，增加乘法次数，存储量。(b)所有谐振器的极点都是在单位圆上,由决定考虑到系数量化的影响，当系数量化时，极点会移动，有些极点就不能被梳状滤波器的零点所抵消。（零点由延时单元决定，不受量化的影响）系统就不稳定了。 (3)将一阶网络合并为二阶网络(a)第k和第N-k个谐振器合并为一个实系数的二阶网络，因为h

10、(n)是实数，它的DFT也是圆周共轭对称的。 (2.5)因此，可以将第k和第N-k个谐振器合并为一个二阶网络。 (2.6)(b)第k和第N-k个谐振器合并为一个二阶网络的极点在单位圆内，而不是在单位圆上，因而从频率响应的几何解释可知，它相当于一个有限Q的谐振器。其谐振频率为：图2-3 二阶网络结构图 除了共轭复根外，还有实根。当N=偶数时，有一对实根，它们分别为两点。 和 (2.7)当N=奇数时，只有一个实根z=r(k=0），即只有H0(z).(c)修正频率抽样结构流图（N=偶数）图2-4修正频率抽样结构流图（N=偶数） (2.8)修正频率抽样结构流图（N=奇数）图2-5修正频率抽样结构流图（

11、N=奇数） (2.9)2.3频率采样法设计思路:这种设计方法是从频域进行设计的一种方法，首先给定一个希望逼近的频率响应。 (2.10) 知道H(k)后, 由IDFT定义, 可以用这N个采样值H(k)来惟一确定有限长序列h(n)，即 (2.11) (2.12) (2.13)内插公式: (2.14)式中, ()是内插函数 (2.15)四种线性相位的FIR滤波器如下表2-1所示。表2-1 四种线性相位的FIR滤波器逼近误差及其改进措施这样设计所得到的频响H(ej)与要求的理想频响Hd(ej)会有一定的差别(1)采样点上滤波器的实际频率响应是严格地和理想频率响应数值相等的。(2)但是在采样点之间的频响

12、则是由各采样点的加权内插函数的延伸叠加而成的, 因而有一定的逼近误差， 误差大小取决于理想频率响应曲线形状。(3)理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小。(4)如果采样点之间的理想频率特性变化越陡，则内插值与理想值的误差就越大，因而在理想频率特性的不连续点附近，就会产生肩峰和起伏。滤波器性能的改善:(1)增加过渡带采样点,它可以大大减少振荡，阻带衰减也可以得到进一步改善。一般一点到二点的过渡带采样即可得到满意的结果。(2)增加采样点密度，过渡带的宽度与采样点数N成反比。但N值意味着或长度的增加，滤波器运算量必然增大4。3设计步骤3.1 设计流程图图3-1 流程图3.2

13、录制语音信号使用电脑录音机录制一段语音信号，属性设置为8000Hz，8 位单声道的音频格式。如图3.1-2 所示：图3-2 语音信号设置然后在MATLAB 软件平台下，利用函数wavread 对语音信号进行采样，源程序为：x,fs,bits=wavread('ht.wav')，记住采样频率和采样点数，matlab 实现得：fs=8000；bits=8。3.3 语音信号的频谱分析在MATLAB中编写加噪之后时域，频域对比图的源程序如下：x,fs,bits=wavread('D:ht.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，

14、fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2500; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x' y=x+0.1*sin(fn*2*pi*t); %加噪sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2

15、); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围 运行结果如下图所示：图3-3 加噪前后时域频域对比图由图3.3-1可以看出，在频域为2500Hz处加入了一个单频噪声，而加入噪声之后，时域的波形出现了明显失真，通过听取原声音信号x和加噪声音信号y，可以明显听到y语音信号中有一明显尖锐噪声。3.4 滤波器设计设计一个带阻滤波器，将单频信号滤出去，源程序如下所示：fpd=2400;fsd=2450;fsu=2550;fpu=2600;Rp=0.5;As=45; % 带阻滤波器设计指标fcd=(fpd+fsd)

16、/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率fcd fcu，和频率间隔dfwp=fcd/fs*2*pi;ws=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率M=41;alpha=(M-1)/2;%设置抽样点数l=0:M-1;wl=2*pi*l/M;T1=0.6025,T2=0.127%过渡带的抽样值NP1=fix(M*fpd/8000);NS1=fix(M*fsd/8000);NS2=ceil(M*fsu/8000);NP2=ceil(M*fpu/8000);%

17、各边界频率的抽样值1Hrs=ones(1,NP1-2),0.6025,0.127,zeros(1,NP2-NP1+1),0.127,0.6025,ones(1,M-2*NP2-5),0.6025,0.127,zeros(1,NP2-NP1+1),0.127,0.6025,ones(1,NP1-3); Hdr=1,1,0,0,1,1;%带阻滤波器的参数设置wdl1=fpd*2/fs;wdl2=fsd*2/fs;wdl3=fsu*2/fs;wdl4=fpu*2/fs;%radwdl=0,wdl1,wdl2,wdl3,wdl4,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2

18、)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,grd,w=freqz_m(h,1);Hr,ww,a,L=hr_type2(h);程序运行结果如下图3.3-2所示：图3-4 滤波器参数图由图3-4可以看出，滤波器的衰减大于设定值As=45,满足性能指标，滤波器的衰减可由增加过度带宽来得到。上图中的抽样值M41，因此阻带的带宽比较宽。可以通过调大M值来使得阻带的带宽变窄，但这样会使得滤波器的阶数相应的增加，使得其结构图标的复杂，为了是滤波器能达

19、到相应的性能要求，又使得其结构变得简单，最终取M值为41。3.5 信号滤波处理原程序如下所示：%进行滤波处理z=filter(h,1,y); % 用设计好的滤波器对y进行滤波Z=fft(z);Z=Z(1:length(Z)/2); % 计算频谱取前一半程序运行结果如下图3.3-3所示：图3-5 滤波前后时域，频域对比图 频谱取前一半由上图3-5可以看出，在滤波之后时域图能得到恢复，频域图中的单频噪声信号也的到滤除，说明了设计的滤波器能滤除加入的噪声信号，因此说达到了设计的要求。3.6 结果分析开始通过分析决定设计一个带阻的滤波器来滤除加入的单频噪声，根据噪声的频率来设计阻带的范围。在采集语音信

20、号后，按照步骤用频率采样法设计频率采样型滤波器。由图3-4可知，设计的滤波器达到要求。我们观察到图3-5 滤波前后语音信号的波形对比图，发现时域波形中加干扰噪声后有明显的变化，不过经过滤波后几乎没有变化，说明设计的滤波器达到要求。再通过听取原始语音信号，加噪信号y，和滤波之后的信号z。对比之后，发现滤波器确实滤除了噪声。从理想的角度考虑，该带阻滤波器的阻带带宽应该可以变得更窄，让滤波效果更好，但是这样的采样值会变得非常大，考虑到实际的情况，通过权衡之后，决定牺牲带宽来使得滤波器的阶数降低，因此在上图3-5中我们可以看出，在噪声频谱左右两边的信号也被滤除了。3.6 滤波器结构设计通过计算得出滤波

21、器的系数如下图所示：图3-6 滤波器系数根据以上求出的系数，利用VISIO软件来绘制滤波器的结构图如下图所示：图3-7 滤波器的结构图4 出现的问题及解决办法在这次课程设计当中，由于基础不扎实，出现了很多问题，即有MATLAB软件操作不当，也有知识掌握程度不够出现的各种问题。1，在一开始录制语音信号时，没有修改参数值，导致频率抽样过高，后在老师的提醒之下将频率修改为8000hz。2，在调用音乐文件时没有将文件放在MATLAB的工作文件夹下面，到时文件找不到，还有程序中用到的各种函数都没有放入正确的位置。3，在绘制加噪前后频率对比图时，留白过多，对比不明显，通过axis函数对横纵坐标进行限定4，

22、在使用visio软件绘制滤波器的结构图时，没有合理规划整个页面，导致重新绘制。5，在使用visio绘制结构图时，前面滤波器参数设置出现错误，所以滤波器的系数发生了改变，只能重新修改参数。6，在设计滤波器的过程中，没有准确理解每一条指令代表的含义，导致程序前后不对应，出现很多错误，通过matlab中的错误提示，准确找到错误的那一行代码，进行修改。7，在设计滤波器的过程当中，滤波器的衰减小于开始所设置的值，通过牺牲过渡带和调节过渡带的采样值，即T1和T2来使得衰减大于所设定值As。8，频谱图没有进行去前半段的处理。9，最后听取滤波后声音，没有将其保存，在请教了同学之后，学会了如何保存滤波后的声音文

23、件。10，在设计滤波器的过程当中，抽样点数过高，是的滤波器的阶数很高，导致绘图困难。通过降低抽样的点数，来平衡滤波器阶数和阻带带宽之间的矛盾。5 结束语 这是第三次进行课程设计了，在前面的课程设计当中我们学习到了如何使用matlab，所以对于matlab软件的使用并没有那么陌生了，尽管如此，在使用matlab的过程当中还是出现了很多错误，比如说忘记添加函数文件，参数前后不对应之类的低级错误。在这次课程设计当中，还学习到了另外一个绘图软件，VISIO，相对而言，是比较容易操作的，在这次课程设计中主要用来绘制各类结构图，以及流程图。此次课程设计，让我更深入的了解到频率抽样法以及频率抽样型的滤波器，

24、由开始的无从下手，再翻阅书上的例子，首先决定使用低通的滤波器将噪声滤除，但由于是加入的是一个单频的信号噪声，就决定使用带阻滤波器，在网上查找相关资料时，大多是低通滤波器，没有带阻滤波器的例子，所以就通过书籍信号与系统上高通的例题来进行修改，当然这过程不是一蹴而就的。经过了一次又一次的错误，才得出来最后的模型。还有这也是我和其他同学一起讨论出来的。这次课程设计让我对滤波器的类型有了一个更加完整了解，在设计中也使我对一些概念有了更深刻的认识。例如：在滤波器分类方面，我深刻的了解了低通，高通滤波器与带通，带阻滤波器的特性区别。还有在课程设计中每一次的数据输入都有其重要意义，用MATLAB 编译程序时

25、，可以根据滤波器指标的要求实时知道对滤波器的影响。通过一次次的调试和权衡使滤波器的性能达到最佳。 课程设计不仅要求对滤波器理论的研究，更重要的是培养一种遇到问题解决问题的思维。因为有了这次课程设计，我懂得了书本知识只是实际应用的理论指导。如果仅仅只学习书本知识，不去在实践中运用，那只是停留在只知其表，不知其因的层面。比如在数学计算上，可以将噪声完全滤除，而在这次设计中，若要完全滤除噪音，滤波器的阶数就会增高，在现实生活中是很难实现的，所以噪声是不能完全滤除的。课程设计结束了，我相信这次课程设计对今后的学习是很有帮助的，它让我将理论更好地和实践相结合，提高了动手的能力，也填补了自己学习上的一些不

26、足。这次课程设计能取得成功，不仅仅是我一个人的努力的结果，更离不开指导老师与同学的帮助，在此向老师和同学们表示衷心的感谢。参考文献1 吴镇扬，数字信号处理M.高等教育出版社，20042 张圣勤.MATLAB7.0 实用教程M.北京：机械工程出版社，20063 程佩青数字信号处理教程M北京：清华大学出版社，20024高西全，丁玉美.数字信号处理M. 第三版.西安：西安科大出版社，1994附录：源程序% 程序名称：ht1.m% 程序功能：设计一频率采样法设计的频率采样型的滤波器，然后对一加单凭频噪声信号去噪% 程序作者：贺婷% 最后修改日期：2015-3-1

27、0 x,fs,bits=wavread('D:ht.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2500; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x' y=x+0.1*sin(fn*2*pi*t); %添加噪声sound(y,fs,bits); % 可以明显听出有尖锐的单频啸

28、叫声X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围%设计滤波器fpd=2400;fsd=2450;fsu=2550;fpu=2600;Rp=0.5;As=45; % 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率fcd fcu，和频率间隔dfwp=

29、fcd/fs*2*pi;ws=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率M=41;alpha=(M-1)/2;%设置抽样点数l=0:M-1;T1=0.6025,T2=0.127%过渡带的抽样值NP1=fix(M*fpd/8000);NS1=fix(M*fsd/8000);NS2=ceil(M*fsu/8000);NP2=ceil(M*fpu/8000);%各边界频率的抽样值1Hrs=ones(1,NP1-2),0.6025,0.127,zeros(1,NP2-NP1+1),0.127,0.6025,ones(1,M-2*NP2-

30、5),0.6025,0.127,zeros(1,NP2-NP1+1),0.127,0.6025,ones(1,NP1-3); k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);%求FIR单位冲击响应db,mag,pha,grd,w=freqz_m(h,1);%幅频特性曲线Hr,ww,a,L=hr_type2(h);%画图subplot(2,2,1);stem(l,h);title('滤波器单

31、位冲击响应');xlabel('n');ylabel('y(n)');axis(0 40 -1 1);grid on;subplot(2,2,2);plot(w/pi,pha);title('滤波器相位响应');xlabel('w/pi');ylabel('相位pha');axis(0 1 -4 4);grid on;subplot(2,2,3);plot(w/pi,db);title('滤波器幅度响应');xlabel('w/pi');ylabel('dB'

32、);axis(0 1 -50 20);grid on;%画线X_l=0,0,wp/pi,ws/pi;1,1,wp/pi,ws/pi;Y_l=-As,-Rp,-50,-50;-As,-Rp,20,20; % 在wp，ws，Rp，As处画线以更直观判断设计是否达标，每列参数是每个线条的端点坐标line(X_l,Y_l,'Color','r','LineWidth',2,'LineStyle','-') % 添加线宽为2的红色线条subplot(2,2,4);plot(w/pi,mag);title('滤波器幅度

33、响应');xlabel('w/pi');ylabel('幅度mag');axis(0 1 0 1.2);grid on;%进行滤波处理z=filter(h,1,y); % 用设计好的滤波器对x1进行滤波Z=fft(z);Z=Z(1:length(Z)/2); % 计算频谱取前一半%比较滤波前后的性能%在同一张大图里分别绘制原始信号x，加噪信号y，滤波去噪信号z的时域波形和频谱，以便比较和分析。figure(2);subplot(321);plot(t,x);axis tight;grid on;title('原始语音信号时间x');xlabel('时间t');ylabel('幅度');% 绘制原始录制的语音信号时间波形图subplot(323);plot(t,y);axis tight;grid on;title('加干扰语音信号时间y');xlabel('时间t');ylabel('幅度');% 绘制加噪语音信号时间波形图subplot(325);plot(t,z);axis tight;grid on;title(