数字信号处理课程设计

1、键入文字陈军 学 院 计算机与通信工程 专 业 通信工程 班 级 通信1303班 学 号 19 学生姓名 陈 军 指导教师 黄红兵 课程成绩 完成日期 2016年3月4日 课程设计成绩评定学 院 计算机与通信工程 专 业 通信工程 班 级 通信1303班 学 号 19 学生姓名 陈 军 指导教师 黄红兵 完成日期 2016 年 3月 4 日 指导教师对学生在课程设计中的评价评分项目优良中及格不及格课程设计中的创造性成果学生掌握课程内容的程度课程设计完成情况课程设计动手能力文字表达学习态度规范要求课程设计论文的质量指导教师对课程设计的评定意见综合成绩 指导教师签字 2016 年 3 月 6 日课

2、程设计任务书计算机与通信工程 院 通信工程 专业 课程名称数字信号处理课程设计时间20152016学年第二学期12周学生姓名 陈军指导老师黄红兵题 目音乐信号滤波去噪使用汉宁窗设计的FIR滤波器主要内容： 用麦克风采集一段音乐信号，绘制波形并观察其频谱，给定相应技术指标，用汉宁窗设计一个满足指标的FIR滤波器，对该音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，绘制出级联型实现结构，根据结果和学过的理论得出合理的结论。要求：（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按

3、要求编写课程设计报告书。应当提交的文件：（1）课程设计学年论文。（2）课程设计附件（主要是源程序）。19陈军 音乐信号滤波去噪使用汉宁窗设计FIR滤波器 第19页 共19页音乐信号滤波去噪使用汉宁窗设计FIR滤波器学生姓名：陈军 指导教师：黄红兵摘 要 本课程设计是对一段加入噪声的音乐信号，用汉宁窗设计出的FIR滤波器进行滤波去噪，分析其前后时域和频域波形。课程设计平台为MATLAB7.0。设计步骤为：首先采集一段音乐信号并观察其频谱，然后设计一个汉宁窗FIR滤波器，最后对该信号进行滤波。信号在进行滤波处理后，观察并记录滤波前后波形和频谱的变化，能够听到滤波后的音乐信号和滤波前相比明显的变得清

4、晰，基本达到了设计目的。最后用visio软件绘制出级联型实现结构，根据结果和学过的理论得出合理的结论。关键词 课程设计；滤波去噪； FIR滤波器；汉宁窗；MATLAB7.01 引 言数字信号处理（Digital Signal Processing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科1。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。目前常用的滤波器设计方法普遍采用Matlab仿真，DSP实现。本课程设计的音乐信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到音乐的读取，音乐信号的抽样、频谱分析，滤波器

5、的设计及音乐信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。1.1课程设计的目的设计一个FIR滤波器，可以有多种方法，窗函数法是设计FIR数字滤波器的最简单也是工程上常用的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、汉明窗(Hamming window)、布拉克曼窗(Blackman windo

6、w)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window) 2。在本次课程设计的目的是如何设计一个汉宁窗FIR滤波器，从而达到对音乐信号滤波的效果。1.2课程设计的要求（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告书。1.3课程设计的平台课程设计的主要设计平台是MATLAB 7.0。MATLAB是矩阵实验室（Matrix Laboratory）的简称5，是美国MathWorks公司出品的

7、商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MathWork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入

8、了对C，FORTRAN，C+ ，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。2 设计原理2.1 FIR滤波器数字滤波是将输入的信号序列3，按规定的算法进行处理，从而得到所期望的输出序列。一个线性位移不变系统的输出序列y(n)和输入x(n)之间的关系，应满足常系数线性差分方程，见公式2-1如下： (2-1)其中，x(n)为输入序列，y(n)为输出序列，为滤波器系数，N是滤波器的阶数。若上式中所有的均为零，则有FIR滤波器的差分方程为： (2-2)对上式进行

9、Z变换得到FIR滤波器的传递函数为： (2-3)由上式可以看出，H(z)是的N-1次多项式，它在z平面内有N-1个零点，同时在原点处有N-1个重极点。N阶滤波器通常采用N个延迟单元、N个加法器与N+1个乘法器，取图2-1 中(a)、(b)两种结构。图2-1 FIR滤波器的一般结构因为FIR滤波器的单位抽样响应是有限长的，所以它永远是稳定的。另外，若对 h(n)提出一些约束条件，那么可以很容易地使 H(z)具有线性相位，这在信号处理的很多领域是非常重要的。FIR滤波器的设计任务，是要决定一个转移函数H(z)，使它的频率响应满足给定的要求。这里所说的要求，除了通带频率、阻带频率及两个带上的最大和最

10、小衰减和外，很重要的一条是保证H(z)具有线性相位。2.2 窗口设计方法窗口函数法也称为傅立叶级数法4。理想的数字滤波器频率特性是无法实现的， FIR的设计就是要寻找一个可以得到的频率特性=来逼近，这相当于用一个可实现的单位脉冲响应 h(n)去逼近一个理想单位脉冲响应.可由理想频率特性 通过傅氏反变换得到， (2-4)一般来说，这样得到的理想单位脉冲响应序列是个无限长序列，因而是非因果的。设有一个截止频率为的理想线性相位低通，延时为，其频率特性是： (2-5)得到： (2-6)这是一个以为中心偶对称的无限长非因果序列，要想用一个有限长的因果序列去逼近它，最简单的方法是截取从0到N-1的一段来表

11、示它，即 h(n)= ；其他N：h(n)=0。同时，为了保证线性相位，还要满足偶对称h(n)=h(N-1-n)。这就好像通过一个窗口观看到的一段，因此h(n)就表示成和一个窗口函数的乘积，这样对h(n)的求解就变为h(n)*，这里的就称为窗口函数, 既然一个频域上的标准的矩形窗口对应于时域是一个无限长的序列, 那么在时域上截取一段势必造成频域的矩形窗口的失真。结果就是截取出的信号也相应失真，为了补偿这种失真，只有改变原来窗口的形状，修正经过时域截取后的窗口失真。窗函数设计方法的基本步骤是：(1) 把作傅里叶逆变换，得；(2) 对自然截短到所需的长度，如2M+1；(3) 将截短后的右移M个采样间

12、隔，得h(n)；(4) 将h(n)乘以合适的窗口，即得所要滤波器的冲击响应，窗函数以n=M对称。利用所求得的单位抽样响应，即可用硬件构成滤波器的转移函数H(z)，也可利用h(n)在计算机上用软件来实现滤波。数据窗在FIR滤波器的窗函数设计中起着重要的作用，它的性能的好坏直接影响着滤波器的过渡带宽度和衰减的大小。对窗函数总的要求，是希望它的频谱中的主瓣尽量窄，边瓣幅度尽量得小，也即它频域的能量主要集中在主瓣内。此外，窗函数还应该满足下列要求以便可以定量地比较各函数的性能。w(n)应是非负的实偶函数，为了使滤波器获得较大的主旁瓣能量比，从对称中心开始w(n)应是非递增的；为了保证滤波器的通带增益为

13、1，应有： (2-7)为了保证滤波器的相位特性不因加窗而改变，一般要求是恒正的；这里给出如下三个频域指标作为窗函数性能的性能参数：(1) 3dB带宽B，它是主瓣归一化幅度()下降到-3dB时的带宽。当数据长度为N时，最大可能的频率分辨率是，则B的单位可以是；(2) 最大旁瓣峰值A(dB)。A越小，由旁瓣引起的振荡幅度越小；(3) 旁瓣峰值渐进衰减速度D(dB/oct)。2.3 汉宁窗汉宁窗（Hanning Window）又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以

14、看，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。 （2-8）根据傅里叶变换的线性性质和调制定理得到 （2-9）当时，则 （2-10）为汉宁窗的幅度响应函数1。3设计步骤3.1流程采集音乐信号后对原始音乐信号进行加噪声，然后设计出合适的滤波器，最后用设计好的滤波器对音乐信号进行滤波去噪，观察并记录结果。设计流程如图3-1所示。 在Windows下截取一段音乐并将格式改为.wav加入单频噪声对音乐信号进行频谱分析，画出时域和频域波形图 用汉宁窗设计FIR滤波器进行滤波画出滤波前后波形和频谱并且进行比较分析在visi

15、o软件上绘制出级联型实现结构开始结束图3-1 设计流程图3.2 录制音乐信号用麦克风录制一段格式.wav的音乐信号，时间为5s左右。音乐信号属性如图3-2所示。图3-2 音乐信号属性在Matlab软件平台下，利用函数wavread对录入音乐信号进行采样，记住采样频率和采样点数，并试听。x,fs,bits=wavread('C:Documents and SettingsAdministrator桌面319陈军.wav');%输入参数为文件的全路径和文件名sound(x,fs,bits);%按指定的采样率和每样本编码位数回放采集完成后，在音乐信号中加入一个单频

16、噪声（频率为1900，振幅为1的正弦信号），试听加噪信号，明显听出有尖锐啸叫声。N=length(x); %计算信号x的长度fn=1900; %单频噪声频率t=0:1/fs:(N-1)/fs;%计算时间范围，样本数除以采样频率y=x'+sin(fn*2*pi*t);%加入一个单频噪声sound(y,fs,bits);%可以明显听出有尖锐的单频啸叫声分别对原始信号和加噪信号进行傅里叶（FFT）变换，对比分析两种信号在时域、频谱的幅度，幅度谱如图3-3所示。图3-3 原始、加噪信号时域图和频谱图根据上图分析，由原始音乐信号加入单频噪声后得到的加噪信号，在时域上，振幅比原始信号稍大；在频域上

17、，在频率为1900处幅度谱多出一条高频冲击，可知加入噪声信号成功。实现程序如下：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;%计算频谱频率范围3.3 滤波器设计根据加入噪声信号的频率，选择合适的滤波器。本次课程设计选用带阻滤波器，它的各项性能指标为：下阻带边缘：fpd=1300;fsd=1800; 上阻带边缘 fsu=2000;fpu=2500;通带波纹及最小衰减：Rp=1;As=37。根据上述指标，使

18、用汉宁窗设计FIR滤波器。Matlab实现程序如下：fpd=1300;fsd=1800;fsu=2000;fpu=2500;%FIR滤波器的上下截止频率Rp=1;As=37;%带阻滤波器设计指标 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu);%计算上下边带中心频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率 wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(6.2*pi/dw)+1;%计

19、算汉宁窗设计该滤波器时需要的阶数 n=0:M-1;%定义时间范围w_han=hanning(M);%产生M阶的汉宁窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应，函数ideal_lp见附录2 h_bs=w_han'.*hd_bs;%用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1);%调用自编函数计算滤波器的频率特性，函数freqz_m见附录2设计好滤波器后，画出其幅度响应、相位响应以及脉冲响应，如图3-4所示。 图3-4滤波器幅度和

20、相位响应图3.4 信号滤波处理设计好滤波器后，我们将要对音乐信号进行滤波去噪。将原始信号与滤波器特征方程相乘，得到滤波后信号。分别绘制原始信号x，加噪信号y，滤波去噪信号y_fil的时域波形和频谱，以便比较和分析。Matlab实现程序如下：y_fil=filter(h_bs,1,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半绘制得原始信号x，加噪信号y，滤波去噪信号Y的时域波形和频谱图如图3-5所示。图3-5 滤波前和滤波后音乐信号的波形及频谱的对比得出滤波前和滤波后音乐信号的波形及频谱图后，画一个

21、纵坐标为对数坐标的去噪前后的幅度频谱图对比图，便于比较滤波效果。Matlab实现程序如下：Y1=20*log10(Y); % 画一个纵坐标为对数坐标的去噪前的幅度频谱图Y2=20*log10(Y_fil); % 画一个纵坐标为对数坐标的去噪后的幅度频谱图以纵坐标为对数坐标的去噪前后的幅度频谱图对比图如图3-6所示图3-6 以纵坐标为对数坐标的去噪前后的幅度频谱图对比图3.5结果分析通过观察图3-5，滤波后的音乐信号发生了衰减，说明滤波器起到了滤波作用，同时通过幅频对比，可以看出滤波器滤掉了一部分频率范围内的信号。观察图3-6，明显发现频率范围内滤波后音乐信号噪声减弱。语音信号经过FIR滤波器的

22、滤除噪声的处理，在Matlab中，函数sound可以对声音进行回放。sound(y_fil,fs,bits); % 按指定的采样率和每样本编码位数回放分别听原始音乐和滤波后的音乐信号，发现滤波后的音乐信号噪声减小了，同时原始信号强度稍有减弱，基本达到了滤波的效果。3.6滤波器结构设计调用函数C,B,A=dir2cas(h,1)计算相应结构的滤波器系数。Matlab实现程序如下：C,B,A=dir2cas(h_bs,1) %将FIR滤波器的直接型转级联型得到的相应结构的滤波器系数如图3-7所示图3-7 相应结构的滤波器系数图绘制的级联型滤波器结构图如图3-8所示图3-8 级联型滤波器结构图4出现

23、的问题及解决方法4.1出现的问题（1）原始信号和加噪声后的波形和频谱图显示不出来；（2）自编函数调用不到；（3）绘制出滤波后的波形，发现FIR滤波器没有滤掉单频噪声。4.2 解决方法（1）重新设置噪声信号的幅度，取到合适的幅度大小即可；（2）把自编函数和源程序放在同一个文件夹里即可调用；（3）单频噪声的频率改动后，FIR滤波器的频率没有改动。所以单频噪声的频率也应该自己先定义，FIR滤波器的截止频率应该以单频噪声的频率为中心，这样重新运行后，结果正确5结束语本次课程设计中，我们对上学期的专业课数字信号处理进行了复习，熟悉了Matlab 软件操作环境以及应用方法，然后通过Matlab软件处理音乐

24、信号，让我们更清楚的认识到数字信号的含义我所做的课程设计主要是是对一段音乐信号，加入单频噪声后，用汉宁窗函数法设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理，并且分析对比前后时域和频域波形的程序设计。虽说是比较的简单，但是也会出现一些错误，分析并总结这些错误主要是语句上的一些用法出现错误，以及没有及时地掌握老师所提醒的部分。课程设计是我们运用所学知识，动手实践的一个很好的机会。它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。在实物

25、的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。在这次课程设计中，胡双红、黄红兵等老师给了我们很大的帮助和启发，随时随地为我们答疑，引导我们思考，而不是一味的接受。在这里向给我们以指导的老师和热心帮助我们的同学表示衷心的敬意和感谢！参考文献1 程佩青数字信号处理教程M北京：清华大学出版社，20022 维纳·K·恩格尔，约翰·G·普罗克斯.数字信号处理M .西安：西安交通大学出版社，2002.63 陈金鹰.DSP技术及应用M. 北京：机械工业出版社.2004.64 孙宗瀛. DSP原理设计与应用M. 北京：清华大学出版社，20

26、025 薛年喜主编.MATLAB 在数字信号处理中的应用M北京：清华大学出版社 附录1：语音信号滤波去噪设计源程序清单% 程序名称：voicenoisefilter.m% 程序功能：采用基于汉宁窗的窗口设计法，设计FIR滤波器对含噪语音进行滤波去噪处理。% 程序作者：陈军% 最后修改日期：2016-3-1x,fs,bits=wavread('C:Documents and SettingsAdministrator桌面319陈军.wav');%输入参数为文件的全路径和文件名sound(x,fs,bits);%按指定的采样率和每样本编码位数回放N=length

27、(x);%计算信号x的长度fn=1900;%单频噪声频率t=0:1/fs:(N-1)/fs;%计算时间范围，样本数除以采样频率y=x'+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);%截取前半部分deltaf=fs/ N;%计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf;%计算频谱频率范围figure(1)subplot(2,2

28、,1);plot(t,x);axis(0 6 -1 1)xlabel('时间（单位：s）');ylabel('幅度');title('原始音乐信号');subplot(2,2,2);plot(f,X);axis(0 4000 0 1000);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('原始音乐信号幅度谱图');subplot(2,2,3);plot(t,y);axis(0 6 -2&

29、#160;2)xlabel('时间（单位：s）');ylabel('幅度');title('加入单频噪声音乐信号');subplot(2,2,4);plot(f,Y);axis(0 4000 0 1000)xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入单频噪声信号幅度谱图');fpd=1300;fsd=1800;fsu=2000;fpu=2500; % FIR滤波器的上下截止频率Rp=1;As=37;%

30、0;带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu);  % 计算上下边带中心频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;  % 将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(6.2*pi/dw)+1;         

31、 % 计算汉宁窗设计该滤波器时需要的阶数n=0:M-1; %定义时间范围w_ham=hanning(M);%产生M阶的汉宁窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_ham'.*hd_bs; %用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性y_fil=filter(h_bs,1,y);%用设计好的滤波器对y进行滤波fi

32、gure(2)subplot(2,2,1);plot(w/pi,db);axis(0.2 0.8 -60 1);grid onxlabel('w/pi');ylabel('db');title('滤波器幅度响应图');line(0,0.8,-As,-As,'color','r','linestyle','-','LineWidth',2);line(0,0.8,-Rp,-Rp,'color','r'

33、;,'linestyle','-','LineWidth',2);line(wsd/pi,wsd/pi,-60,10,'color','r','linestyle','-','LineWidth',2);line(wsu/pi,wsu/pi,-60,10,'color','r','linestyle','-','LineWidth',2);subplot(2,2,2);plot(w/pi,

34、mag);axis(0.2 0.8 0 1.15);grid onxlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图')subplot(2,2,3);plot(w/pi,pha);axis(0 1 -3 3);grid onxlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图');subplot(2,2,4);plot(n,h_bs);gr

35、id onxlabel('n');ylabel('h_bs(n)');title('滤波器脉冲响应图')Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2);% 计算频谱取前一半Y1=20*log10(Y);Y2=20*log10(Y_fil);figure(3)subplot(3,2,1);plot(t,x);axis(0 5 -2 2)xlabel('时间（单位：s）');ylabel('幅度');title('原始音乐信号x'

36、);subplot(3,2,2);plot(f,X);axis(0 5000 0 1000)xlabel('频率（单位：Hz）');ylabel('幅度谱');title('原始音乐信号幅度谱图X');subplot(3,2,3);plot(t,y);axis(0 5 -2 2)xlabel('时间（单位：s）');ylabel('幅度');title('加入单频噪声音乐信号x');subplot(3,2,4);plot(f,Y);axis(0 5000 0 1000)xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入单频噪声信号幅度谱图X');subplot(3,2,5);plot(t,y_fil);axis(0&#