语音信号处理作业概要

华南理工大学

姓名： 学号：

班级：10级电信5班 日期：2013年5 月24日

《语音信号处理》作业报告1.实验要求

编程实现：

作业1、提取一段语音信号的短时能量、过零率、短时平均幅度差。

作业2、提取一段语音的傅里叶变换幅度谱、线性倒谱、梅尔频率倒谱（MFCC）。 作业3、 提取一段语音的LPC参数。 作业4、 估计一段语音的基音频率。

作业5、 估计一段语音的前3个共振峰频率。

作业1：

1、实验原理

（1）、短时能量

语音和噪声的区别可以体现在它们的能量上，语音段的能量比噪声段能量大，语音段的能量是噪声段能量叠加语音声波能量的和。在信噪比很高时，那么只要计算输入信号的短时能量或短时平均幅度就能够把语音段和噪声背景区分开。这是仅基于短时能量的端点检测方法。

信号{x(n)}的短时能量定义为:

语音信号的短时平均幅度定义为:

其中w(n)为窗函数。 （2）、短时平均过零率

短时过零表示一帧语音信号波形穿过横轴(零电平)的次数。过零分析是语音时域分析中最简单的一种。对于连续语音信号，过零意味着时域波形通过时间轴；而对于离散信号，如果相邻的取样值的改变符号称为过零。过零率就是样本

改变符号次数。

信号{x(n)}的短时平均过零率定义为:

式中，sgn为符号函数，即:

过零率有两类重要的应用:第一，用于粗略地描述信号的频谱特性;第二，用于判别清音和浊音、有话和无话。从上面提到的定义出发计算过零率容易受低频干扰，特别是50Hz交流干扰的影响。解决这个问题的办法，一个是做高通滤波器或带通滤波，减小随机噪声的影响；另一个有效方法是对上述定义做一点修改，设一个门限T，将过零率的含义修改为跨过正负门限。

于是，有定义:

2、实验结果及讨论

本次实验选取语音文件phrase.WAV，运行程序，结果如下图：

3、实验代码

[x,fs,nbits]=wavread('E:\\yuuyin\\phrase.WAV'); x = x / max(abs(x));%幅度归一化到[-1,1] %参数设置

FrameLen = 256; %帧长

inc = 90; %未重叠部分 amp1 = 10; %短时能量阈值 amp2 = 2;

zcr1 = 10; %过零率阈值 zcr2 = 5;

%计算过零率

tmp1 = enframe(x(1:end-1), FrameLen,inc); tmp2 = enframe(x(2:end) , FrameLen,inc); signs = (tmp1.*tmp2)<0; diffs = (tmp1 -tmp2)>0.02; zcr = sum(signs.*diffs,2);

%计算短时能量

amp = sum((abs(enframe(filter([1 -0.9375], 1, x), FrameLen, inc))).^2, 2);

subplot(3,1,1) plot(x)

axis([1 length(x) -1 1])

xlabel('帧数');ylabel('Speech');

subplot(3,1,2) plot(amp);

axis([1 length(amp) 0 max(amp)]) xlabel('帧数');ylabel('Energy'); subplot(3,1,3) plot(zcr);

axis([1 length(zcr) 0 max(zcr)]) xlabel('帧数');ylabel('ZCR');

作业2、3：

1、 提取一段语音的傅里叶变换幅度谱

[x]=wavread('E:\\yuuyin\\monologue speech_male.wav'); y=fft(x); %傅里叶变换函数 plot(abs(y)); %振幅频率 title('傅里叶变换幅度谱');

2、 提取一段语音的线性倒谱和LPC参数

基本原理：

由于频率响应H(ejw)反映声道的频率响应和被分析信号的谱包络，因此用

log|H(ejw)|做反傅里叶变换求出的LPC倒谱系数。

通过线性预测分析得到的合成滤波器的系统函数为H(z)1/(1aizi)，其

i1p冲激响应为h(n)。h(n)的倒谱为h(n)，H(z)h(n)zn1^^^^n就是说H(z)的逆变换

^h(n)是存在的。设h(0)0，将式H(z)h(n)zn两边同时对z1求导，得

n1^^^^1^nn1log[]h(n)z得到，于是有 nh(n)zi1pp11zzn1n11aiz11aiziiazii1pi1i1(1aiz)nh(n)z1i1n1^p^n1iaizi1令其左右两边z的各次幂前系数分别相

n1等，得到h(n)和ai间的递推关系

^h(1)a1^n1^ih(n)an(1)aih(nk),1npn ，按其可直接从预测系数｛ai｝求得i1p^^ih(n)(1)aih(n1),npni1倒谱h(n)。这个倒谱是根据线性预测模型得到的，又称为LPC倒谱。LPC倒谱由于利用线性预测中声道系统函数H（z）的最小相位特性，因此避免了一般同态处理中求复对数的麻烦。

^实验结果分析：

选取语音文件monologue speech_female，先读出原始语音文件的波形：

原始语音帧倒谱和预测语音帧倒谱如下：

LPC系数： ai =

Columns 1 through 8

1.0000 -0.5313 -0.1624 -0.2005 -0.0726 0.1481 0.0192 -0.1592

Columns 9 through 16

-0.1248 0.1857 0.0636 0.0848 -0.1475 -0.0562 0.0129 -0.0196

通过计算LPC系数可以很好的利用先行预测中声道系统函数的最小相位特性来提取特征函数

实验代码：

I = wavread('E:\\yuuyin\\monologue speech_female.wav');%读入原始语音 %subplot(3,1,1), plot(I);

title('原始语音波形') %对指定帧位置进行加窗处理 Q = I';

N = 256; % 窗长

Hamm = hamming(N); % 加窗 frame = 60;%需要处理的帧位置

M = Q(((frame - 1) * (N / 2) + 1):((frame - 1) * (N / 2) + N)); Frame = M .* Hamm';%加窗后的语音帧

[B,F,T] = specgram(I,N,N/2,N); [m,n] = size(B); for i = 1:m

FTframe1(i) = B(i,frame); end

P =input('请输入预测器阶数 = '); ai = lpc(Frame,P); % 计算lpc系数

LP = filter([0 -ai(2:end)],1,Frame); % 建立语音帧的正则方程 FFTlp = fft(LP);

E = Frame - LP; % 预测误差 pause

fLength(1 : 2 * N) = [M,zeros(1,N)]; Xm = fft(fLength,2 * N); X = Xm .* conj(Xm); Y = fft(X , 2 * N);

Rk = Y(1 : N);

PART = sum(ai(2 : P + 1) .* Rk(1 : P)); G = sqrt(sum(Frame.^2) - PART);

A = (FTframe1 - FFTlp(1 : length(F'))) ./ FTframe1 ; pause

%求出预测误差的倒谱

pitch = fftshift(rceps(E));

M_pitch = fftshift(rceps(Frame)); subplot(2,1,1),plot(M_pitch);grid; xlabel('语音帧');ylabel('/dB'); title('原始语音帧倒谱');

subplot(2,1,2),plot(pitch);grid; xlabel('语音帧');ylabel('/dB'); title('预测误差倒谱'); pause

3、 提取一段语音信号的梅尔频率倒谱（MFCC）

基本原理： MFCC：

语音识别和说话人识别中，常用的语音特征是基于Mel频率的倒谱系数（即MFCC）。MFCC参数是将人耳的听觉感知特性和语音的产生机制相结合。

Mel频率可以用如下公式表示：

fMel2595log(1f/700)

在实际应用中，MFCC倒谱系数计算过程如下;

① 将信号进行分帧，预加重和加汉明窗处理，然后进行短时傅里叶变换并得到其频谱。

② 求出频谱平方，即能量谱，并用M个Mel带通滤波器进行滤波；由于每一个频带中分量的作用在人耳中是叠加的。因此将每个滤波器频带内的

'x能量进行叠加，这时第k个滤波器输出功率谱(k)。

③ 将每个滤波器的输出取对数，得到相应频带的对数功率谱；并进行反

离散余弦变换，得到L个MFCC系数，一般L取12~16个左右。MFCC系数为

Cnlogx'(k)cos[(k0.5)n/M]k1M， n=1,2,...,L

④ 将这种直接得到的MFCC特征作为静态特征，再将这种静态特征做一阶和二阶差分，得到相应的动态特征。 实验结果分析：

采用语音文件monologue speech_female.wav，运行程序，得到的MFCC特征提取图像为：

通过计算MFCC参数，获得了声纹识别的特征参数。由于MFCC参数是对人耳听觉特征的描述，因此，可以认为，不同声纹的MFCC参数距离，能够代表人耳对两个语音听觉上的差异，可以为声纹的识别提供可靠的依据。

作业4：估计一段语音的基音周期 1、 实 （1）、基音周期

基音是发浊音时声带震动所引起的周期性，而基音周期是指声带震动频率的倒数。基音周期是语音信号的重要的参数之一，它描述语音激励源的一个重要特征，基音周期信息在多个领域有着广泛的应用，如语音识别、说话人识别、语音分析与综合以及低码率语音编码，发音系统疾病诊断、听觉残障者的语音指导等。因为汉语是一种有调语言，基音的变化模式称为声调，它携带着非常重要的具有辨意作用的信息，有区别意义的功能，所以，基音的提取和估计对汉语更是一个十分重要的问题。

由于人的声道的易变性及其声道持征的因人而异，而基音周期的范围又很宽，而同—个人在不同情态下发音的基音周期也不同，加之基音周期还受到单词发音音调的影响，因而基音周期的精确检测实际上是一件比较困难的事情。基音提取的主要困难反映在：①声门激励信号并不是一个完全周期的序列，在语音的头、尾部并不具有声带振动那样的周期性，有些清音和浊音的过渡帧是很难准确地判断是周期性还是非周期性的。②声道共振峰有时会严重影响激励信号的谐波结构，所以，从语音信号中直接取出仅和声带振动有关的激励信号的信息并不容 易。③语音信号本身是准周期性的(即音调是有变化的)，而且其波形的峰值点或过零点受共振峰的结构、噪声等的影响。④基音周期变化范围大，从老年男性的50Hz到儿童和女性的450Hz，接近三个倍频程，给基音检测带来了一定的困难。由于这些困难，所以迄今为止尚未找到一个完善的方法可以对于各类人群(包括男、女、儿童及不向语种)、各类应用领域和各种环境条件情况下都能获得满意的检测结果。

尽管基音检测有许多困难，但因为它的重要性，基音的检测提取一直是一个研究的课题，为此提出了各种各样的基音检测算法，如自相关函数(ACF)法、峰值提取算法(PPA)、平均幅度差函数(AMDF)法、并行处理技术、倒谱法、SIFT、谱图法、小波法等等。 （2）、自相关函数

对于离散的语音信号x(n)，它的自相关函数定义为：

R(k)=Σx(n)x(n-k)，

如果信号x(n))具有周期性，那么它的自相关函数也具有周期性，而且周期与信号x(n)的周期性相同。自相关函数提供了一种获取周期信号周期的方法。在周期信号周期的整数倍上，它的自相关函数可以达到最大值，因此可以不考虑起始时间，而从自相关函数的第一个最大值的位置估计出信号的基音周期，这使自

相关函数成为信号基音周期估计的一种工具。 （3）、短时自相关函数

语音信号是非平稳的信号，所以对信号的处理都使用短时自相关函数。短时自相关函数是在信号的第N个样本点附近用短时窗截取一段信号，做自相关计算所得的结果

Rm(k)=Σx(n)x(n-k)

式中，n表示窗函数是从第n点开始加入。

2、 实验结果分析 读取wav文件

function pitch

x=wavread('E:\\yuuyin\\isolated word.WAV');%读取声音文件 figure(1);

stem(x,'.'); %显示声音信号的波形 得到的波形如下：

利用自相关法进行基音周期估计

n=160; %取20ms的声音片段，即160个样点 for m=1:length(x)/n; %对每一帧求短时自相关函数 for k=1:n; Rm(k)=0;

for i=(k+1):n;

Rm(k)=Rm(k)+x(i+(m-1)*n)*x(i-k+(m-1)*n); end end

p=Rm(10:n); %防止误判，去掉前边10个数值较大的点 [Rmax,N(m)]=max(p); %读取第一个自相关函数的最大点 end %补回前边去掉的10个点 N=N+10;

T=N/8; %算出对应的周期

figure(2);stem(T,'.');axis([0 length(T) 0 11]);

xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');

结果如下图：

由图中可以看出基音周期大约为10ms，但是图中存在太多的野点，为此，需要对此进行进一步的处理，即去除野点。

去除野点

T1= medfilt1(T,5); %去除野点

figure(3);stem(T1,'.');axis([0 length(T1) 0 11]);

xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');

得到的结果如下：

3、 实验代码 function pitch

x=wavread('E:\\yuuyin\\isolated word.WAV');%读取声音文件 figure(1);

stem(x,'.'); %显示声音信号的波形

n=160; %取20ms的声音片段，即160个样点 for m=1:length(x)/n; %对每一帧求短时自相关函数 for k=1:n; Rm(k)=0;

for i=(k+1):n;

Rm(k)=Rm(k)+x(i+(m-1)*n)*x(i-k+(m-1)*n); end end

p=Rm(10:n); %防止误判，去掉前边10个数值较大的点 [Rmax,N(m)]=max(p); %读取第一个自相关函数的最大点 end %补回前边去掉的10个点 N=N+10;

T=N/8; %算出对应的周期

figure(2);stem(T,'.');axis([0 length(T) 0 11]);

xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');

T1= medfilt1(T,5); %去除野点

figure(3);stem(T1,'.');axis([0 length(T1) 0 11]);

xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');

作业5：估计一段语音的前3个共振峰频率

1、 实验原理 （1）共振峰的概念

共振峰是反映声道谐振特性的重要特征，它代表了发音信息的最直接来源，而且人在语音感知中利用了共振峰信息。所以共振峰是语音信号处理中非常重要的特征参数，已经广泛的应用于语音识别的主要特征和语音编码传输的基本信息。共振峰信息包含在频率包络之中，共振峰参数提取的关键是估计自然语音频谱包络，一般认为谱包络中最大值就是共振峰。 共振峰参数包括共振峰频率，频带宽度和幅值，共振峰信息包含在频率包络之中，并认为谱包络中最大值就是共振峰，利用语音频谱傅里叶变换相应的低频部分进行逆变换 ，就可以得到语音频谱的包络曲线。依据频谱包络线各峰值能量的大小确定出第一到第四共振峰。 （2）提取共振峰的方法

基于线性预测（LPC）的共振峰提取方法。一种有效的频谱包络估计方法是从线性预测分析角度推导出声道滤波器，根据这个声道滤波器找出共振峰。虽然线性预测法也有一定的缺点，如其频谱林灵敏度于人耳不想匹配。但对于许多应用来说，它仍然是一种行之有效的方法。线性预测共振峰通常有两种途径可供选择：一种途径是利用一种标准的寻找复根的程序计算预测误差滤波器的根，称为求根法；另一种途径是找出由预测其导出的频谱包络中的局部极大值，称为选峰法。

2、 实验结果分析

选取语音文件monologue speech_female，运行程序，得到的结果为： formants =

1.0e+003 *

0.686271282561051 1.7586031632567 2.592660479271902 4.000000000000000

所以该段语音的前三个共振峰频率大约为：686.27HZ、1758.6HZ、2592.7HZ 3、 实验代码

Fs=8000;

x=wavread('E:\\yuuyin\\monologue speech_female.wav') x1 = x.*hamming(length(x)); preemph = [1 0.63];

x1 = filter(1,preemph,x1); A = lpc(x1,8); rts = roots(A);

rts = rts(imag(rts)>=0);

angz = atan2(imag(rts),real(rts));

[frqs,indices] = sort(angz.*(Fs/(2*pi))); bw = -1/2*(Fs/(2*pi))*log(abs(rts(indices)));

nn = 1;

for kk = 1:length(frqs)

if (frqs(kk) > 90 && bw(kk) <400) formants(nn) = frqs(kk); nn = nn+1; end end

formants