FIRIIR时域滤波滤除高频噪声.docx

1、FIRIIR时域滤波滤除高频噪声 数字信号处理综合实验报告题 目： FIR-IIR-时域滤波滤除高频噪声 姓 名： 张宝元 学 号： 20141060040 年 级： 2014级 专 业： 电子信息工程 时 间： 2016年12月25日 摘要 数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。IIR滤波器的首要优点是可在相同阶数时取得更好的滤波效果。但是IIR滤波器设计方法的一个缺点是无法控制滤波器的相位特性。与IIR滤波器相比，FIR的实现是非

2、递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。本次课程设计根据信号的特性，在时域上设计滑动平均滤波器，在频域上分别设计FIR和IIR数字滤波器，对采集的音乐信号进行滤波去噪处理，并绘制出处理前后的时域波形图和频谱图。最后根据处理前后的图形和音乐回放来分析滤波去噪的效果。关键词：滤波去噪 滑动平均滤波器 FIR滤波器 IIR滤波器 MATLAB GUI一、设计基本原理（一）MATLAB软件设计平台简介MATLAB是矩阵实验室（Matrix Laborator

3、y）的简称，是美国MathWorks公司研发的商业软件，用于数据可视化、算法开发、数值计算数据分析以及数值计算的交互式环境和高级技术计算语言中，其中主要包含Simulink和MATLAB两大部分。MATLAB是由美国mathworks公司发布的主要面对可视化、交互式程序设计以及科学计算的高科技计算环境。它将矩阵计算、科学数据、可视化、非线性动态系统的建模和仿真以及数值分析等一系列强大功能集成在一个简单方便使用的可视窗口中，为工程设计、科学研究以及那些必须从事有效数值计算的一系列科学领域提供了全方面的解决办法，并在极大程度上舍弃了C、Fortran等传统非交互式程序设计语言的编辑模式，从而体现了

4、当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称三大数学软件。MATLAB可以进行绘制函数和数据、矩阵运算、连接其他编程语言的程序、创建用户界面、实现算法等，主要应用于控制设计、工程计算、图像处理、信号处理与通讯、金融建模设计与分析、信号检测等领域。同时利用附加的工具箱来扩展 MATLAB 环境，其中专用的MATLAB函数集可以解决一些应用领域特定类型内无法解决的问题。MATLAB的主要特点如下：(1) 程序的可移植性良好应用于其他程序。(2) 程序限制宽泛，程序设计自由。有大量已经系统定义的函数可直接应用，并且能够用户自定义函数。(3) 语言简洁，使用灵活

5、方便，库函数相当丰富。(4) 源程序向大众开放。用户可灵活的对源文件进行修改以及加入自己的设计语音构成新的工具箱。(5) 最后MATLAB的一个重要特点是功能强大的工具箱。MATLAB包含两个重要的部分：核心部分和各种可选的工具箱。（二）FIR滤波器设计的基本原理1.2.1数字滤波器的概念 数字滤波器(DigitalFilter，简称为DF)是指用来对输入信号进行滤波的硬件和软件。所谓数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器和模拟滤波器相比，因为信号的形式和实现滤波的方式不同，数字滤波器具有比模拟滤波器精度

6、高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。一般用两种方法来实现数字滤波器:一是采用通用计算机，把滤波器所要完成的运算编成程序通过计算机来执行，也就是采用计算机软件来实现;二是采用实际专用的数字处理硬件。1 1.2.2 IIR和FIR滤波器数字滤波器在数字信号处理的各种应用中发挥着十分重要的作用。它是通过对采样数据信号进行数学运算处理来达到滤波的目的。数字滤波器从实现的网络结构或者从单位脉冲响应可分为无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。FIR滤波器的设计方法和IIR滤波器的设计方法有很大的不同，FIR滤波器设计任务是选择有限长度的h(n),使传输函数H（ejw

7、）满足技术要求，FIR数字滤波器设计的方法有三种，第一种是窗函数法，第二种是频率采样法，第三种是切比雪夫等波逼近法。从性能上说，IIR滤波器以非线性相位为代价以较低的阶数获得较高的选择性。而FIR滤波器想要获得相同的选择性阶数是IIR滤波器的5-10倍，结果成本较高、信号时延也较大：从结构上说，IIR采用递归结构，FIR采用非递归结构；从设计工具上说；IIR可以借助于模拟滤波器的成果，FIR滤波器一般采用没有封闭形式的设计公式；从使用场合上来看，在对相位要求不敏感的场合，如语音通讯等，选用IIR较为合适，可以充分发挥经济高效的特点。对图像处理、数据传输等以波形携带信息的系统，使用FIR较好。1

8、.2.3设计IIR数字滤波器的基本思想设计IIR数字滤波器的方法主要有基于冲激响应不变法的IIR数字滤波器设计，基于双线性Z变换法的IIR数字滤波器设计，数字高通、带通及带阻IIR滤波器设计，基于MATLAB函数直接设计IIR数字滤波器。本实验中采用双线性变换法变换的巴特沃思数字滤波器。1.2.3.1巴特沃思低通数字滤波器(1)选择来自于D盘的“ding.wav”声音作为语音信号（用如下语句调用x,FS,bits=wavread(D:ding.wav)）。(2)给信号加一个大频率的噪声（取噪声频率远大于语音信号的最大频率），产生污染信号。(3)设计一个巴特沃思低通滤波器，通带范围包括语音信号，

9、阻带频率设定为小于噪声信号频率。(4)将设计好的巴特沃思低通滤波器滤除被噪声污染后的语音信号。还原语音信号。1.2.3.2巴特沃思高通数字滤波器(1)选择来自于D盘的“ding.wav”声音作为语音信号（用如下语句调用x,FS,bits=wavread(D:ding.wav)）。(2)给信号加一个小频率的噪声（取噪声频率远小于语音信号的最小频率），产生污染信号。(3)设计一个巴特沃思高通滤波器，通带范围包括语音信号，阻带频率设定为大于噪声信号频率。(4)将设计好的巴特沃思低通滤波器滤除被噪声污染后的语音信号。还原语音信号。1.2.3.3巴特沃思带通数字滤波器(1)选择来自于D盘的“ding.w

10、av”声音作为语音信号（用如下语句调用x,FS,bits=wavread(D:ding.wav)）。(2)给信号加一个小频率或大频率的噪声（取噪声频率远小于语音信号的最小频率或大于语音信号的最大频率），产生污染信号。本实验取小频率的噪声信号。(3)设计一个巴特沃思带通滤波器，通带范围包括语音信号，阻带频率设定为不包括噪声信号频率。(4)将设计好的巴特沃思带通滤波器滤除被噪声污染后的语音信号。还原语音信号。1.2.4设计FIR滤波器的基本思想FIR滤波器通常采用窗函数方法来设计。正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以

11、下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Bartlett window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。本实验中选用布拉克曼窗(Bartlett window)设计滤波器。1.2.4.1 凯泽窗低通滤波器(1)选择来自于桌面的“qwe.wav”声音作为语音信号（用如下语句调用y,fs=audioread(C:UsersAdministrato

12、rDesktopqwe.wav);）。(2)给信号加一个大频率的噪声（取噪声频率远大于语音信号的最大频率），产生污染信号。(3) 设计一个凯泽窗低通滤波器，通带范围包括语音信号，阻带频率设定为小于噪声信号频率。(4)将设计好的凯泽窗低通滤波器滤除被噪声污染后的语音信号。还原语音信号。（三）语音信号的采样理论依据1采样频率采样频率是计算机每秒钟采集的声音样本数，是描述声音文件的音调和音质，是用来衡量声卡和声音文件的质量标准。采样频率越高，对声音波形的表示也越精确。同时采样频率与声音频率也有一定的关系，由奈奎斯特定理可知只有采样频率大于声音最高频率的2倍时，才能把用数字信号来表示的声音还原成原始语

13、音信号。因此采样频率是用来衡量声卡采集、记录和还原声音文件的质量标准。2采样位数采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。采样频率是指录音设备在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。3采样定理在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中，最高频率fmax的2倍时，即：fs.max=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的510倍；采样定理又称奈奎斯特定理。二、语音信号去噪实现框图本次课程设计先完成语音信号的

14、采集，并对所采集的语音信号加入不同的干扰噪声，对加入噪声的信号进行频谱分析，针对受干扰语音信号的特点设计不同的滤波器，然后利用窗函数法设计低通，高通，带通等滤波器对采集到的语音信号进行滤波处理，分析语音信号各频率段的特性。对加噪信号进行滤波，恢复原信号。把原始语音信号、加噪语音信号和滤波后的信号进行时域变换和频域变换，画出它们的时域波形和频域波形图，从视觉角度比较分析滤波的效果。 图2-1 整体设计流程图三、语音信号去噪的详细设计3.1 语音信号的采集y,fs=audioread(C:UsersAdministratorDesktopqwe.wav);%打开音频文件（格式为wav的音频文件），

15、所得y为采样数据，fs为采样率y=y(:,1);y1=fft(y,2048); %对信号做1024点FFT变换f=fs*1:1024/2048;figure(1)subplot(211);plot(y); %绘制原始语音信号的时域波形图title(原始语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f(1:1024),abs(y1(1:1024);title(原始语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on;利用audioread函数

16、对语音信号进行采样，采集出原始信号波形与频谱，y,fs=audioread(C:UsersAdministratorDesktopqwe.wav)，用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)。播放语音信号可以调用函数sound（x,fs）。原始语音信号的时域波形和频谱如图3-1所示：图3-1 原始信号时域波形和频谱3.2 加噪语音信号的频谱分析% 加噪Au=0.3; t=0:1/8000:(size(y)-1)/8000; d=Au*cos(2*pi*2000*t); %加噪2K HZ m,n=size(y);%查看y的大小，【此处y是m行，n列的数据】m1,n1=size(d

17、);%查看y6的大小，【此处y6是m6行，n6列的数据】z=zeros(max(m,m1)-min(m,m1),n);%生成0矩阵，用于加在时间较短的那么音频的后面if length(y)length(d);y1=y;z;y2=y1+d;% sound(y8,fs)else y1=d;z;y2=y1+y; % sound(y8,fs)end; % sound(y2,fs)y3=fft(y2,2048); %对信号做1024点FFT变换f3=fs*1:1024/2048;figure(2)subplot(211);plot(y2); %绘制原始语音信号的时域波形图title(加噪语音信号时域波形

18、图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f3(1:1024),abs(y3(1:1024);title(加噪语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on; % sound(y2,fs)函数中的d是加的高频噪声，语句y2=y1+d;实现了两个信号的相加，然后绘制加噪后的语音信号时域波形和频谱图并回放加噪后的语音信号如图3-2。图3-2加噪信号时域波形图和频谱图3.3 语音信号的滤波去噪3.3.1 FIR数字滤波器的滤波效果 % FIR低通滤波

19、fp1=1000;fs1=3400;As1=100;wp1=2*pi*fp1/fs; %ws1=2*pi*fs1/fs; %BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;%按凯泽窗计算滤波器阶数N1=M1+1;beta1=0.1102*(As1-8.7);Window=(kaiser(N1,beta1); %求凯泽窗窗函数b1=fir1(M1,wc1/pi,Window);% wc1/pi为归一化,窗函数法设计函数 figure(3);freqz(b1,1,512);title(FIR低通滤波器的频率响应);x1_low

20、 = filter(b1,1, y2);%对信号进行低通滤波sound(x1_low,fs); x3=fft(x1_low,2048); %对信号做1024点FFT变换f6=fs*1:1024/2048;figure(4);subplot(211);plot(x1_low); %绘制原始语音信号的时域波形图title(经过FIR低通滤波后语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f6(1:1024),abs(x3(1:1024);title(经过FIR低通滤波后语

21、音信号频谱图)xlabel(Hz);ylabel(fudu); grid on;图3-3凯泽窗频率响应图图3-4FIR滤波后时域波形图与频谱图3.3.2 IIR数字滤波器的滤波效果 % 经过IIR低通滤波fp1i=3000;fs1i=3500;wp1i=2*pi*fp1i/fs;ws1i=2*pi*fs1i/fs;Rp1i=5;Rs1i=50;Ts=1/fs;Wp1i=2/Ts*tan(wp1i/2);Ws1i=2/Ts*tan(ws1i/2); %按频率转换公式进行转换，预畸变N1i,Wn1i=cheb1ord(Wp1i,Ws1i,Rp1i,Rs1i,s); %计算模拟滤波器的最小阶数B1i

22、,A1i=cheby1(N1i,Rp1i,Wn1i,s);%设计模拟原型滤波器bz1i,az1i=bilinear(B1i,A1i,fs); %运用双线性变换法得到数字滤波器传递函数figure(5);freqz(bz1i,az1i,1024,fs);title(切比雪夫1型低通滤波器的频率响应);x1= filter(bz1i,az1i, y2);sound(x1,fs); x2=fft(x1,2048); %对信号做1024点FFT变换f7=fs*1:1024/2048;figure(6);subplot(211);plot(x1); %绘制原始语音信号的时域波形图title(经过IIR低

23、通滤波后语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f7(1:1024),abs(x2(1:1024);title(经过IIR低通滤波后语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on;图3-5切比雪夫1型频率响应图图3-6 IIR滤波后时域波形图与频谱图 由所得结果可知，所设计的滤波器符合要求。 总结此次的语音信号的数字滤波处理课程设计，使我更进一步地熟悉了Matlab软件的使用与滤波器的设计方法。在这次课程设计中，我不仅学到了Ma

24、tlab的软件知识，更增加了个人发现问题解决问题的能力，可谓获益匪浅。通过这次毕业设计，我对语音信号的滤波功能有了全面的认识，对数字信号处理的知识点有了更深层次的理解，进一步了解信号的产生、频谱分析的方法，学会了分析滤波器的优劣和性能，提高了分析问题和解决问题的能力。附录附录一 y,fs=audioread(C:UsersAdministratorDesktopqwe.wav);%打开音频文件（格式为wav的音频文件），所得y为采样数据，fs为采样率y=y(:,1);y1=fft(y,2048); %对信号做1024点FFT变换f=fs*1:1024/2048;figure(1)subplot

25、(211);plot(y); %绘制原始语音信号的时域波形图title(原始语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f(1:1024),abs(y1(1:1024);title(原始语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on; 附录二% 加噪Au=0.3; t=0:1/8000:(size(y)-1)/8000; d=Au*cos(2*pi*2000*t); %加噪2K HZ m,n=size(y);%查看y的大小，【此处y

26、是m行，n列的数据】m1,n1=size(d);%查看y6的大小，【此处y6是m6行，n6列的数据】z=zeros(max(m,m1)-min(m,m1),n);%生成0矩阵，用于加在时间较短的那么音频的后面if length(y)length(d);y1=y;z;y2=y1+d;% sound(y8,fs)else y1=d;z;y2=y1+y; % sound(y8,fs)end; % sound(y2,fs)y3=fft(y2,2048); %对信号做1024点FFT变换f3=fs*1:1024/2048;figure(2)subplot(211);plot(y2); %绘制原始语音信号

27、的时域波形图title(加噪语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f3(1:1024),abs(y3(1:1024);title(加噪语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on; % sound(y2,fs) 附录三 % 经过IIR低通滤波fp1i=3000;fs1i=3500;wp1i=2*pi*fp1i/fs;ws1i=2*pi*fs1i/fs;Rp1i=5;Rs1i=50;Ts=1/fs;Wp1i=2/Ts*tan(

28、wp1i/2);Ws1i=2/Ts*tan(ws1i/2); %按频率转换公式进行转换，预畸变N1i,Wn1i=cheb1ord(Wp1i,Ws1i,Rp1i,Rs1i,s); %计算模拟滤波器的最小阶数B1i,A1i=cheby1(N1i,Rp1i,Wn1i,s);%设计模拟原型滤波器bz1i,az1i=bilinear(B1i,A1i,fs); %运用双线性变换法得到数字滤波器传递函数figure(5);freqz(bz1i,az1i,1024,fs);title(切比雪夫1型低通滤波器的频率响应);x1= filter(bz1i,az1i, y2);sound(x1,fs); x2=ff

29、t(x1,2048); %对信号做1024点FFT变换f7=fs*1:1024/2048;figure(6);subplot(211);plot(x1); %绘制原始语音信号的时域波形图title(经过IIR低通滤波后语音信号时域波形图);xlabel(time n);ylabel(fuzhi n); grid on; subplot(212); %绘制原始语音信号的频率响应图 plot(f7(1:1024),abs(x2(1:1024);title(经过IIR低通滤波后语音信号频谱图)xlabel(Hz);ylabel(fudu); grid on; 附录四 % FIR低通滤波 fp1=10

30、00;fs1=3400;As1=100;wp1=2*pi*fp1/fs; %ws1=2*pi*fs1/fs; %BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;%按凯泽窗计算滤波器阶数N1=M1+1;beta1=0.1102*(As1-8.7);Window=(kaiser(N1,beta1); %求凯泽窗窗函数b1=fir1(M1,wc1/pi,Window);% wc1/pi为归一化,窗函数法设计函数 figure(3);freqz(b1,1,512);title(FIR低通滤波器的频率响应);x1_low = filter(b1,1,