滤波去噪课程设计

音乐信号滤波去噪—使用汉宁窗设计FIR滤波器学生姓名：陈浩指导教师：胡双红摘要本课程设计是对一段加入噪声的音乐信号，用汉宁窗设计出的FIR滤波器进行滤波去噪，分析其前后时域和频域波形。课程设计平台为MATLAB7.0。设计步骤为：首先采集一段音乐信号并观察其频谱，然后设计一个汉宁窗FIR滤波器，最后对该信号进行滤波。信号在进行滤波处理后，观察并记录滤波前后波形和频谱的变化，能够听到滤波后的音乐信号和滤波前相比明显的变得清晰，基本达到了设计目的。关键词课程设计；滤波去噪；FIR滤波器；汉宁窗；MATLAB7.01引言数字信号处理（DigitalSignalProcessing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。目前常用的滤波器设计方法普遍采用Matlab仿真，DSP实现。本课程设计的音乐信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到音乐的读取，音乐信号的抽样、频谱分析，滤波器的设计及音乐信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。1.1课程设计的目的设计一个FIR滤波器，可以有多种方法，窗函数法是设计FIR数字滤波器的最简单也是工程上常用的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗（Rectangularwindow）＞三角窗（Triangularwindow)、汉宁窗(Hanningwindow)、汉明窗(Hammingwindow)、布拉克曼窗(Blackmanwindow)、切比雪夫窗(Chebyshevwindow)、巴特里特窗(Bartlettwindow)及凯塞窗(Kaiserwindow)o在本次课程设计的目的是如何设计一个汉宁窗FIR滤波器，从而达到对音乐信号滤波的效果。1.2课程设计的要求滤波器指标必须符合工程实际。设计完后应检查其频率响应曲线是否满足指标。处理结果和分析结论应该一致，而且应符合理论。独立完成课程设计并按要求编写课程设计报告书。1.3课程设计的平台课程设计的主要设计平台是MATLAB7.0。MATLAB是矩阵实验室(MatrixLaboratory)的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MathWork也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。2设计原理2.1FIR滤波器

数字滤波是将输入的信号序列，按规定的算法进行处理，从而得到所期望的输出序列。一个线性位移不变系统的输出序列y(n)和输入x(n)之间的关系，应满足常系数线性差分方程，见公式2-1如下：(2-1)y(n)=习bx(n一i)-l^a,y(n-i)(2-1)i=0i=1其中，x(n)为输入序列，y(n)为输出序列，气和匕为滤波器系数，N是滤波器的阶数。若上式中所有的b若上式中所有的b均为零，

k则有FIR滤波器的差分方程为:(2-2)y(n)=0一ax(n-k)(2-2)k=0对上式进行Z变换得到FIR滤波器的传递函数为:(z)=(z)=华=0

Xk)i=0b^z-k(2-3)由上式可以看出，H(z)是z-1的N-1次多项式，它在/平面内有N-1个零点，同时在原点处有N-1个重极点。N阶滤波器通常采用N个延迟单元、N个加法器与N+1个乘法器，取图2-1中(a)、(b)两种结构。图2-1图2-1FIR滤波器的一般结构因为FIR滤波器的单位抽样响应是有限长的，所以它永远是稳定的。另外，若对h(n)提出一些约束条件，那么可以很容易地使H(z)具有线性相位，这在信号处理的很多领域是非常重要的。FIR滤波器的设计任务，是要决定一个转移函数H(z)，使它的频率响应满足给定的要求。这里所说的要求，除了通带频率°八阻带频率及两个带上的最大和最小衰减6p和°s外，很重要的一条是保证H(z)具有线性相位。窗口函数法也称为傅立叶级数法。理想的数字滤波器频率特性H(ejw)是无法实现的，FIR的设计就是要寻找一个可以得到的频率特性H(ejw)=Vh(n)e-w来逼近n=0H(eJW)，这相当于用一个可实现的单位脉冲响应h(n)去逼近一个理想单位脉冲响应h(n).h(n)可由理想频率特性H{ew)通过傅氏反变换得到，dddhd(n)=2-^^七(ejw)ewd3(2-4)一般来说，这样得到的理想单位脉冲响应序列勺(n)是个无限长序列，因而是非因果的。设有一个截止频率为3c的理想线性相位低通，延时为T，其频率特性是:H侦)="0书<3(2-5)d[03<p<-得到：sin3(n-T)]h(n)=(*)——-8<n<8(2-6)这是一个以n=T为中心偶对称的无限长非因果序列，要想用一个有限长的因果序列去逼近它，最简单的方法是截取n从。到N-1的一段来表示它，即h(n)=hd(n)(0<n<N-1)；其他N：h(n)=0。同时，为了保证线性相位，还要满足偶对称h(n)=h(N-1-n)。这就好像通过一个窗口观看到的一段h'n)，因此h(n)就表示成h「n)和一个窗口函数的乘积，这样对h(n)的求解就变为h(n)=h*)*W，这里的W就称为窗口函数，既然一个频域上的标准的矩形窗口对应于时域是一个无限长的序列，那么在时域上截取一段势必造成频域的矩形窗口的失真。结果就是截取出的信号也相应失真，为了补偿这种失真，只有改变原来窗口的形状，修正经过时域截取后的窗口失真。窗函数设计方法的基本步骤是：把七^e>扁傅里叶逆变换，得h'n)；对、(n)自然截短到所需的长度，如2M+1；将截短后的气(n)右移M个采样间隔，得h(n)；

(4)将h(n)乘以合适的窗口，即得所要滤波器的冲击响应，窗函数以n=M对称。利用所求得的单位抽样响应，即可用硬件构成滤波器的转移函数H(z)，也可利用h(n)在计算机上用软件来实现滤波。数据窗在FIR滤波器的窗函数设计中起着重要的作用，它的性能的好坏直接影响着滤波器的过渡带宽度和衰减的大小。对窗函数总的要求，是希望它的频谱中的主瓣尽量窄，边瓣幅度尽量得小，也即它频域的能量主要集中在主瓣内。此外，窗函数还应该满足下列要求以便可以定量地比较各函数的性能。w(n)应是非负的实偶函数，为了使滤波器获得较大的主旁瓣能量比，从对称中心开始w(n)应是非递增的；为了保证滤波器的通带增益为1，应有：①(0)=—1W(沁)w=1(2-7)2冗-1为了保证滤波器的相位特性不因加窗而改变，一般要求W(加)是恒正的；这里给出如下三个频域指标作为窗函数性能的性能参数:⑴3dB带宽⑴3dB带宽B，它是主瓣归一化幅度(20lo)下降到-3dB时的带宽。当数据长度为N时，最大可能的频率分辨率是△①=21N，则B的单位可以是Ao；最大旁瓣峰值A(dB)。A越小，由旁瓣引起的振荡幅度越小；旁瓣峰值渐进衰减速度D(dB/oct)。2.3汉宁窗汉宁窗(HanningWindow)又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个sin(t)型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了n/T从而使旁瓣互相抵消，消去高频干扰和漏能。可以看，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗.但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。*(n)=*(n)=0-51-(21)cosIn-1J=0.5]R。)C°Sl兽j(2-8)根据傅里叶变换的线性性质和调制定理得到WCm)=FTW(n也WCo)e-ja(N-1)WCm)=FTW(n也WCo)e-ja(N-1)/2hn0.5WQ)+0.25WRgRghng(2兀)①+IN-1)+0.25WRg(2兀)①一——NJe-Mn-1)/2(2-9)WnQ)=0.5WrQ)+0.25Wr(2号(2兀）①++0.25W①-VN)RgVN-1)（2-10）whn（^）为汉宁窗的幅度响应函数[1]3设计步骤3.1流程采集音乐信号后对原始音乐信号进行加噪声，然后设计出合适的滤波器，最后用设计好的滤波器对音乐信号进行滤波去噪，观察并记录结果。设计流程如图3-1所示。结束结束图3-1设计流程图3.2录制音乐信号用麦克风录制一段格式.wav的音乐信号，时间为2s左右。音乐信号属性如图3-2所示。日另存为序)..一|名称通:日另存为序)..一|名称通:图3-2音乐信号属性在Matlab软件平台下，利用函数wavread对录入音乐信号进行采样，记住采样频率和采样点数，并试听。>>[x,fs,bits]=wavread('C:\DocumentsandSettings\student\桌面\cuinaio.wav');>>sound(x,fs,bits);%按指定的采样率和每样本编码位数回放采集完成后，在音乐信号中加入一个单频噪声(频率为2000，振幅为1的正弦信号)，试听加噪信号，明显听出有尖锐啸叫声。>>N=length(x);%计算信号x的长度>>fn=2000;%单频噪声频率>>t=0:1/fs:(N-1)/fs;%计算时间范围，样本数除以采样频率>>x=x';y=x+sin(fn*2*pi*t);>>sound(y,fs,bits);%明显听出有尖锐的单频啸叫声分别对原始信号和加噪信号进行傅里叶(FFT)变换，对比分析两种信号在时域、频谱的幅度，幅度谱如图3-3所示。

图3-3原始、加噪信号时域图和频谱图根据上图分析，由原始音乐信号加入单频噪声后得到的加噪信号，在时域上，振幅比原始信号稍大；在频域上，在频率为2000处幅度谱多出一条高频冲击，可知加入噪声信号成功。实现程序如下：>>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=1710;fsd=1755;上阻带边缘：fsu=1845;fpu=1890;通带波纹及最小衰减：Rp=1;As=44。根据上述指标，使用汉宁窗设计FIR滤波器。Matlab实现程序如下：>>fpd=1710;fsd=1755;fsu=1845;fpu=1890;Rp=1;As=44;%带阻滤波器设计指标>>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;;%计算汉宁窗设计该滤波器时需要的阶数>>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所示。Bd-4600DH5602滤波器幅度响应图-1.24860.10.20.30.4滤波器幅度响应图15Ogm度幅20-2Bd-4600DH5602滤波器幅度响应图-1.24860.10.20.30.4滤波器幅度响应图15Ogm度幅20-20.5w/pi位相w/pi滤波器相位响应图00w/pi滤波器脉冲响应图n图3-4滤波器幅度和相位响应图设计完成滤波器后，求出滤波器的频率特性>>As=-round(max(db(wsd/delta_w+1:1:wsu/delta_w+1)))%实际阻带波As=46>>Rp=-(min(db(1:1:wcd/delta_w+1)))%实际通带波Rp=0.0238由上式可知实际设计滤波器相应参数符合要求，滤波器设计成功。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_fil的时域波形和频谱图如图3-5所示。1r-111c01234时间(单位：s)原始音乐信号x谱度幅O度幅010002000300040005000频率(单位：Hz)加入单频噪声信号幅度谱图Y原始音乐信号幅度谱图X10000051r-111c01234时间(单位：s)原始音乐信号x谱度幅O度幅010002000300040005000频率(单位：Hz)加入单频噪声信号幅度谱图Y原始音乐信号幅度谱图X100000501I,b加入单频噪声音乐信号y11000滤波后音乐信号yii谱度幅42OO度幅-10123时间(单位：s)4谱度幅O度幅500-0—:010002000300040005000频率（单位：Hz）滤波后音乐信号Yfil100000510002000300040005000频率（单位：Hz）3.5结果分析通过观察图3-5，滤波后的音乐信号发生了衰减，说明滤波器起到了滤波作用，同时通过幅频对比，可以看出滤波器滤掉了一部分频率范围内的信号。语音信号经过FIR滤波器的滤除噪声的处理，在Matlab中，函数sound可以对声音进行回放。>>sound(y_fil,fs,bits);%按指定的采样率和每样本编码位数回放分别听原始语音和滤波后的语音信号，发现滤波后的语音信号噪声减小了，同时原始信号强度稍有减弱，基本达到了滤波的效果。4出现的问题及解决方法4.1出现的问题录入的wav格式音乐信号无法用Matlab中Sound函数播放；原始信号和加噪声后的波形和频谱图图显示不出来；绘制出滤波后的波形，发现FIR滤波器没有滤掉单频噪声。4.2解决方法把音乐信号改为单声道信号后能放出音乐；重新设置噪声信号的幅度，取到合适的幅度大小即可；单频噪声的频率改动后，FIR滤波器的频率没有改动。所以单频噪声的频率也应该自己先定义，FIR滤波器的截止频率应该以单频噪声的频率为中心，这样重新运行后，结果正确。5结束语本次课程设计中，我们对上学期的专业课《数字信号处理》进行了复习，熟悉了Matlab软件操作环境以及应用方法，然后通过Matlab软件处理音乐信号，让我们更清楚的认识到数字信号的含义我所做的课程设计主要是是对一段音乐信号，加入单频噪声后，用汉宁窗函数法设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理，并且分析对比前后时域和频域波形的程序设计。虽说是比较的简单，但是也会出现一些错误，分析并总结这些错误主要是语句上的一些用法出现错误，以及没有及时地掌握老师所提醒的部分。课程设计是我们运用所学知识，动手实践的一个很好的机会。它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。在实物的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。在这次课程设计中，胡双红老师给了我们很大的帮助和启发，随时随地为我们答疑，引导我们思考，而不是一味的接受。在这里向给我们以指导的老师和热心帮助我们的同学表示衷心的敬意和感谢！参考文献[1]程佩青.数字信号处理教程[M].北京：清华大学出版社，2002⑵维纳・K・恩格尔，约翰・G・普罗克斯.数字信号处理[M].西安：西安交通大学出版社，2002.6陈金鹰.DSP技术及应用[M].北京：机械工业出版社.2004.6孙宗瀛.DSP原理设计与应用[M].北京：清华大学出版社，2002薛年喜主编.MATLAB在数字信号处理中的应用】M].北京：清华大学出版社附录1：Matlab源程序>>[x,fs,bits]=wavread('C:\DocumentsandSettings\student\桌面\cuiniao.wav');>>sound(x,fs,bits);%按指定的采样率和每样本编码位数回放>>N=length(x);%计算信号x的长度>>fn=1800;%单频噪声频率>>t=0:1/fs:(N-1)/fs;%计算时间范围，样本数除以采样频率>>x=x';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;%计算频谱频率范围>>subplot(2,2,1);plot(t,x);axis([06-11])xlabel('时间(单位：s)');ylabel('幅度');title('原始音乐信号');subplot(2,2,2);plot(f,X);axis([0800008000]);xlabel('频率(单位：Hz)');ylabel('幅度谱');title('原始音乐信号幅度谱图');subplot(2,2,3);plot(t,y);axis([06-22])xlabel('时间(单位：s)');ylabel('幅度');title('加入单频噪声音乐信号');subplot(2,2,4);plot(f,Y);axis([0800008000])xlabel('频率(单位：Hz)');ylabel('幅度谱');title('加入单频噪声信号幅度谱图');>>fpd=1900;fsd=1950;fsu=2050;fpu=2100;Rp=1;As=50;%带阻滤波器设计指标>>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;%计算汉宁窗设计该滤波器时需要的阶数>>n=0:M-1;%定义时间范围>>w_han=hanning(M);%产生M阶的汉宁窗>>hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);%调用自编函数计算理想带阻滤波器的脉冲响应>>h_bs=w_han'.*hd_bs;%用窗口法计算实际滤波器脉冲响应>>[db,mag,pha,grd,w]=freqz_m(h_bs,1);%调用自编函数计算滤波器的频率特性>>subplot(2,2,1);plot(w/pi,db);axis([00.5-501]);gridonxlabel('w/pi');ylabel('db');title('滤波器幅度响应图')；subplot(2,2,2);plot(w/pi,mag);axis([00.501.15]);gridonxlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图')subplot(2,2,3);plot(w/pi,pha);axis([01-33]);gridonxlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图');subplot(2,2,4);plot(n,h_bs);gridonxlabel('n');ylabel('h_bs(n)');title('滤波器脉冲响应图')axis([01500-0.11]);xlabel('n');ylabel('h(n)');grid>>As=-round(max(db(wsd/delta_w+1:1:wsu/delta_w+1)))%%实际阻带波As=46y_fil=filter(h,1,y);%用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:N/2);%计算频谱取前一半subplot(3,2,1);plot(t,x);axis([04-11])xlabel('时间(单位：s)');ylabel('幅度');title('原始音乐信号x');subplot(3,2,2);plot(f,X);axis([05000010000])xlabel('频率(