语音信号滤波去噪——使用双线性变换法设计的切比雪夫II型滤波器
一、引言 (3)
课程设计目的 (3)
课程设计要求 (3)
二、(
(4)
三、设计原理
IIR滤波器 (4)
切比雪夫I型滤波器 (4)
双线性变换法 (7)
四、设计步骤 (8)
设计流程图 (8)
语音信号的采集 (9)
语音信号的频谱分析 (10)
？
滤波器设计 (12)
完整的滤波程序及滤波效果图 (15)
结果分析 (18)
五、出现的问题及解决方法 (20)
六、课程设计心得体会 (21)
七、参考文献 (22)
【
语音信号滤波去噪
——使用双线性变换法设计的切比雪夫I型滤波器
摘要随着信息技术的发展，现代信号处理正向着数字化、软件化方向发展。滤波器设计是信号处理的重要组成部分，而研究语音信号的滤波设计是现代信息处理的基本内容。本设计利用计算机WINDOWS下的录音机录入一句语音信号，用MATLAB软件对其进行频谱分析，然后加入一干扰信号，利用设计好的滤波器将干扰信号去除，最后对各部分的频谱进行分析比较。
关键词滤波设计； MATLAB；
一、引言
用麦克风采集一段8000Hz，8k的单声道语音信号，绘制波形并观察其频谱，给定通带截止频率为2000Hz，阻带截止频率为2100Hz，通带波纹为1dB，阻带波纹为60dB，用双线性变换法设计的一个满足上述指标的切比雪夫II型IIR滤波器，对该语音信号进行滤波去噪处理。
*
课程设计目的
《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。
课程设计的要求
（1）学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法；
（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；（3）采用双线性变换法，设计满足上述性能指标要求的ChebyshevII型数字低通滤波
器；
@
（4）设计完后应检查其频率响应曲线是否满足指标；
（5）处理结果和分析结论应该一致，而且应符合理论；
（6）独立完成课程设计并按要求编写课程设计报告书；
二、设计原理
用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用脉
冲响应不变法设计的一个满足指标的巴特沃斯IIR滤波器，对该语音信号进行滤波去
噪处理，比较滤波前后的波形和频谱并进行分析。
IIR滤波器
【
从离散时间来看，若系统的单位抽样(冲激)响应延伸到无穷长，称之为“无限长
单位冲激响应系统”，简称为IIR系统。
无限长单位冲激响应（IIR）滤波器有以下几个特点:
(1)系统的单位冲激响应h(n)是无限长；
(2)系统函数H(z)在有限z平面（0 (3)结构上存在着输出到输入的反馈，也就是结构上是递归型的。 IIR滤波器采用递归型结构，即结构上带有反馈环路。同一种系统函数H（z） 可以有多种不同的结构，基本网络结构有直接Ⅰ型、直接Ⅱ型、级联型、并联型 四种，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤 波器的成果，巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭 圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些典型的滤波器各有特点。有 现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 > 切比雪夫I型滤器 2.2.1 切比雪夫滤波器简介 切比雪夫滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。 这种滤波器来自切比雪夫多项式，因此得名，用以纪念俄罗斯数学家巴夫尼提·列波维其·切比雪夫 2.2.2切比雪夫滤波器原理 巴特沃兹滤波器在通带内幅度特性是单调下降的，如果阶次一定，则在靠近截止处，幅度下降很多，或者说，为了使通带内的衰减足够小，需要的阶次N很高，为了克服这一缺点，采用切比雪夫多项式来逼近所希望的。切比雪夫滤波器的在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。 切比雪夫滤波器的振幅平方函数为 ! （1） 式中 Ωc —有效通带截止频率 —与通带波纹有关的参量， 大,波纹大 0< <1 V N （x ）—N 阶切比雪夫多项式 （2） |x|≤1时,|V N (x)|≤1 】 |x|>1时, |x|↗, V N (x)↗ 切比雪夫滤波器的振幅平方特性如图所示,通带内， 的变化范围为 1(max) → 2 11 ε+ (min) 时，|x|>1,随 ↗， →0 (迅速趋于零) 当 =0时， （3） N 为偶数，cos 2( )=1，得到min ， 、 ， （4） N 为奇数，cos 2( ，得到max ， （5） 切比雪夫滤波器的振幅平方特性如图1所示。 % 图1 切比雪夫滤波器的振幅平方特性 双线性变换法 双线性变换法是使数字信号滤波器的频率响应与模拟滤波器的频率响应相似的一 种变换方法。为了客服多值映射这一缺点，我们首先把整个s 平面压缩变换到某一中介的s1平面的一条横带里（宽度为 T π2,即从T π-到T π ),其次再通过上面讨论过的标准变换关系T s e z 1=将此横带变换到这个z 平面上去，这样就使s 平面与z 平面式一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。 将s 平面整个Ωj 轴压缩变换到s1平面Ωj 轴上的T π-到T π 一段，可以采用以下变换关系： )2 tan( 1T Ω=Ω （6） - 这样，±∞=Ω变换到T π ± =Ω1，0=Ω变到01=Ω可将（6）式写成 2 2 221111T j T j T j T j e e e e j Ω-ΩΩ-Ω--= Ω （7） 解析延拓到整个s 平面和s1平面，令s j =Ω，11s j =Ω，则得 T s T s T j T j T j T j e e T s th e e e e j 11 111111212 2 22--Ω-ΩΩ-Ω+-=?? ????=--= Ω （8） 再将s1平面通过以下标准变换关系映射到z 平面： T s e z 1= （9） 从而得到s 平面和z 平面的单值映射的关系为 1 1 11--+-=z z s （10） ) s s z -+= 11 （11） 一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应的关系，可以引入待定常数c ，使（6）式和（7）式变换成 )2 tan( 1T c Ω=Ω （12） T s T s e e c T s cth s 111121--+-=??? ???= （13） 仍将T s e z 1=代入（13）式，可得 1 1 11--+-=z z c s （14） s c s c z -+= （15） （14）式和（15）式是s 平面与z 平面之间的单值映射关系，这种变换就称为双线性变换。 . 三、设计步骤 设计流程图 语音信号滤波去噪——使用双线性变换法设计的切比雪夫I 型滤波器的设计流程如图3.1.1所示： 【 . 双线性变换法切比雪夫II型滤波器对语音信号去噪流程图 语言信号的采集 点击windows系统桌面的“开始”按钮，点击开始菜单栏里的“附件”，选择“录音机”选项，点击录音机“文件”选项，进入“声音选定”设置，把属性一栏设置成“8000Hz，8位，单声道，7KB/秒”(见图3.2.1)。点击确定，然后开始语言信号的采集，采集时间为1秒左右为最佳。采集的声音文件以“.wav”格式存储（见图）。 % 图3.2.1 采集声音的参数设置 图3.2.2 采集声音 、 语音信号的频谱分析 在MATLAB中编辑m函数，使用wavread函数读取采集的声音文件（.wav)将它赋值给某一向量，再对其进行采样,然后使用plot语句画出相关的频谱图形在 figure(1）上。 （1）Wavread函数调用格式： [y,Fs,nbits]=wavread(file) ^ 功能说明：采样值放在向量y中，Fs表示采样频率（Hz），nbits表示采样位数。 （2）快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下： Xk=fft(x,n) 参数x为被变换的时域序列向量，N是DFT变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。，当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。 在本次课程设计中，我们利用fft函数对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。 （3）声音采样文件读取的程序（文件名： %用麦克风采集一段8000Hz，8k的单声道语音信号，绘制波形并观察其频谱 Fs=8000; %语音信号采样频率为8000 | x1=wavread(''); %读取语音信号的数据，赋给变量x1 sound(x1,8000); %播放语音信号 y1=fft(x1,8000); %对信号做FFT变换 subplot(1,2,1); plot(x1) %做原始语音信号的时域图形 title('原始语音信号时域图'); xlabel('时间 n'); ylabel('音量 n'); - subplot(1,2,2); plot(abs(y1(1:8000))) %做原始语音信号的FFT频谱图 axis([0,4000,0,10]) title('原始语音信号FFT频谱') figure freqz(x1) %绘制原始语音信号的频率响应图 title('原始语音信号频率响应图') 所得图形 ` 图3.3.1 原始语音信号的时域图和FFT频谱 图3.3.2原始语音信号频率响应图 滤波器设计 设计指标：通带截止频率为1100Hz，阻带截止频率为1200Hz，通带波纹为1dB，阻带波纹为20dB，用脉冲（冲激）响应不变法设计的一个满足上述指标的切比雪夫I 型滤波器 双线性变换法设计切比雪夫I型滤波器 Fs=8000; %采样频率 wp=1100*2/Fs;ws=1200*2/Fs;%根据采样频率将滤波器边界进行转化 Rp=1; %通带波纹 — Rs=20;%阻带波纹 滤波程序如下： fp=2000;fs=2100; %给出数字滤波器的参数 Rp=1; As=60; wp=fp/Fs*2*pi; %将模拟指标转换成数字指标 ws=fs/Fs*2*pi; } Ts=1/Fs; wp1=2/Ts*tan(wp/2); %频率预畸 ws1=2/Ts*tan(ws/2); [N,Wn]=cheb1ord(wp1,ws1,Rp,As,'s'); %选择滤波器的最小阶数 %创建chebyshef模拟滤波器 [b,a]=cheby1(N,1,Wn,'low','s'); [bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 （ [H,W]=freqz(bz,az); %绘制数字滤波器频率响应幅度图 figure; plot(W*Fs/(2*pi),abs(H)) axis([0,5000,0,1]) grid xlabel('频率／Hz') ylabel('滤波器频率响应幅度') 、 title('切比雪夫型') figure; freqz(bz,az,1024,Fs); title('数字滤波器的频率响应图') xlabel('Hz'); ylabel('fuzhi'); sound(f1); %播放滤波之后的信号… figure(9) plot(angle(f1(1:8000))) 所得图形有： 图3.4.1 滤波器频率响应幅度（切比雪夫型） % 图3.4.2 数字滤波器的频率响应图 图3.4.3 数字滤波器的相角曲线图 完整的滤波程序及滤波效果图%%%%%%%%%%%%%%%%%%%%% 读入声音，对其进行FFT变换 %%%%%%%%%%%%%%%%%%%%%% Fs=8000; %语音信号采样频率为8000 … x1=wavread(''); %读取语音信号的数据，赋给变量x1 sound(x1,8000); %播放语音信号 y1=fft(x1,8000); %对信号做FFT变换 subplot(1,2,1); plot(x1) %做原始语音信号的时域图形 title('原始语音信号时域图'); xlabel('时间 n'); ylabel('音量 n'); ? subplot(1,2,2); plot(abs(y1(1:8000))) %做原始语音信号的FFT频谱图 axis([0,4000,0,10]) title('原始语音信号FFT频谱') figure freqz(x1) %绘制原始语音信号的频率响应图 title('原始语音信号频率响应图') 》 fp=2000;fs=2100; %给出数字滤波器的参数 Rp=1; As=60; wp=fp/Fs*2*pi; %将模拟指标转换成数字指标 ws=fs/Fs*2*pi; $ Ts=1/Fs; wp1=2/Ts*tan(wp/2); %频率预畸 ws1=2/Ts*tan(ws/2); [N,Wn]=cheb1ord(wp1,ws1,Rp,As,'s'); %选择滤波器的最小阶数 %创建chebyshef模拟滤波器 [b,a]=cheby1(N,1,Wn,'low','s'); ( [bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,W]=freqz(bz,az); %绘制数字滤波器频率响应幅度图figure; plot(W*Fs/(2*pi),abs(H)) axis([0,5000,0,1]) grid xlabel('频率／Hz') ' ylabel('滤波器频率响应幅度') title('切比雪夫型') figure; freqz(bz,az,1024,Fs); title('数字滤波器的频率响应图') 、 %开始滤波 f1=filter(bz,az,x1); figure subplot(2,2,1) plot(x1) %画出滤波前的时域图 title('滤波前的时域波形'); subplot(2,2,2) [ plot(f1); %画出滤波后的时域图 title('滤波后的时域波形'); f=fs*(0:7999)/8000; subplot(2,2,3); plot(f,abs(y1(1:8000))); %画出滤波前的频谱图axis([0,4000,0,10]) ' title('滤波前的频谱') xlabel('Hz'); ylabel('fuzhi'); F0=fft(f1,8000); subplot(2,2,4) F1=plot(f,abs(F0(1:8000))); %画出滤波后的频谱图> axis([0,4000,0,10]) title('滤波后的频谱') xlabel('Hz'); ylabel('fuzhi'); sound(f1); %播放滤波后的信号xlabel('Hz'); ylabel('fuzhi'); sound(f1); %播放滤波之后的信号 。 figure(9) plot(angle(f1(1:8000))) 所得图形： 图3.5.1 滤波前后的图形 结果分析 ) 在MATLAB中，经sound函数，对经过切比雪夫I型滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。 通过以下语句来进行语音信号回放比较： >>sound（x，fs，bits） >>sound（z，fs，bits） 所得结果证明了切比雪夫I型滤波器去噪设计成功。 } 四、出现的问题及解决方法 在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。 （1）在进行语音信号提取时，进过多次录取才得到理想的语音信号，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的语音信号的时域图和频谱图 （2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的语音信号wav文件放到Matlab的work文件夹中。 （3）还要在滤波器性能曲线的wc处画一根竖线，这样更方便看出结果，其中wc处线的确定还需计算出wb/pi的值。 （4）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。