语音信号的短时时域分析

语音信号的时域分析就是分析和提取语音信号的时域参数。语音信号本身就是时域信号，因而时域分析是最早使用，也是应用最广泛的一种分析方法，这种方法直接利用语音信号的时域波形。时域分析通常用于最基本的参数分析及应用，如语音的分割、预处理、分类等。语音信号的时域参数有短时能量、短时过零率、短时自相关函数和短时平均幅度差函数等。这些最基本的短时参数在各种语音信号数字处理技术中都有重要的应用。

1、短时能量与短时平均幅度

设第n帧语音信号xn（m）的短时能量用En表示，则其计算公式如下：

En是一个度量语音信号幅度值变化的函数，但它有一个缺陷，即它对高电平非常敏感（因为它计算时用的是信号的平方）。为此，可采用另一个度量语音信号幅度值变化的函数，即短时平均幅度函数Mn，它定义为

Mn也是一帧语音信号能量大小的表征，它与En的区别在于计算时小取样值和大取样值不会因取平方而造成较大差异，在某些应用领域中会带来一些好处。

2、短时过零率

短时过零率表示一帧语音中语音信号波形穿过横轴（零电平）的次数。对于连续语音信号，过零即意味着时域波形通过时间轴；而对于离散信号，如果相邻的取样值改变符号则称为过零。因此，过零率就是样本改变符号的次数。定义语音信号xn（m）的短时过零率Zn为

式中，sgn[·]是符号函数，即

在实际中求过零率参数时，需要注意的一个问题是，如果输入信号中包含有50Hz的工频干扰或者A/D转换器的工作点有偏移（这等效于输入信号有直流偏移），往往会使计算的过零率参数很不准确。为了解决前一个问题，A/D转换器前的防混叠带通滤波器的低端截频应高于50Hz，以有效地抑制电源干扰。对于后一个问题除了可以采用低直流漂移器件外，也可以在软件上加以解决，这就是算出每一帧的直流分量并予以滤除。
这里应注意，在MATLAB编程中，实际上求短时平均过零率并不是按照上述公式计算，而是使用了另外一种方法。按上述过零的描述，即离散信号相邻的取样值改变符号，那它们的乘积一定为负数， 即

3、短时自相关

自相关函数具有一些性质（如它是偶函数；假设序列具有周期性，则其自相关函数也是同周期的周期函数等）。对于浊音语音可以用自相关函数求出语音波形序列的基音周期。此外，在进行语音信号的线性预测分析时，也要用到自相关函数。语音信号xn（m）的短时自相关函数Rn（k）的计算式如下：

这里K是最大的延迟点数。
短时自相关函数具有以下性质：
1）如果x（m）是周期的（设周期为N），则自相关函数是同周期的周期函数，即Rn(k)=Rn(k+Np)。
2)Rn（k)是偶函数，即Rn（k）=R（-k)。
3）当k=0时，自相关函数具有最大值，即Rn（0）≥|Rn（k）|，并且Rn（0）等于确定性信号序列的能量或随机性序列的平均功率。

4、短时平均幅度差

短时自相关函数是语音信号时域分析的重要参量。但是，计算自相关函数的运算量很大，其原因是乘法运算所需要的时间较长。利用快速傅里叶变换（FFT）等简化计算方法都无法避免乘法运算。为了避免乘法，一个简单的方法就是利用差值。为此常常采用另一种与自相关函数有类似作用的参量，即短时平均幅度差函数。

平均幅度差函数能够代替自相关函数进行语音分析的原因在于：如果信号是完全的周期信号（设周期为Np），则相距为周期的整数倍的样点上的幅值是相等的，差值为零。即：

对于实际的语音信号，d（n）虽不为零，但其值很小。这些极小值将出现在整数倍周期的位置上。为此，可定义短时平均幅度差函数：

分帧函数enframe代码如下：

function frameout=enframe(x,win,inc)

nx=length(x(:));            % 取数据长度
nwin=length(win);           % 取窗长
if (nwin == 1)              % 判断窗长是否为1，若为1，即表示没有设窗函数
   len = win;               % 是，帧长=win
else
   len = nwin;              % 否，帧长=窗长
end
if (nargin < 3)             % 如果只有两个参数，设帧inc=帧长
   inc = len;
end
nf = fix((nx-len+inc)/inc); % 计算帧数
frameout=zeros(nf,len);            % 初始化
indf= inc*(0:(nf-1)).';     % 设置每帧在x中的位移量位置
inds = (1:len);             % 每帧数据对应1:len
frameout(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));   % 对数据分帧
if (nwin > 1)               % 若参数中包括窗函数，把每帧乘以窗函数
    w = win(:)';            % 把win转成行数据
    frameout = frameout .* w(ones(nf,1),:);  % 乘窗函数
end

 分帧后计算每帧对应的时间函数为FrameTimeC代码如下：

function frameTime=FrameTimeC(frameNum,framelen,inc,fs)
% 分帧后计算每帧对应的时间
frameTime=(((1:frameNum)-1)*inc+framelen/2)/fs;

短时自相关函数STAc代码如下：

function para=STAc(X)
    para=zeros(size(X));
    fn=size(X,2);                            % 求出帧数
    wlen=size(X,1);                         %求帧长
    for i=1 : fn
        u=X(:,i);                               % 取出一帧
        R=xcorr(u);                         %短时自相关计算
        para(:,i)=R(wlen,end);          %只取k为正值的自相关函数
    end
 end

 短时平均幅度差函数STAmdf代码如下：

function para=STAmdf(X)
    para=zeros(size(X));
    fn=size(X,2);              % 求出帧数
    wlen=size(X,1);
    for i=1 : fn
        u=X(:,i);              % 取出一帧
        for k=1:wlen
            para(:,k)=sum(abs(u(k:end)-u(1:end-k+1)));      %求每个样本点的幅度差
        end
    end
 end

 短时能量计算函数STEn代码如下：

function para=STEn(x,win,inc)
    X=enframe(x,win,inc)';     % 分帧
    fn=size(X,2);              % 求出帧数
    for i=1 : fn
        u=X(:,i);              % 取出一帧
        u2=u.*u;               % 求出能量
        para(i)=sum(u2);         % 对一帧累加求和
    end
 end

短时平均幅度计算函数STMn代码如下：

function para=STMn(x,win,inc)
    X=enframe(x,win,inc)';     % 分帧
    fn=size(X,2);              % 求出帧数
    for i=1 : fn
        u=X(:,i);              % 取出一帧
        para(i)=sum(abs(u))/200;         % 对一帧累加求和
    end
 end

短时过零率计算函数STZcr代码如下：

function para=STZcr(x,win,inc)
    X=enframe(x,win,inc)';        % 分帧
    fn=size(X,2);                       % 求出帧数
    wlen=length(win);               % 求出帧长
    para=zeros(1,fn);                 % 初始化
    for i=1:fn
        z=X(:,i);                           % 取得一帧数据
        for j=1: (wlen- 1) ;            % 在一帧内寻找过零点
             if z(j)* z(j+1)< 0         % 判断是否为过零点
                 para(i)=para(i)+1;   % 是过零点，记录1次
             end
        end
    end
 end

 案例、 读入语音信号数据，计算其短时平均幅度、短时能量、短时过零率、短时幅度差以及短时自相关，并绘制成相应的图形。程序如下：

clear all; clc; close all;

[x,Fs]=audioread('C3_2_y.wav');       % 读入数据文件

wlen=200; inc=100;          % 给出帧长和帧移
win=hanning(wlen);          % 给出海宁窗
N=length(x);                    % 信号长度
time=(0:N-1)/Fs;                % 计算出信号的时间刻度
En=STEn(x,win,inc);             %短时能量
Mn=STMn(x,win,inc);             %短时平均幅度
Zcr=STZcr(x,win,inc);              %短时过零率
%此处和上述3个参数不同，返回的不是向量而是矩阵，因为一帧信号得到的不是一个数值
X=enframe(x,win,inc)';     % 分帧
xn=X(:);
Ac=STAc(X);                         %计算短时自相关
Ac=Ac(:);
Amdf=STAmdf(X);             %计算短时幅度差
Amdf=Amdf(:);

fn=length(En);             % 求出帧数

figure(1)
subplot 311; plot(time,x,'b'); axis tight% 画出时间波形 
title('(a)语音波形');
ylabel('幅值'); xlabel(['时间/s' 10 ]); 
frameTime=FrameTimeC(fn,wlen,inc,Fs);   % 求出每帧对应的时间
subplot 312; plot(frameTime,Mn,'b')     % 画出短时幅度图
title('(b)短时幅度');
ylabel('幅值'); xlabel(['时间/s' 10 ]);  
subplot 313; plot(frameTime,En,'b')     % 画出短时能量图
title('(c)短时能量');
 ylabel('幅值'); xlabel(['时间/s' 10 '(b)']);

 figure(2)
subplot 211; plot(time,x,'b'); axis tight% 画出时间波形 
title('(a)语音波形');
ylabel('幅值'); xlabel(['时间/s' 10 ]); 
subplot 212; plot(frameTime,Zcr,'b')     % 画出短时过零率图
title('(b)短时过零率');
ylabel('幅值'); xlabel(['时间/s' 10 ]);  

figure(3)
subplot 211; plot(xn,'b'); % 画出时间波形 
title('(a)语音波形');
ylabel('幅值'); xlabel(['点数' 10 ]); 
subplot 212; plot(Ac,'b')                    % 画出短时自相关图
title('(b)短时自相关');
ylabel('幅值'); xlabel(['点数' 10 ]);  

figure(4)
subplot 211; plot(xn,'b'); % 画出时间波形 
title('(a)语音波形');
ylabel('幅值'); xlabel(['点数' 10 ]); 
subplot 212; plot(Amdf,'b')                     % 画出短时幅度差
title('(b)短时幅度差');
ylabel('幅值'); xlabel(['点数' 10 ]);  

 运行结果如下：

所使用到的语音数据下载链接如下：

https://mp.csdn.net/mp_download/manage/download/UpDetailed

参考文献：语音信号处理实验教程；梁瑞宇、赵力、魏昕（编著）