语音信号滤波去噪使用双线性变换法设计的级联型椭圆滤波器.docx
- 文档编号:18087259
- 上传时间:2023-08-13
- 格式:DOCX
- 页数:30
- 大小:534.50KB
语音信号滤波去噪使用双线性变换法设计的级联型椭圆滤波器.docx
《语音信号滤波去噪使用双线性变换法设计的级联型椭圆滤波器.docx》由会员分享,可在线阅读,更多相关《语音信号滤波去噪使用双线性变换法设计的级联型椭圆滤波器.docx(30页珍藏版)》请在冰点文库上搜索。
语音信号滤波去噪使用双线性变换法设计的级联型椭圆滤波器
语音信号滤波去噪
——使用双线性变换法设计的级联型椭圆滤波器
学生姓名:
Su指导老师:
摘要本课程设计主要内容是利用双线性变换法设计一个级联型的椭圆IIR滤波器,对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。
本课程设计仿真平台为MATLAB7.0,开发工具是M语言编程。
首先在windows下用录音机工具录制一段语音信号,并人为加入一单频噪声,然后对信号进行频谱分析以确定所加噪声频率,并设计滤波器进行滤波去噪处理,最后比较滤波前后的波形和频谱并进行分析。
由分析结果可知,滤波器后的语音信号与原始信号基本一致,即设计的IIR滤波器能够去除信号中所加单频噪声,达到了设计目的。
关键词滤波去噪;IIR椭圆滤波器;双线性变换法;级联型;MATLAB
1引言
信号处理是科学研究和工程技术许多领域都需要进行的一个重要环节,传统上对信号的处理大都采用模拟系统实现。
随着人们对信号处理要求的日益提高,以及模拟信号处理中一些不可克服的缺点,对信号的许多处理而采用数字的方法进行。
数字信号处理系统无论在性能、可靠性、体积、耗电量、成本等诸多方面都比模拟信号处理系统优越的多,使得许多以往采用模拟信号处理的系统越来越多地被数字处理系统所代替,数字信号处理技术在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。
滤波器是一种用来消除干扰杂讯的器件,凡是可以使信号中特定的频率成分通过,而极大地衰减或抑制其他频率成分的装置或系统称之为滤波器。
在数字信号处理中,数字滤波器设计在电子工程、应用数学和计算机科学领域都是非常重要的内容。
设计滤波器的方法有多种,在各种滤波器中,椭圆滤波器相比其他类型的滤波器,在阶数相同的条件下有着最小的通带和阻带波动,它在通带和阻带的波动相同。
本课程设计主要解决在含噪情况下对语音信号的滤波去噪处理,处理时采用的是利用双线性变换法设计的级联型的椭圆IIR滤波器。
1.1课程设计目的
本课程设计主要是用录音机采集一段语音信号,利用MATLAB,加入一个带外单频噪声,使用双线性变换法设计的级联型椭圆滤波器,对该含噪语音信号进行滤波去噪处理,画出级联型滤波器结构图。
比较滤波前后的波形和频谱并进行分析,根据结果和学过的理论得出合理的结论,从而加深对所学知识的理解,加强我们的动手操作能力。
另一方面也是对课堂所学的理论知识应用于实践中,熟悉和巩固滤波器的设计原理与方法,掌握双线性变换法级联型椭圆滤波器的设计方法和有关滤波器设计的相关经典算法,熟练地利用MATLAB语言进行设计滤波器。
通过本次课程设计,可以运用所学的东西解决一些实际的问题,提高自身对所学知识的综合运用能力以及分析和解决问题的能力。
1.2课程设计的要求
(1)滤波器指标必须符合工程实际。
(2)设计完后应检查其频率响应曲线是否满足指标。
(3)处理结果和分析结论应该一致,而且应符合理论。
(4)独立完成课程设计并按要求编写课程设计报告书。
1.3设计平台
本次课程设计中所用到的软件平台是是由美国mathworks公司发布的MATLAB软件。
MATLAB是矩阵实验室(Matrix Laboratory)的简称,它和Mathematical、Maple并称为三大数学软件。
除具备卓越的数值计算能力外,它还提供了专业水平的符号计算,文字处理,可视化建模仿真和实时控制等功能。
MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等,主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB的应用范围非常广,包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。
附加的工具箱(单独提供的专用MATLAB函数集)扩展了MATLAB环境,可以解决这些应用领域内特定类型的问题。
现在MATLAB已经成为国际上最流行的科学与工程计算软件工具,它已不仅仅是一个“矩阵实验室”了,它已经成为一种最具有广泛应用前景的全新的计算机高级语言,在世界上很多国家,MATLAB已经成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具;成为攻读学位的大学生、硕士生、博士生必须掌握的基本技能。
在设计研究单位和工业部门,MATLAB也被广泛用于科学研究和解决各种具体问题当中[1]。
2设计原理
2.1IIR滤波器
IIR滤波器,又名“无限脉冲响应滤波器”,或“递归滤波器”。
递归滤波器,也就是IIR数字滤波器,顾名思义,具有反馈性,一般认为具有无限的脉冲响应,因此IIR滤波器能够与模拟滤波器相匹敌[2]。
一个IIR滤波器的系统函数给出为
(2.1)
其中
和
是滤波器系数。
不失一般性已假定
。
如果
。
N就是这个IIR滤波器的阶。
一个IIR滤波器的差分方程表示是[1]
(2.2)
IIR滤波器具有系统的冲击响应
是无穷长的、系统函数
在有限z平面
上和结构上存在递归的特点。
一个IIR滤波器可以用三种结构实现,这三种结构包括直接型、级联型和并联型,其中直接型又可分为直接Ⅰ型和直接Ⅱ型两种。
级联型中,将系统函数H(z)写成具有实系数的二阶节的乘积
(2.3)
式中
,
,
,
和
都是代表实数的二阶节系数。
这些二阶节是
(2.4)
称为第k个双二阶节。
每个级联型结构的滤波器都可以由k个双二阶节用直接Ⅱ型实现,双二阶节结构图如图2-1所示。
图2-1双二阶节结构
一般说来,所有的模拟滤波器都有无限长的脉冲响应,因此IIR滤波器的设计基本方法是利用复值映射将模拟滤波器变换为数字滤波器。
目前,已有许多成熟的模拟滤波器应用于信号处理中,如巴特沃兹滤波器、切比雪夫(Ⅰ型和Ⅱ型)滤波器、椭圆滤波器等,我们将这些各有特色的模拟滤波器成为原型滤波器。
在设计IIR滤波器时,可以直接借用这些典型的模拟滤波器进行转换,这种方法称为A/D(模拟—数字)滤波器变换。
因为AFD表格仅对低通滤波器适用,但在实际应用中也会需要设计其他频率选择性滤波器(高通、带通、带阻滤波器等等)[3]。
为此需要对低通滤波器实行频带变换,基本方法存在两种途径,如图2-2所示。
途径1:
途径2:
图2-2IIR滤波器的两种设计途径
在MATLAB中采用第一种途径设计IIR滤波器。
2.2椭圆低通滤波器
椭圆滤波器又称考尔滤波器,是通带和阻带等波纹的一种滤波器。
在阶数相同的条件下,与巴特沃兹滤波器和切比雪夫滤波器相比,椭圆滤波器有着最小的通带和阻带波动。
就这点而言,椭圆滤波器是最优的。
从传递函数来看,椭圆函数滤波器在有限频率上既有零点又有极点。
极零点在通带内产生等波纹,阻带内的有限传输零点减少了过渡区,可获得极为陡峭的衰减曲线。
但它陡峭的过渡带特性是用通带和阻带的起伏为代价来换取的,并且在通带和阻带的波动相同,这一点区别于在通带和阻带都平坦的巴特沃斯滤波器,以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。
椭圆滤波器的传输函数是一种较复杂的逼近函数,利用传统的设计方法进行电路网络综合要进行繁琐的计算,还要根据计算结果进行查表,整个设计,调整都十分困难和繁琐。
为此,常常需要用一些程序和表格来设计椭圆滤波器,用MATLAB设计它们可以大大简化设计过程。
椭圆滤波器的幅度平方响应给出为
(2.5)
式中
是通带波纹(它与
有关),
是N阶雅可比(Jacobian)椭圆函数,N是阶次。
对于奇数和偶数N,可得到不同的椭圆滤波器的幅度响应特性图,具体如下图2-3所示。
图2-3椭圆滤波器的幅度特性
其中椭圆滤波器阶次N的计算公式给出为
(2.6)
式中
,
和
。
当确定了阶次N,便可设计椭圆滤波器了。
MATLAB中提供了丰富的函数,我们可将MATLAB中一个afd_elip函数用于设计一个已知设计指标的模拟低通滤波器[3]。
2.3双线性变换法
双线性变换法保留的是从模拟到数字域的系统函数的表示,它是一种稳定的变换。
双线性变换的主要优点是S平面与Z平面是一种一一对应的单值关系,这种影射是最好的变换方法,具体变换公式为
(2.7)
这里T为采样间隔。
式(2.7)亦可转换为
(2.8)
从式(2.8)可以看出,只要其中一个变量固定,这在每个变量上都是线性的,也就是说在s和z上是双线性的。
在式(2.7)的制约下的复平面影射如图2-4所示。
从图中可以看出,
图2-4双线性变换法中的复平面影射
第一步先将整个S平面压缩影射到S1平面的
~
之间的一条横带里;第二步再通过标准变换关系将此横带转换到整个Z平面上去[3]。
由图2-4可以看出,S1平面的整个左半平面影射到单位圆内,S1平面的虚轴(整个jΩ)对应于Z平面单位圆的一周,S1平面的Ω=0处对应于Z平面的ω=0处,所以双线性变换不存在混迭效应。
为了将S平面的jΩ轴压缩到S1平面jΩ1轴上的
到
段上,可以通过如下正切变换实现
(2.9)
由此可求得Ω1。
当Ω1由
经过0变化到
时,Ω由-∞经过0变化到+∞,也即影射了整个jΩ轴。
将式(2.9)写成
(2.10)
令s=jΩ,s1=jΩ1,则可得到
(2.11)
最后通过标准变换关系
将S1平面映射到Z平面上去,便可得到式(2.7)所示的S平面和Z平面的单值影射关系。
这便是双线性变换[4]。
3设计步骤
3.1设计流程图
语音信号滤波去噪——使用双线性变换法设计的椭圆滤波器的设计流程如图3-1所示。
图3-1设计流程图
3.2录制语音信号及加噪
利用Windows平台下的录音机,录制一段语音信号“大家好,我是***”,时间在1-2s左右。
录制时,单击如图3-2所示的录音机界面的圆圈,便可录制语音信号。
录音完毕,点击如图3-2所示的方形,录制结束。
图3-2录音机界面
通过点击录音界界面上的文件,选择下拉菜单中的属性,在弹出的对话框中点击立即转换,设置录音的格式为PCM编码,属性8kHz,8位,单声道7kB/s,这样的设置在相同的时间下数据最少,处理时间也就最短,如图3-3所示。
图3-3格式设置对话框
格式设置好后,点击文件,在下拉菜单中选择另存为su.wav,并将其放入到MATLAB安装盘下的work文件夹中。
然后在MATLAB软件平台下,利用函数wavread将语音信号读入workspace,并对其进行频谱分析,画出时域和频域波形图。
函数wavreadd的调用格式为
[x,fs,nbits]=wavread(file)
其中x为采样值;fs为采样频率,单位为Hz;bits表示采样位数;file表示待读取的文件名和全路径。
采集完成后在信号中加入一频率为2200Hz的单频噪声,然后设计滤波器将信号中加入的单频噪声滤除,还原出原始信号。
下面是语音信号采集及加入单频噪声的部分程序。
[x,fs,bits]=wavread('e:
\su.wav');%输入参数为文件的全路径和文件名
N=length(x);%计算信号x的长度
fn=2200;%单频噪声频率
t=0:
1/fs:
(N-1)/fs;%计算时间范围,样本数除以采样频率
x=x';y=x+0.1*sin(fn*2*pi*t);
wavwrite(y,fs,bits,'su1');%将加噪后的语音信号写成wav文件,并命名为su1
在MATLAB平台上执行,可以通过sound(y,fs,bits)回放该语音信号,加入单频干扰后,在明显听出“大家好,我是***”的声音中伴随较尖锐的干扰啸叫声。
加入单频干扰信号后,画出干扰前后语音信号的时域波形,然后进行快速傅里叶变换,得到信号的频谱特性图。
快速傅里叶变换算法FFT计算DFT的函数为fft,其调用格式为
X=fft(x,n)
其中x为被变换的时域序列向量,n是DFT的变换区间长度,当n大于x的长度时,fft函数自动在x后面补零。
当n小于xn的长度时,fft函数计算x的前n个元素,忽略其后面的元素。
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;%计算频谱频率范围
在MATLAB平台中执行后,得到原始语音信号和加入干扰后的语音信号的时域和频域波形图如图3-4所示。
由图可以看出,加入单频干扰信号后,在时域波形图中,加入干扰后的语音信号的幅度明显比原始语音信号的大;在频域波形图中,在2200Hz处明显出现一尖锐单峰,与干扰频率刚好一致。
图3-4加噪前后语音信号时域和频域波形比较
3.3滤波器设计
滤波器的设计指标是通带截至频率fp=2000Hz,阻带截至频率fc=2150Hz,通带波纹为1dB,阻带波纹为60dB。
这里用双线性变换法设计一个满足上述指标的椭圆低通滤波器。
MATLAB中提供了一个称为[z,p,k]=ellipap(N,Rp,AS)的函数用于设计一个阶次为N,通带波纹为Rp和阻带波纹为As的归一化椭圆模拟原型滤波器。
在这里用一个不是MATLAB系统自带的称为afd_elip的函数用于设计一个已知设计指标时的模拟椭圆低通滤波器,然后再利用双线性进行模—数转换得到所要求的IIR滤波器。
afd_elip的调用格式为
[b,a]=afd_elip(Wp,Ws,Rp,As)
其中b为滤波器系统函数的分子,a为滤波器系统函数的分母,Wp为进行了预畸变的通带截至频率,Ws为进行了预畸变的阻带截至频率。
设计滤波器的部分主要程序如下
fn=2200;fs=8000;Rp=1;As=60;%低通滤波器设计指标
fp=fn-200;fc=fn-50;%定义通带和阻带截止频率
wp=fp/fs*2*pi;ws=fc/fs*2*pi;%计算对应的数字频率
T=1;%定义采样间隔
Omegap=2/T*tan(wp/2);Omegas=2/T*tan(ws/2);%截止频率预畸变
[c,d]=afd_elip(Omegap,Omegas,Rp,As)%计算滤波器系统函数分子分母系数
[b,a]=bilinear(c,d,1/T)%双线性变换得到数字滤波器系统函数分子分母系数
[db,mag,pha,grd,w]=freqz_m(b,a);%计算椭圆低通滤波器的幅度,幅度响应,脉冲响应和相位响应等数据图,验证滤波器是否达到指定性能
delta=[1,zeros(1,99)];ha=filter(b,a,delta);%计算脉冲响应
得到的椭圆滤波器的幅度、幅度响应、相位响应和脉冲响应图如图3-5所示。
图3-5滤波器的幅度、幅度响应、相位响应和脉冲响应图
图3-5中,第一幅图中的两根横线上面的横线是Rp,下面的一根横线表示As;两根竖线中,左边的一根是Wp,右边的一根表示Ws,把这四根线通过函数line标出,可以更直观得判断设计是否达标。
由图中可看出,设计的滤波器为理想的滤波器。
3.4信号滤波处理
用设计好了的滤波器对采样信号进行滤波处理,将加入的干扰信号滤除,还原出原始语音信号,并将还原出的语音信号用wavwrite函数写成wav文件。
在MATLAB,IIR滤波器利用函数filter对信号进行滤波,其调用格式为
y_fil=filter(b,a,y)
其中b为进行双线性变换后得到的数字滤波器系统函数的分子系数,a为进行双线性变换后得到的数字滤波器系统函数的分母系数。
滤波器对信号进行滤波的主要程序如下
y_fil=filter(b,a,y);%用设计好的IIR滤波器对信号进行滤波处理
Y_fil=fft(y_fil);Y_fil=Y_fil(1:
N/2);%计算频谱取前一半
wavwrite(y_fil,fs,bits,'su2');%将滤波后的语音信号写成wav文件
得到原始信号、加噪信号和滤波去噪信号的时域和频谱图如图3-6所示。
由图3-6可以
图3-6滤波前后信号的时域和频谱图
明显地看出,无论时域波形图还是频谱图,滤波后,加入的噪声被有效地抑制,达到了滤波去噪的效果,能够较好地还原出原始信号。
3.5结果分析
在MATLAB中,经sound函数对经过椭圆滤波器滤波之后的语音信号进行回放,可以听出滤波之后的信号比原始信号还清晰一些,因为周围环境的声音被滤除了一些,加入的干扰信号也被滤除了。
Sound函数的调用格式为
Sound(x,fs,bits)
其中x表示语音信号文件,fs表示采样频率,bits表示文件编码时用到的样本编码位数。
sound(x,fs,bits);%按指定的采样率和每样本编码位数回放
sound(y,fs,bits);%可以明显听出有尖锐的单频啸叫声
sound(y_fil,fs,bits);%按指定的采样率和每样本编码位数回放
在MATLAB平台上,利用前面的程序,分别插入上面的三个语句并运行程序,我们可以依次听到清晰的原始录音、伴有尖锐的啸叫声的语音和清晰的与原始语音差别不大的滤波后的语音,所得的结果证明了椭圆滤波器设计成功。
下面与两个同学的使用双线性变换设计的级联型的切比雪夫Ⅰ型滤波器和巴特沃兹滤波器的滤波效果进行比较。
加入干扰后的语音信号经这三位同学的滤波器滤波后,得到的原始信号、加噪信号和滤波去噪信号的时域和频谱图如图3-7所示。
(a)切比雪夫Ⅰ型滤波器滤
(b)巴特沃兹滤波器
图3-7双线性变换设计的不同滤波器滤波效果图
将图3-6和图3-7对比就可以看出,相同的设计指标fp=2000Hz,fc=2150Hz,Rp=1dB,As=60dB,巴特沃兹的滤波效果相对于其他三种滤波器来说滤波效果是最差的,而且需要的阶数也最高,为65阶。
但从图中,很难看出椭圆滤波器和切比雪夫Ⅰ型滤波器滤这两种滤波器的优劣,我们可以通过MATLAB的运行结果可以知道,契比雪夫Ⅰ、Ⅱ型滤波器所需要的阶数都是17阶,说明它们具有相同的性能;而椭圆滤波器具有最小的阶数8,说明它的性能在这个意义上是最优的。
3.6滤波器结构设计
利用前面的3.3节滤波器的设计中经过双线性变换得到的数字滤波器系统函数的分子分母,调用函数[C,B,A]=dir2cas(b,a),将IIR滤波器的直接型转换成级联型,运行的结果如下
[C,B,A]=dir2cas(b,a)
C=
0.0249
B=
1.00001.75521.0000
1.00000.82621.0000
1.00000.35941.0000
1.00000.21061.0000
A=
1.0000-0.00130.9731
1.0000-0.12670.8885
1.0000-0.47760.6957
1.0000-1.03240.3982
根据以上参数,将滤波器结构用visio绘图软件画好,并填上系数。
画滤波器结构图的具体操作步骤首先是打开visio绘图软件,点击文件—新建—新建绘图,进入绘图的界面。
点击绘图对话框
中的线条工具,可以跟据自己的喜好选择绘制结构图线条的大小和线段,在绘图界面中间白色的网格中便可以绘制滤波器的结构图了。
画好的级联型椭圆滤波器的结构图如图3-8所示。
图3-8级联型椭圆滤波器的结构图
4出现的问题及解决方法
从一开始接触MATLAB已经有一年多的时间了,而且已经运用MATLB软件做过了一次课程设计,但是到真正再次操作起来还是感觉很生疏,对MATLAB的提供的函数的调用还是不熟悉。
所以,这次又是在MATLAB平台上操作的课程设计,从一开始就遇到了许多问题,但在老师和同学的帮助下,当然也在自己的努力下,问题都一一得到了解决。
刚开始,因为没有电脑的录音机的麦克风坏掉了,一连换了好多台都录不进去,最后在一个麦克风是好的的同学那里好不容易录了一段音,结果因为背景有太多杂音,运行后得到的语音信号用自己设计的滤波器滤波效果不是很好,最后在下课后机房没有什么人的时候又重录了一段,这次的录音质量终于好了一些。
还有就是在画语音信号的时域和频域波形图的时候,没有对波形的坐标进行修正,观察比较的效果图不是很好。
在老师的指导下,用axis函数进行修正,便可以得到比较好看的信号时域和频谱波形图,下面图4-1是修正前和修正后的波形图。
图4-1用函数axis修正前后的波形图
在调用fld_elip等函数时,因为没有把老师提供的function文件夹放进到MATLAB工作下的work文件夹中,在运行程序的时候老是出错,提示调用的函数不存在。
最后问同学,才知道系统中没有这些函数,调用时要把这些函数的脚本放到MATLAB工作的路径下。
设计滤波器时,因为参数没有设置好,得到的滤波器滤波效果很不好,最后通过修改阻带波纹,改善了滤波器的滤波效果。
在设计滤波器的时候还遇到了一个问题,就是在画滤波器的幅度谱图的时候,需要在Wp,Ws,Rp和As处画线以便更直观判断设计得滤波器是否达标。
因为当时还不太理解标注的原理,结果因为第一次用的是下面的程序
xx=[01];yy=[00];line(xx,yy,'Color','r','LineWidth',2,'LineStyle','--');
xx=[01];yy=[-4-4];line(xx,yy,'Color','r','LineWidth',2,'LineStyle','--');
xx=[wp/piwp/pi];yy=[-9090];line(xx,yy,'Color','r','LineWidth',2,'LineStyle','--');
xx=[ws/piws/pi];yy=[-9090];line(xx,yy,'Color','r','LineWidth',2,'LineStyle','--');
运行得到的结果得到如下图4-2所示的滤波器波形图,这个图是不正确的。
图4-2-As标错的滤波器波形图
后来在老师批改了作业,并指出了自己的问题,把上面的程序改为
X_l=[0,0,wp/pi,ws/pi;1,1,wp/pi,ws/pi];Y_l=[-As,-Rp,-90,-90;-As,-Rp,3,3];
%在wp,ws,Rp,As处画线以更直观判断设计是否达标,每列参数是每个线条的端点坐标
line(X_l,Y_l,'Color','r','LineWidth',2,'LineStyle','--')%添加线条,红色,线宽为2
运行得到正确的滤波器波形图如下图4-3所示。
图4-3–As标错的滤波器波形图
5结束语
在这短短的两个星期的数字信号处理课程设计中,因为一直在不停的做老师布置的任务,感觉时间过得飞快而且充实。
这次课程设计,因为在第一个星期内就得把所有设计任务完成,同时还要把论文初稿完成提交给老师批改,所以其实也不能算是两个星期。
但是,时间虽然短,但感觉比自己一个学期学到的关于数字信号处理的知识还要多,而且动手能力也得到了大大的提高。
首先,因为之前上数字信号处理的课程,都是只听课,而没有动手操作,即使后来有实验课也因为某种原因没有将自己所学的东西在MATLAB平台上进行实践操作,致使到做这个课程设计之前对数字信号处理的意义都是不太明朗的。
如今,因为
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语音 信号 滤波 使用 双线 变换 设计 级联 椭圆 滤波器