语音信号处理系统设计

课题六语音信号处理系统设计

摘要：语音信号处理是研究用数字信号处理技术对语音信号进行处理的一门学科。语音信号处理的目的是得到某些参数以便高效传输或存储,或者是用于某种应用，如人工合成出语音、辨识出讲话者、识别出讲话内容、进行语音增强等。本文简要介绍了语音信号采集与分析以及语音信号的特征、采集与分析方法，并在采集语音信号后，在MATLAB 软件平台上进行频谱分析,并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。利用MATLAB来读入（采集）语音信号，将它赋值给某一向量，再将该向量看作一个普通的信号，对其进行FFT变换实现频谱分析，再依据实际情况对它进行滤波，然后我们还可以通过sound命令来对语音信号进行回放，以便在听觉上来感受声音的变化。

关键词：Matlab，语音信号，傅里叶变换，滤波器

1课程设计的目的和意义

本设计课题主要研究语音信号初步分析的软件实现方法、滤波器的设计及应用。通过完成本课题的设计，拟主要达到以下几个目的：

1.1．了解Matlab软件的特点和使用方法。

1.2．掌握利用Matlab分析信号和系统的时域、频域特性的方法；

1.3．掌握数字滤波器的设计方法及应用。

1.4．了解语音信号的特性及分析方法。

1.5．通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。

2 设计任务及技术指标

设计一个简单的语音信号分析系统，实现对语音信号时域波形显示、进行频谱分析，利用滤波器滤除噪声、对语音信号的参数进行提取分析等功能。采用Matlab设计语言信号分析相关程序，并且利用GUI设计图形用户界面。具体任务是：

2.1．采集语音信号。

2.2．对原始语音信号加入干扰噪声，对原始语音信号及带噪语音信号进行时频域分析。

2.3．针对语音信号频谱及噪声频率，设计合适的数字滤波器滤除噪声。

2.4．对噪声滤除前后的语音进行时频域分析。

2.5.对语音信号进行重采样，回放并与原始信号进行比较。

2.6．对语音信号部分时域参数进行提取。

2.7．设计图形用户界面（包含以上功能）。

3 设计方案论证

3.1语音信号的采集

使用电脑的声卡设备采集一段语音信号，并将其保存在电脑中。

3.2语音信号的处理

语音信号的处理主要包括信号的提取播放、信号的重采样、信号加入噪声、信号的傅里叶变换和滤波等，以及GUI图形用户界面设计。

Ⅰ.语音信号的时域分析

语音信号是一种非平稳的时变信号，它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在与方便有效的提取并表示语音信号所携带的信息。语音信号分析可以分为时域和变换域等处理方法，其中时域分析是最简单的方法。

Ⅱ.语音信号的频域分析

信号的傅立叶表示在信号的分析与处理中起着重要的作用。因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅立叶分析方法能完善地解决许多信号分析和处理问题。另外，傅立叶表示使信号的某些特性变得更明显，因此，它能更深入地说明信号的各项红物理现象。

由于语音信号是随着时间变化的，通常认为，语音是一个受准周期脉冲或随机噪声源

激励的线性系统的输出。输出频谱是声道系统频率响应与激励源频谱的乘积。声道系统的频率响应及激励源都是随时间变化的，因此一般标准的傅立叶表示虽然适用于周期及平稳随机信号的表示，但不能直接用于语音信号。由于语音信号可以认为在短时间内，近似不变，因而可以采用短时分析法。 Ⅲ.模拟滤波器原理

（1）模拟巴特沃思滤波器原理

巴特沃斯滤波器具有单调下降的幅频特性：在小于截止频率c Ω的范围内，具有最平幅度的响应，而在c Ω>Ω后，幅频响应迅速下降。 巴特沃思低通滤波器幅度平方函数为：

221()1()

a N c

H j Ω=

Ω+Ω (2-1) 式中N 为滤波器阶数，c Ω为3dB 截止角频率。将幅度平方函数写成s 的函数：

21

()()1()

a a N

c

H s H s s j -=

+Ω (2-2) 该幅度平方函数有2N 个等间隔分布在半径为c Ω的圆上的极点121()

22 k j N

k c s e

π++=Ω，

0,1,...21k N =- 为了形成稳定的滤波器，取左半平面的N 个极点构成()a H s ，即：

1

()()N N a c

k

k H s s s -==Ω

-∏ (2-3)

为使设计统一，将频率归一化，得到归一化极点121()

22k j N

k p e π++=，相应的归一化系

统函数为：

1

()1

()N a k

k H p p p -==-∏ (2-4)

多项式形式为： 01()1(...)N a H p b b p p =+++ (2-5)

（2）模拟切比雪夫滤波器原理

切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是I 型滤波器，在通带内单调、在阻带内等波纹的是II 滤波器。以I 型滤波器为例。

切比雪夫滤波器的幅度平方函数为：

2222

1

()()1(

)a

N p

A H j C εΩ=Ω=

Ω

+Ω (2-6)

ε为小于1的正数，表示通带内幅度波动的程度。Ωp 称为通带截止频率。令λ=Ω/Ωp ，称为对Ωp 的归一化频率。C N (x)为N 阶切比雪夫多项式。幅度平方函数的极点是分布在b Ωp 为长半轴，a Ωp 为短半轴的椭圆上的点。同样取s 平面左半平面的极点构成()a H s ：

1

1()2

()N

N N a p

i i H s s s ε-==Ω

??-∏ (2-7)

进行归一化，得到：1

1

()12

()N

N a i i H p p p ε-==??-∏ (2-8)

其中(21)(21)sin[

]cos[]22k k p ch jch i N N

ππ

ξξ--=-?+? ，11()Arsh N ξε=

（3）模拟滤波器数字化原理

将模拟滤波器转化为数字滤波器在工程上常用的有脉冲响应不变法和双线性变换法。 脉冲响应不变法时一种时域上的转换方法，它是数字滤波器的单位取样响应在抽样点上等于模拟滤波器的单位冲激响应，即：

()()h n h nT a = (2-9)

设模拟滤波器只有单阶极点，其系统函数为：

1

()N

i

a i i

A H s s s ==-∑

(2-10) 对()a H s 进行拉氏反变换得到()a h t ，对()a h t 进行等间隔抽样，得到()()h n h nT a =，对

()h n 进行Z 变换，得到数字滤波器系统函数：

11()1i N

i

s T i A H z e

z -==-∑

(2-11) 这种方法s 和z 的关系是：sT z e =。该方法的优点是频率坐标变换时线性的切数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域特性逼近好；缺点是会产生频谱混叠现象，适合低通、带通滤波器的设计，不适合高通、带阻滤波器的设计。 双线性变换法为了克服频谱混叠现象，采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到/T π±之间，再用sT z e =转换到Z 平面上。

这种方法s 和z 的关系是：11(2/)(1/1)s T z z --=-+。该方法克服了频谱混叠现象，但带来了频率坐标变换的非线性：(2/)tan(/2)T ωΩ=，由模拟滤波器系统函数转换为数字滤波器系统函数公式为：

1211()()|

a z z T z H z H s --=

+= (2-12)

3.3语音信号的效果显示

图形用户界面（Graphical User Interface ，GUI ）是由窗口、按键、菜单、文字说明等对象（Objects ）构成的一个用户界面。用户通过一定的方法（如鼠标、键盘）选择激活这些图形对象，实现计算、绘图等。 创建图形用户界面须具有三类基本元素：

（1）组件：图形化控件（如按钮、编辑框、列表框等）、静态元素（如文本字符串）、菜单和坐标系。

（2）图形窗口：GUI 的每一个组件都须安排在图形窗口中。

（3）回应：如用户用鼠标单击或用键盘输入信息后，程序要有相应的动作。 3.4短时能量及短时平均幅度

短时能量函数和短时平均幅度函数是基于语音信号幅度的变化。清音段幅度小，其能量集中于高频段；浊音段幅度较大，其能量集中于低频段。短时能量函数对信号电平值过于敏感。由于需要计算信号样值的平方和，在实际应用中（如定点设备）很容易溢出，所以可以用平均幅度函数来代替短时能量函数。 短时能量函数定义：

∑-==1

2

)(N m n n m x E （6-1）

短时平均幅度函数定义：

∑-==

1

|)(|N m n n m x M （6-2）

（2）作用

(a) 区分清/浊音：En 、Mn 大，对应浊音； En 、Mn 小，对应清音。

(b) 在信噪比高的情况下，能进行有声/无声判决：无声时，背景噪声的En 、Mn 小；有声时，En 、Mn 显著增大。判决时可设置一个门限。 3.5短时平均过零率

（1）过零率定义：信号跨越横轴的情况。

对于连续信号，观察语音时域波形通过横轴的情况；对于离散信号，相邻的采样值具有不同的代数符号，也就是样点改变符号的次数。 （2）短时平均过零率

对于语音信号，是宽带非平稳信号，考察其短时平均过零率。

∑-=--=1

)]1(sgn[)](sgn[21N m n n n m x m x Z （6-3）

其中sgn[.]为符号函数，

?

?

?<=≥=0)(1))(sgn(0

)(1))(sgn(n x n x n x n x （6-4） （3）作用

(a) 区分清/浊音：浊音平均过零率低，集中在低频端；清音平均过零率高，集中在高频端。

(b) 从背景噪声中找出是否有语音，以及语音的起点。

3.6短时自相关函数

相关函数用于确定两个信号在时域内的相似性。常用的物理量为自相关函数和互相关函数。当两个信号的互相关函数大时，则说明一个信号可能是另一个信号的时间滞后或提前；当互相关函数为0时，则两个信号完全不同。

自相关函数用于研究信号本身，如波形的同步性和周期性。 （1）自相关函数定义

K k k m x m x k R k

N m n n n ≤<+=

∑--=0)

()()(10

（5-5）

（2）自相关函数特点

(a) 当时域信号为周期信号时，自相关函数也是周期性函数，两者具有同样的周期。 (b) Rn(k)为偶函数， Rn(k)＝Rn(－k)

(c) Rn(0)最大， Rn(0)≥ |Rn(k)|, Rn(0)=En,对于确定信号， Rn(0)是信号能量；对于随机信号或周期信号， Rn(0)是平均功率。 （3）作用

(a) 区分清/浊音。浊音语音的自相关函数具有一定的周期性。清音语音的自相关函数不具有周期性，类似噪声。

(b) 估计浊音语音信号的周期，即估计基音周期。 3.7系统初步流程图

下图列出了整个语音信号处理系统的工作流程：

4 设计内容

4.1语音信号的打开

程序：

[x,fs,nbits]=wavread('mei.wav'); %打开语音信号

sound(x,fs,nbits); %播放语音信号

N=length(x); %长度

n=0:N-1;

w=2*n*pi/N;

y1=fft(x); %对原始信号做FFT变换subplot(2,1,1);

plot(n,x) %做原始语音信号的时域波形图title('原始语音信号时域图');

xlabel('时间t'); ylabel('幅值');

subplot(2,1,2); %做原始语音信号的频谱图 plot(w/pi,abs(y1)); title('原始语音信号频谱') xlabel('频率Hz'); ylabel('幅度'); 程序结果：

原始语音信号时域图

时间t

幅值

原始语音信号频谱

频率Hz

幅度

-----------------------------------------------------------------------------------------------------------------------------------------

4.2 信号重采样 程序：

[x,fs,nbits]=wavread('mei.wav'); %打开语音信号 x1=resample(x,1,2); %进行8KHZ 重采样 sound(x1,fs/2,nbits); %读重采样后的信号

x2=resample(x,1,4); %进行4KHZ重采样

sound(x2,fs/4,nbits); %读重采样后的信号

n1=0:N/2-1;

n2=0:N/4-1;

w1=4*n1*pi/N;

w2=8*n2*pi/N;

y1=fft(x1); %对重采样信号做FFT变换

y2=fft(x2);

figure(1)

subplot(2,1,1);

plot(n1,x1) %做8khz重采样语音信号的时域波形图title('8KHZ重采样语音信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,1,2); %做8KHZ重采样语音信号的频谱图plot(w1/pi,abs(y1));

title('8KHZ重采样语音信号频谱')

xlabel('频率Hz');

ylabel('幅度');

figure(2)

subplot(2,1,1);

plot(n2,x2) %做4KHZ 重采样语音信号的时域波形图 title('4KHZ 重采样语音信号时域图'); xlabel('时间t'); ylabel('幅值');

subplot(2,1,2); %做4KHZ 重采样语音信号的频谱图 plot(w2/pi,abs(y2));

title('4KHZ 重采样语音信号频谱') xlabel('频率Hz'); ylabel('幅度');

程序结果：

8KHZ 重采样语音信号时域图

时间t

幅值

0.2

0.4

0.6

0.8

1 1.2

1.4

1.6

1.8

2

8KHZ 重采样语音信号频谱

频率Hz

幅度

时间t

幅值

频率Hz

幅度

-----------------------------------------------------------------------------------------------------------------------------------------

4.3加噪音 程序：

x=wavread('mei.wav')'; %打开语音信号 N=length(x); %长度 n=0:N-1; w=2*n*pi/N;

z=0.02*cos(10*n); %噪音 sound(z,fs,nbits); %读噪音 y=x+z; %加噪音

sound(y,fs,nbits); %读加噪音后的信号 y1=fft(z); %对噪音做FFT 变换 y2=fft(y); %对加噪信号做FFT 变换 figure(1)

subplot(2,1,1); %做噪音信号的时域波形图

plot(n,z);

axis([0,50,-0.02,0.02]);

title('噪音信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,1,2); %做噪音信号的频谱图plot(w/pi,abs(y1));

title('噪音信号频谱')

xlabel('频率Hz');

ylabel('幅度');

figure(2)

subplot(2,1,1); %做加噪音信号的时域波形图

plot(n,y);

axis([0,50,-0.02,0.02]);

title('加噪音信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,1,2); %做加噪音信号的频谱图plot(w/pi,abs(y2));

title('加噪音信号频谱')

xlabel('频率Hz'); ylabel('幅度');

程序结果：

噪音信号时域图

时间t 幅值

噪音信号频谱

频率Hz

幅度

加噪音信号时域图

时间t 幅值

加噪音信号频谱

频率Hz

幅度

-----------------------------------------------------------------------------------------------------------------------------------------

4.4滤波器

function [db,mag,pha,w]=freqz_m(b,a);

[H,w]=freqz(b,a,1000,'whole'); %在0-2*pi 之间选取N 个点计算频

率响应

H=(H(1:501))'; %频率响应 w=(w(1:501))'; %频率 mag=abs(H); %响应幅度 db=20*log10((mag+eps)/max(mag)); %增益 pha=angle(H);

-----------------------------------------------------------------------------------------------------------------------------------------

function [b,a]=afd_butt(Wp,Ws,Rp,As)

N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws))); %滤波器阶数N 为整数

fprintf('\n Butterworth Filter Order=%2.0f\n',N)

OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N))) %求对应于N 的3db 截止频率

[b,a]=u_buttap(N,OmegaC);

-----------------------------------------------------------------------------------------------------------------------------------------

function [b,a]=imp_invr(c,d,T)

[R,p,k]=residue(c,d); %部分分式展开

p=exp(p*T); %从模拟到数字极点对应关系sT e z ,部分分式系数相同

[b,a]=residuez(R,p,k); %将部分分式的形式变换成多项式之比的形式

b=real(b'); %求出数字滤波器系数 a=real(a');

----------------------------------------------------------------------------------------------------------------------------------------- 程序：

wp=0.6*pi;ws=0.75*pi;Rp=1;As=30;T=1; %定义滤波器参数

OmegaP=wp/T;OmegaS=ws/T;

[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);

[b,a]=imp_invr(cs,ds,T);

[db,mag,pha,w]=freqz_m(b,a);

figure(1)

subplot(2,1,1);plot(w/pi,mag); %作出滤波器的图形

title('digital filter Magnitude Response')

axis([0,1,0,1.1])

subplot(2,1,2);plot(w/pi,db);

title('digital filter Magnitude in DB')

axis([0,1,-40,5]);

x=wavread('mei.wav')'; %打开语音信号

N=length(x); %长度

n=0:N-1;

w=2*n*pi/N;

z=0.02*cos(10*n); %噪音

y=x+z;

y1=filter(b,a,y); %用滤波器滤波

y2=fft(y1); %对滤波后信号做FFT变换sound(y1); %读滤波后信号

figure(2)

subplot(2,1,1);

plot(n,y1) %做滤波后信号的时域波形图title('滤波后信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,1,2); %做滤波后信号的频谱图plot(w/pi,abs(y2));

title('滤波后信号频谱')

xlabel('频率Hz');

ylabel('幅度');

digital filter Magnitude Response

digital filter Magnitude in DB

滤波后信号时域图

时间t 幅值

滤波后信号频谱

频率Hz

幅度

-----------------------------------------------------------------------------------------------------------------------------------------

4.5 GUI 的函数 1 开语音

global x;global fs; global nbits; global w;global N; H={'*.wav'};

[x,fs,nbits]=wavread('mei.wav'); y=fft(x);

N=length(x); n=0:N-1; w=2*n*pi/N; subplot(2,1,1);

plot(n,x) title('原始语音信号时域图'); xlabel('时间t');

ylabel('幅值');

subplot(2,1,2);

plot(w/pi,abs(y));

title('原始语音信号频谱')

xlabel('频率Hz');

ylabel('幅度');

点开语音键结果：

----------------------------------------------------------------------------------------------------------------------------------------- 2 读语音

global fs; global x;

sound(x,fs);

-----------------------------------------------------------------------------------------------------------------------------------------

3 重采样

global x;

global fs;

global nbits;

N=length(x);

x1=resample(x,1,2);

sound(x1,fs/2,nbits);

x2=resample(x,1,4);

sound(x2,fs/4,nbits);

n1=0:N/2-1;

n2=0:N/4-1;

w1=4*n1*pi/N;

w2=8*n2*pi/N;

y11=fft(x1);

y12=fft(x2);

subplot(2,2,1);

plot(n1,x1) %8kh重采样语音信号时域图

title('8KH重采样语音信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,2,2); %8KH重采样语音信号频域图

plot(w1/pi,abs(y11));

title('8KH重采样语音信号频域图')

xlabel('频率Hz');

ylabel('幅度');

subplot(2,2,3);

plot(n2,x2) %4KH重采样语音信号时域图

title('4KH重采样语音信号时域图');

xlabel('时间t');

ylabel('幅值');

subplot(2,2,4); %4KH重采样语音信号频域图

plot(w2/pi,abs(y12));

title('4KH重采样语音信号频域图')

xlabel('频率Hz');

ylabel('幅度');

点重采样键结果：