1、第 1 页语音信号滤波去噪使用 TUKEYWIN窗设计的 FIR 滤波器摘 要 本课程设计主要使用 TUKEYWIN 窗设计的 FIR 滤波器对语音信号进行滤波去噪。课程设计的的平台为 MATLAB。采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用窗函数法设计一个满足指标的 FIR 滤波器，对该语音信号进行滤波去噪处理。根据滤波后的时域图和原始语音信号时域图的比较，以及滤波后信号的频谱图和原始音乐信号频谱图的比较，最后回放滤波后音乐信号，滤波后的音乐信号与原始音乐信号一样清晰，成功地实现了滤波，达到了设计的要求。关键词 课程设计；滤波去噪；FIR 滤波器；TUKEYWIN 窗；MA
2、TLAB1 引 言本课程设计是采用 TUKEYWIN 窗设计的 FIR 滤波器对语音信号进行滤波去噪。通过课程设计了解 FIR 滤波器设计的原理和步骤，掌握用 Matlab 语言设计滤波器的方法，了解 DSP 对 FIR 滤波器的设计及编程方法。通过观察音乐信号滤波前后的时域波形的比较，加深对滤波器作用的理解。通过对比滤波前后波形图的比较和放滤波前后音乐信号的对比 ，可以看出滤波器对有用信号无失真放大具有重大意义。1.1 课程设计目的熟悉 Matlab 语言环境，掌握 Matlab 语言的编程规则，利用 TUKEYWIN 窗函数设计法来设计符合要求的 FIR 滤波器来实现语音信号的滤波去噪。并
3、绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用 TUKEYWIN 窗函数法设计 FIR 滤波器的过程。增强自己独立解决问题的能力，提高第 2 页自己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。1.2 课程设计要求录制一段语音信号，绘制观察波形及频谱图。根据 TUKEYWIN 的性能指标合理设计 FIR 滤波器，对语音信号加入干扰，再用滤波器对干扰语音信号进行滤波去噪，比较滤波前后的频谱图并进行分析。再回放语音信号对比原语音信号。看滤波器是否对语音信号进行了滤波去噪。1.3 课程设计平台MATLAB 名字由 MATrix
4、 和 LABoratory 两词的前 3 个字母组合而成。20 世纪 70年代后期，时任美国新墨西哥大学计算机科学系主任的 Cleve Moler 教授出于减轻学生编程负担的动机，为学生设计了一组调用 LINPACK 和 EISPACK 库程序的“通俗易用”的接口，此即用 FORTRAN 编写的萌芽状态的 MATLAB1。MATLAB 软件包括五大通用功能：数值计算功能（ Nemeric）；符号运算功能（Symbolic）；数据可视化功能（Graphic）；数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。其中，符号运算功能的实现是通过请求 MAPLE 内核
5、计算并将结果返回到 MATLAB 命令窗口。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30 多个应用工具箱。MATLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。MATLAB 在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而 MATLAB 借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求
6、。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z 正反变换等。MATLAB 在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的 S 域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用第 3 页MATLAB 进行信号处理领域的各种分析和实际应用打下基础 2。2 基本理论2.1 FIR 滤波器FIR 滤波器：有限长单位冲激响应滤波器，是数字信号系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的
7、，因而滤波器是稳定的系统。因此，FIR 滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR 数字滤波器设计的基本步骤如下：（1）确定指标在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对指标。它提供对幅度响应函数的要求，一般应用于 FIR 滤波器的设计。第二种指标是相对指标。它以分贝值的形式给出要求。（2）逼近确定了技术指标后，就可以建立一个目标的数字滤波器模型。通常用理想的数字滤波器模型。之后，利用数字滤波器的设计方法，设
8、计出一个实际滤波器模型来逼近给定的目标。（3）性能分析和计算机仿真上两步的结果是得到以差分或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。FIR 滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。2.2 窗口设计法窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器第 4 页的技术指标。在大多数实际应用中，数字滤波器常
9、常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。窗口设计法步骤如下：（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度 N。窗函数的类型可根据最小阻带衰减 AS 独立选择。（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应 hd(n)。（3）由性能指标确定窗函数 W(n)和长度 N。（4）求得实际滤波器的单位脉冲响应 h(n)， h(n)即为所设计 FIR 滤波器系数向量b(n)。()().dhnWn (2.1)常见的窗函数性能表如下图1所示：图1 常见窗函数性能表名称 滤波器过渡带宽 最小阻带 衰减 名称 滤波器过渡带宽 最小阻带衰
10、减矩形 1.8/M 21dB PARZENWIN 6.6/M 56dB巴特利特 6.1/M 25dB FLATTOPWIN 19.6/M 108dB汉宁 6.2/M 44dB GAUSSWIN 5.8/M 60dB汉明 6.6/M 51dB BARTHANNWIN 3.6/M 40dB布莱克曼 11/M 74dB BLACKMANHARRIS 16.1/M 109dBBOHMANWIN 5.8/M 51.5dB CHEBWIN 15.2/M 113dBNUTTALLWIN 15.4/M 108dB TUKEYWIN 2.4/M 22dB第 5 页2.3 TUKEYWIN 窗W=tukeywin
11、(M)在数组 w 中产生 M 点的 Tukeywin 窗函数。由上表可知：Tukeywin 窗的滤波器的过滤带宽为 2.4/M，最小阻带衰减 22db。利用 wvtool 可绘制出 Tukeywin 窗函数波形如图 2：N=128;t=tukeywin(N);Wvtool(t);图 2 Tukeywin 窗函数波形3 设计步骤3.1 设计流程图根据设计的要求，首先自己录制一段语音信号，对语音信号加入噪声干扰，再利用 TUKEYWIN 窗设计合理的 FIR 滤波器。最后用滤波器对干扰后的语音信号进行滤波去噪。具体设计流程图如下图 3 所示：开始录制段语音信号，命名为“speech_dft.wav
12、”，绘制出其时域波形和频谱图。第 6 页NOYES图 3 设计流程图 3.2 录制音乐信号从电脑上录制一段语音信号，并命名为“speech_dft.wav” 。在 Matlab 平台下调用wavread 函数可采集到语音信号。运行程序后得到的波形图如图 4 所示：加入噪声干扰，绘制干扰后的时域波形和频谱图。利用 TUKEYWIN 窗设计合理的 FIR 滤波器。将干扰后的语音信号通过自己设计的滤波器，进行滤波去噪。比较滤波前后的时域波形和频谱图，并回放语音信号，验证是否达到去噪效果。结束是否达到去噪 效果？第 7 页图 4 原始音乐信号时域波形图画出语音信号的时域波形；然后对语音号进行快速傅里叶
13、变换，得到信号的频谱特性，并将原始音乐信号的波形图与加干扰后的波形图进行比较运行程序后我们所得的波形图如图 5 所示：图 5 干扰前后信号时域图与频谱图的比较3.3 滤波器设计滤波器的设计就是要找到一组能满足特定滤波要求的系数向量 a 和 b，其中 b 为系统函数的分子系数，a 为系统函数分母系数。在本次的课程设计中我所采用的就是利用TUKEYWIN 窗函数来设计 FIR 滤波器如图 6 所示。第 8 页图 6 设计的滤波器图形3.4 信号滤波处理滤波器设计完成后，在 MATLAB 平台上用函数 filter 实现滤波， 得到的滤波前后音乐信号的时域波形图和频谱图对比图如图 7 所示：第 9
14、页图 7 滤波前后语音信号的时域波形图和频谱图对比图3.5 结果分析在 MATLAB 中，经过 sound(x,fs,bits)函数，对经过 TUKEYWIN 窗设计的 FIR 滤波器之后的音乐信号进行回放，可以听出滤波之后的信号跟原始信号一样清晰，完全滤除掉了噪声的干扰，通过下面语句来进行语音信号回放比较： sound (y_fil,fs,bits); 播放滤波后的音乐信号 sound(x,fs,bits);所得结果证明了用 TUKEYWIN 窗设计的 FIR 滤波器和音乐信号去噪设计是成功的。4 出现的问题及解决方法在本次课程设计中我遇到的问题如下：第 10 页1、不知道如何将音乐信号截断
15、成自己所需的长度。2、对利用 TUKEYWIN 窗函数设计 FIR 滤波器的设计步骤很生疏。3、在采用 TUKEYWIN 窗函数设计的 FIR 滤波器时得不到理想的滤波器，因而信号的恢复不是特别理想。针对以上问题，相应的解决方案如下：1、向同学学习如何截取音乐。2、自己上网查阅资料，或向图书馆借阅资料，掌握利用 TUKEYWIN 窗函数设计 FIR滤波器的方法和步骤。3、通过不断设置参数的值，最终达到最理想的值，设计出理想的滤波器，使信号得到理想恢复。5 结束语本次的课程设计，我的任务是利用 TUKEYWIN 窗函数设计 FIR 滤波器对音乐信号滤波去噪。在本次课程设计之前，我对 TUKEYW
16、IN 窗函数完全没有了解，因此在看到这个题目时，我是一头雾水。但是通过自己翻阅资料和询问同学，我掌握了用TUKEYWIN 窗函数设计 FIR 滤波器的方法步骤，了解了窗函数的基本设计流程。经过这两周的课程设计，我学会了很多东西。我们通信工程专业是个实践性很强的专业，而我们在校大部分的学习时间都是花在理论学习上面，实践的机会很少。因而我对很多所学的理论知识如何跟实践联系的概念很模糊，这次的课程设计给了我这个机会，加深了我对理论联系实际的理解，增强了自己独立分析问题和解决问题的能力，开阔了自己的思维。还有让我看到了自己的不足，自己对本专业的相关知识掌握的还很少，还有很多知识都没掌握，还让我认识到解
17、决问题的方法、途径很多，做事要开阔自己的思维，看待问题要从多个角度看。在此我要感谢学校为我们提供这次课程设计的机会，感谢老师对我的悉心指导，也感谢同学对我的帮助。这次的课程设计让我理论联系实际，不仅巩固了我们的理论知识，还提高了我的动手能力，在这次课程设计中我所学到的知识是我的财富，让我终身受益。第 11 页参考文献1张志涌精通 MATLAB 6.5 版M北京：北京航空航天大学出版社，20032约翰 普罗克斯数字信号处理M西安：西安交通大学出版社，20093张小虹信号系统与数字信号处理M第 1 版西安：西安电子科技出版社，20024谢德芳数字信号处理M北京：科学出版社，20055郝建军数字通信
18、M 第 2 版北京：北京邮电大学出版社，20106张威MATLAB 基础与编程入门 M西安：西安电子科技大学出版社，2010第 12 页附录一： 语音信号滤波去噪使用 TUKEYWIN 滤波器%程序功能：在 Matlab 中，用窗口设计法设计 FIR 滤波器%程序作者： %最后修改日期：%程序一：x,fs,bits=wavread(speech_dftwav);%读取音乐信号sound(x,fs,bits);%播放音乐信号N=length(x);%计算音乐信号的长度fn=2136; % 单频噪声频率t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); %
19、 将双声道转为单声道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); % 截取前半部分figure(1)subplot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title( 原始语音信号);axis(0,2.5,-0.7,0.7);subplot(2,2,2);plot(f,X);xlabel(频率
20、(f) );ylabel(幅度谱 );title(原始语音信号幅度谱);axis(0,4000,0,300);subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title( 加噪声后的语音信号);axis(0,2.5,-0.7,0.7);subplot(2,2,4);plot(f,Y);xlabel(频率(f) );ylabel(幅度谱 );title(加噪声后的语音信号幅度谱);axis(0,4000,0,300);%程序三：fpd=2000;fsd=2050;fsu=2150;fpu=2200;% 带阻滤波器设计指标Rp=1;As=22; %
21、 通带波纹和阻带衰减 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率，和频率间隔第 13 页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(2.4*pi/dw)+1; % 计算 Tukeywin 窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_tuk=tukeywin(M); % 产生 M 阶的 Tuk
22、eywin 窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_tuk.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性figure(2)subplot(2,2,1);plot(w/pi,db);xlabel(w/pi);ylabel(db);title(FIR 滤波器的幅度响应图 );axis(0.3 0.7 -70 10);line(0,0.7,-1,-1,Color,
23、r,LineWidth,2,LineStyle,-);line(0,0.7,-22,-22,Color,r,LineWidth,2,LineStyle,-);line(wsd/pi,wsd/pi,-70,10,Color,r,LineWidth,2,LineStyle,-);line(wsu/pi,wsu/pi,-70,10,Color,r,LineWidth,2,LineStyle,-);grid on;subplot(2,2,2);plot(w/pi,mag);xlabel(w/pi);ylabel(幅度 mag);title(FIR 滤波器的幅度响应图);axis(0 1 -0.2 1.
24、2);grid on;subplot(2,2,3);plot(w/pi,pha);xlabel(w/pi);ylabel(相位 pha);title(滤波器相位响应图);axis(0 1 -4 4);grid on;subplot(2,2,4);stem(n,h_bs);xlabel(n);ylabel(h(n);title(滤波器脉冲响应图 );grid on;%程序四：y_fil=fftfilt(h_bs,y); % 用设计好的滤波器对 y 进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半sound (y_fil,fs,bits
25、);wavwrite(y_fil,fs,bits,1111 恢复.wav ); %保存恢复的语音信号wavwrite(y,fs,bits,1111 加噪.wav ); %保存加噪的语音信号figure(3)subplot(3,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title( 原始语音信号);axis(0,2.5,-0.7,0.7);subplot(3,2,2);plot(f,X);xlabel(频率(f) );ylabel(幅度谱 );title(原始语音信号幅度谱);axis(0,4000,0,300);subplot(3,2,3);plot(t,
26、y);xlabel(时间(t);ylabel(幅度);title( 加噪声后的语音信号);axis(0,2.5,-0.7,0.7);subplot(3,2,4);plot(f,Y);xlabel(频率(f) );ylabel(幅度谱 );title(加噪声后的语音信号幅度谱);axis(0,4000,0,300);subplot(3,2,5);plot(t,y_fil);xlabel(时间(t);ylabel(幅度); title(滤波后语音信号);axis(0,2.5,-0.7,0.7);subplot(3,2,6);plot(f,Y_fil);xlabel(频率(f) );ylabel(幅
27、度谱 );title(滤波后语音信号幅度谱);axis(0,4000,0,300);第 14 页附录二：函数 FREQZ_M.M 定义：function db,mag,pha,grd,w = freqz_m(b,a);% freqz 子程序的改进版本% -% db,mag,pha,grd,w = freqz_m(b,a);% db = 0 到 pi 弧度 区间内的相对振幅(db)% mag = 0 到 pi 弧度区间内的绝对振幅% pha = 0 到 pi 弧度区间内的相位响应% grd = 0 到 pi 弧度区间内的群迟延% w = 0 到 pi 弧度 区间内的 501 个频率样本向量% b
28、 = Ha(z)的分子多项式系数(对 FIR b=h)% a = Ha(z)的分母多项式系数(对 FIR: a=1)%H,w = freqz(b,a,1000,whole);H = (H(1:1:501); w = (w(1:1:501);mag = abs(H);db = 20*log10(mag+eps)/max(mag);pha = angle(H);% pha = unwrap(angle(H);grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0 grd(1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;附录三：第 15 页函数 IDEAL_LP.M 定义：function hd = ideal_lp(wc,M);% 理想低通滤波器计算% -% hd = ideal_lp(wc,M)% hd = 0 to M-1 之间的理想脉冲响应% wc = 截止频率(弧度) % M = 理想滤波器的长度%alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha + eps;hd = sin(wc*m) ./ (pi*m);
