FIR IIR 时域滤波滤除高频噪声.doc

数字信号处理课程实验设计报告实验五时域与频域对语音信号进行去噪处理学院：计算机与电子信息学院专业班级：通信工程2009级1班组员及分工：谢锋 0907200224 GUI及FIR设计李德志0934110104 IIR设计及程序调试吕翩翩0826211134 时域去噪设计及报告撰写目录摘要1一实验目的2二实验步骤及分析21制作含高频噪声的音乐信号21.1录制语音信号21.2加噪处理22滤波器设计32.1设计FIR数字滤波器32.1.1FIRDF基本设计原理32.1.2用窗函数设计FIRDF的基本方法42.1.3凯塞窗实现FIRDF62.2设计IIR数字滤波器82.2.1IIRDF基本设计原理82.2.2脉冲响应不变法基本原理92.2.3巴特沃斯滤波器设计基本原理102.2.4用脉冲响应不变法实设计IIRDF112.3设计时域滤波器132.3.1滑动滤波器的基本设计原理133GUI设计18三设计总结181录音中所采集的音乐信号在MATLAB环境中无法打开182GUI中的问题19四结束语19参考文献19图表目录图表 1声音的录制2图表 2原始信号的波形3图表 3原始语音信号的频谱3图表 4幅频特性7图表 5相频特性7图表 6FIR滤波后的波形8图表 7FIR滤波后的频谱8图表 8脉冲响应不变法的映射关系9图表 9幅频特性12图表 10相频特性12图表 11IIR滤波后波形13图表 12IIR滤波后频谱13图表 13滑动平均滤波器的频响特性14图表 14幅频特性16图表 15相频特性16图表 16滤波后的波形16图表 17滤波后的频谱图17图表 18滤波后的波形17图表 19滤波后的频谱图17图表 20本设计的界面总览18摘要 数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。IIR滤波器的首要优点是可在相同阶数时取得更好的滤波效果。但是IIR滤波器设计方法的一个缺点是无法控制滤波器的相位特性。与IIR滤波器相比，FIR的实现是非递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。本次课程设计根据信号的特性，在时域上设计滑动平均滤波器，在频域上分别设计FIR和IIR数字滤波器，对采集的音乐信号进行滤波去噪处理，并绘制出处理前后的时域波形图和频谱图。最后根据处理前后的图形和音乐回放来分析滤波去噪的效果。同时为了更好地表现本设计实现的功能，采用GUIDE设计GUI进行展示。关键词：滤波去噪 滑动平均滤波器 FIR滤波器 IIR滤波器 MATLAB GUI一 实验目的采用MATLAB设计相应滤波器，实现对一段音乐信号进行去噪（在频域分别采用FIR和IIR滤波器，在时域采用书上介绍的方法，即采用滑动平均滤波器设计。）实验要求二 实验步骤及分析1 制作含高频噪声的音乐信号1.1 录制语音信号本实验通过windows自带的录音机采集音频信号，并设置合适的参数，将得到的信号保存为“noisy.wav”。图表 1声音的录制1.2 加噪处理本实验添加的噪声为余弦高频噪声，实现程序如下：Au=0.3; t=0:1/8000:(size(x)-1)/8000; d=Au*cos(2*pi*fn*t); %fn为噪声频率，可在文本框中进行设置x = x+d;得到的复合的原始语音信号的波形和频谱如下所示：图表 2原始信号的波形图表 3原始语音信号的频谱由频谱图可知，语音信号主要分布在2KHz以下的部分，而所加的噪声的频谱是3200Hz。因此滤波器设计必须基于此而设定滤波器的指标参数。2 滤波器设计2.1 设计FIR数字滤波器2.1.1 FIRDF基本设计原理有限长单位脉冲响应（FIR）数字滤波器很容易获得严格的线性相位特性，这对于语音信号处理和数据传输是很重要的。另外，它总是可实现的和稳定的，既可按非递归型结构实现，也可按递归型实现。非递归FIR滤波器的运算量化噪声可做的比较小。FIR数字滤波器的单位脉冲响应仅含有有限个（N个）非零值，是因果的有限长序列。该序列的Z变换为 系统的频率响应为 上式中，叫幅度函数，是含正负符号的实函数，且；叫相位函数2。一个数字滤波器的设计分三步完成：一是技术要求。在设计滤波器之前，必须要有某些技术要求。这些技术要求是由用途决定的；二是近似。一旦技术要求确定后，就要用已学过的各种概念何数学提供一种滤波器的表述，它接近于所给出的一组技术要求；三是实现。上面一步的结果是一个滤波器的表述，它可能是一个差分方程的形式，或者是某一系统函数,或者是某一脉冲响应。根据这个表述要用硬件实现这个滤波器3。2.1.2 用窗函数设计FIRDF的基本方法 FIR滤波器的设计问题，就是要使所设计的滤波器的频率响应逼近所要求的理想滤波器的响应。从单位采样响应序列来看，就是使所设计的滤波器的逼近理单位采样响应序列。窗函数设计法就是从时域出发，根据问题的不同性质，理想滤波器可以分为理想低通、高通、带通和带阻，然后用窗函数截取理想滤波器的单位脉冲响应，得到；以有限长的近似理想的，这样得到的频率响应逼近于理想的。 窗口设计法对窗函数的要求是窗谱的主瓣宽度应尽可能窄，以使设计的滤波器有较陡的过渡带。窗谱的最大副瓣相对于主瓣尽可能小，使设计的滤波器幅频特性中肩峰和余振较小，阻带衰减较大4。理想低通、高通、带通和带阻滤波器满足线性相位的频域特点的频率响应函数在主值区间的表达式如下：理想低通为 理想高通为 理想带通为 理想带阻为 其中为理想滤波器的截止频率，和分别为理想滤波器的通带下截止频率和上截止频率；确保线性相位适于条件,且为实序列。窗函数截断的影响：和的差异是由加窗而引起的。时域加窗对频率响应的影响表现在三方面。一是使理想频率特性不连续点处边沿加宽，形成一个过渡带。过渡带的宽度等于窗谱的主瓣宽度。二是过渡带两旁产生肩峰和阻尼余振，其振荡幅度取决于旁瓣的相对幅度。振荡的多少，取决于旁瓣的多少。三是改变窗函数的长度，将改变窗函数谱的绝对大小和主瓣、副瓣宽度，但不能改变主瓣与副瓣的相对比例，这个相对比列直接决定滤波器通带内的平稳和阻带的衰减。窗函数设计法设计流程如下所示比较，合格则设计完毕，否则重新设计序列加窗截断序列傅氏变换2.1.3 凯塞窗实现FIRDF本实验设计FIR采用凯撒窗设计，这是一种可以调整的窗函数，在实际中得到广泛应用。这个窗是由J.F.Kaiser提出来的，式子如下： 是修正的零阶贝塞尔（Bessel）函数，式子如下： 对全部实数x，它是正的。参数控制最小阻带衰减,而且对于接近最佳可以将它选择成产生不同的过渡带.这种窗函数对同一个M值能够给出不同的过渡带，而这个正是其他的窗函数所不能的。Kaiser窗提供了可以变化的过渡带宽度。由于涉及贝塞尔（Bessel）函数的复杂性，这种窗函数的设计方程不是那么容易导出的。下面是一些经验设计式子：以知，和，参数M和给出过渡带带宽=- 滤波器长度 参数 越大，窗函数越窄，频谱的旁瓣越小，但主瓣宽度增加。根据原始信号的频谱分布和噪声频率，选取合适的fp，fs，rp，rs指标设计凯撒窗，在用GUI实现的过程中，添加了文本框控件，实现了指标参数可调的功能。本实验中滤波器的默认参数为fp=1500,fs= 2000;rp= 5,rs=50,fn=3200;Matalb部分实现代码如下：wp=2*pi*fp;ws=2*pi*fs;wp = wp/Fs;ws = ws/Fs;DB = ws-wp; %计算带宽beta = 0.5842*(rs-21)0.4+0.07886*(rs-21); %控制参数N1 = ceil(rs-8)/2.285/DB); %计算出阶数wc = (wp+ws)/2/pi;B = fir1(N1,wc,kaiser(N1+1,beta); 以下为滤波器特性(指标参数为Fp=1500,Fs=2000,rp=5,rs=50)图表 4幅频特性图表 5相频特性从图中可以看出，凯撒窗的相位特性是线性的，能实现线性相位滤波。原始信号经过该滤波器处理后的结果如下：图表 6FIR滤波后的波形图表 7FIR滤波后的频谱从图中可以看出，经过滤波器后，高频噪声被消除，声音较为纯净。结果分析：对比滤波前后的语音信号，可以明显感觉出对噪声已被滤除，因此，本实验设计的FIR滤波器对语音信号的处理符合设计指标要求。最后使用sound函数对处理后的信号进行回放，与前面的原始音频信号相比，明显没有刺耳的噪声了。2.2 设计IIR数字滤波器2.2.1 IIRDF基本设计原理 IIR数字滤波器的设计方法有两类：间接设计法和直接设计法。简洁设计法是借助模拟滤波器设计方法进行设计的，先根据数字滤波器设计指标设计相应的过渡模拟滤波器，再将过渡模拟滤波器转换成数字滤波器。直接设计法是在时域或频域直接设计数字滤波器。由于模拟滤波器设计理论非常成熟，而且有多种性能优良的典型滤波器可供选择（如巴特沃斯、切比雪夫和椭圆滤波器等），设计公式和图表完善，而且许多实际应用需要模拟滤波器的数字仿真。所以间接设计法得到广泛应用。而直接设计法要求解联立方程组，必须采用计算机辅助设计。我们采用间接设计法来完成对语音信号的滤波处理，具体的先根据语音信号的频谱得到通带截止频率、阻带截止频率、通带最大衰减、阻带最大衰减。然后设计巴特沃斯模拟低通滤波器，采用脉冲响应不变法将模拟滤波器转换为数字滤波器。让语音信号通过滤波器，得到滤波后的效果。2.2.2 脉冲响应不变法基本原理 如果从模拟到数字滤波器我们想要保留脉冲响应的形状，那么就得到一种方法称为脉冲不变响应法的变换方法。 脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序h(n)模仿模拟滤波去的冲击响应ha(t)，使h(n)正好等于ha(t)的采样值，即h(n)=ha(nT) T为采样周期。如以Ha(s)及H(z)分别表示ha(t)的拉式变换及h(n)的z变换，即Ha(s)=Lha(t) H(z)=Zh(n) 则根据采样序列z变换与模拟信号拉式变换的关系，得： 上式表明，采样脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是以前讨论的拉式变换到Z变换的标准变换，即首先对Ha(s)作周期严拓，然后再经过Z=est的映射关系映射到Z平面上。应当指出，Z=est的映射关系表明，S平面上每一条宽为2pi/T的横带部分，都将重叠地映射到Z平面的整个全部平面上。每一横带的左半部分映射到Z平面单位圆以内，每一横带的右半部分映射到Z平面单位圆以外，j轴映射在单位圆上，但j轴上的每一段2pi/T都应于绕单位圆一周，如下图所示：图表 8脉冲响应不变法的映射关系Z=est的映射关系反映的是Ha(s)的周期严拓与H(z)的关系，而不是Ha(s)b本身与H(z)的关系，因此，使用脉冲响应不变法时，从Ha(s)到H(z)并没有一个由S平面到Z平面的简单代数映射关系，即没有一个s=f(z)的代数关系式。另外，数字滤波器的频响也不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期严拓，周期为S=2/T=2fs，即 2.2.3 巴特沃斯滤波器设计基本原理这种滤波器最先由英国工程师斯替芬巴特沃斯(Stephen Butterworth)在1903发表英国无线电工程期刊的一篇论文中提出的。巴特沃斯滤波器的特性是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波得图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。N阶巴特沃斯模拟低通滤波器的振幅平方函数为： 式中，N为整数，称为滤波器的阶数，N越大，通带和阻带的近似性越好，过渡带也越陡。c为3dB截止频率。如下图所示： 巴特沃兹filter 振幅平方函数2.2.4 用脉冲响应不变法实设计IIRDF本设计用间接法设计IIRDF，将数字滤波器的指标转换为模拟滤波器的指标参数，为了获得模拟滤波器的系统函数，选择巴特沃斯滤波器实现。脉冲不变响应法设计IIRDF的流程如下结束回放语音信号比较滤波前后语音信号的波形及频谱设计好的滤波器进行滤波处理验证并进行频谱分析脉冲响应不变法法设计巴特沃斯滤波器设定滤波器性能指标，通带截止频率fp，阻带截止频率fc，通带波纹Ap，阻带波纹As快速傅里叶变换，并且画频谱图语音信号的采集（wavread函数），画时域图开始 Matlab部分实现代码:fp=1500;fs=2000;rp=5;rs=50; %四项指标wp=2*pi*fp;ws=2*pi*fs;N1,wc=buttord(wp,ws,rp,rs,s); %得出巴特沃斯的阶数和3dB截止频率B,A=butter(N1,wc,s); %计算N阶巴特沃斯数字滤波器的系统函数H（z） %的分子和分母多项式的系数向量Bz,Az=impinvar(B,A,Fs); %用脉冲响应不变法设计IIR，将模拟转数字wk=0:pi/512:2*pi;Hz,w=freqz(Bz,Az,wk); %得出频率响应y=filter(Bz,Az,x); %根据系统函数Bz和Az对序列x（n）进行滤波处理Y=fft(y,N); %对滤波信号进行FFT处理绘图部分程序省略。以下为滤波器特性(指标参数为Fp=1500,Fs=2000,rp=5,rs=50)图表 9幅频特性图表 10相频特性运行结果如下图表 11IIR滤波后波形图表 12IIR滤波后频谱从图中可以看出，经过滤波器后，高频噪声被消除，声音较为纯净。结果分析：对比滤波前后的语音信号的频谱图，发现高频噪声已被滤除，因此，本实验设计的IIR滤波器对语音信号的处理符合指标要求。最后使用sound函数对处理后的信号进行回放，与前面的原始音频信号相比，噪声已被消除。2.3 设计时域滤波器在实现时域去噪的过程中，实现的思想是设计一个低通滤波器来滤除高频噪声，本设计采用的是滑动滤波器。2.3.1 滑动滤波器的基本设计原理滑动平均滤波器常用于时域滤波，尽管它很简单，但是对于抑制随机噪声并保留陡峭边沿来说是最优的。 这种滤波器起对输入信号进行平滑的作用，相当于一个低通滤波器，滤除高频分量，而保留低通分量。该滤波器是取输入信号的最近的几个值，进行算术平均，即是差分方程的一个重要应用。一个M项平均的滑动平均滤波器的差分方程： 将代替，便可以得到该滤波器的单位脉冲响应，即 如果让一个快速变化的信号通过M项滑动平均滤波器，得到的是相对变化缓慢的输出信号。平均的项数越多，则得到变化越缓慢的输出信号。而可用卷积法和递推法求解差分方程。可以看出滑动平均滤波的冲击响应就是一个矩形脉冲，通过将输入信号和矩形脉冲进行卷积实现滤波。其频响特性如下： 图表 13滑动平均滤波器的频响特性滑动平均滤波器的设计流程开始语音信号的采集，画时域图 快速傅里叶变换，并且画频谱图设定滤波器，列差分方程验证并进行频谱分析设计好的滤波器进行滤波处理比较滤波前后语音信号的波形及频谱回放语音信号结束Matlab部分实现的代码：hn1=ones(1,5);yn1=conv(hn1,x);N1=length(yn1); %产生横坐标n1=0:N1-1;Hz = fft(hn1,N1); %传输函数Yn1=fft(yn1,N1); %输出函数的傅里叶f1=n1*Fs/N1; 得到的滤波器的特性如下：图表 14幅频特性图表 15相频特性高频噪声频率3200的滤波效果：图表 16滤波后的波形图表 17滤波后的频谱图高频噪声频率3000的滤波效果：图表 18滤波后的波形图表 19滤波后的频谱图结果分析：从幅频特性曲线可知，在某些孤立的频点上衰减很大（例如高频噪声为3200Hz），效果很好，但是在某些频带上衰减小（如高频噪声3000），滤波效果不理想。3 GUI设计本设计中采用GUIDE设计GUI，进行对实验现象的观察和分析。本GUI中用一个按钮来添加需要处理的音频文件，三个按钮来选择对原始信号的处理方式，有IIR，FIR和时域处理三种；为了方便对指标参数的更改，从而更好地对比在同一种滤波处理方式下，设计指标不同，得到不同的滤波效果。在界面的美化方面，在相关控件的CreateFcn函数中添加了图片。整个GUI的效果如下：图表 20本设计的界面总览三 设计总结1 录音中所采集的音乐信号在MATLAB环境中无法打开 出现错误：cannot o