基于短时傅里叶变换的语音信号时频分析

摘要：视频分析是近年来信号处理的新热点。本文首先介绍了语音信号STFT的相关知识，随后利用MATLAB将采集到的语音信号进行处理，并进行了信号时域和频域的相关分析。

关键词：语音信号 STFT 时频分析

语音信号的短时傅里叶变换

傅里叶变换是一种信号的整体变换，要么完全在时域，要么完全在频域进行分析处理，无法给出信号的频谱如何随时间变化的规律。而有些信号，例如语音信号，它具有很强的时变性，在一段时间内呈现出周期性信号的特点，而在另一段时间内呈现出随机信号的特点，或者呈现出两个混合的特性。对于频谱随时间变化的确定性信号以及非平稳随机信号，利用傅里叶变换分析方法有很大的局限性，或者说是不合适的。傅里叶变换无法针对性的分析相应时间区域内信号的频率特征。可以用一个窗函数与时间信号相乘积，当该窗函数的时宽足够窄，使取出的信号可以被看成是平稳信号时，就可以对乘积信号进行傅里叶变换，从而反映该时宽中的信号频谱变化规律。

早在1946年，Gabor就提出了短时傅立叶变换（Short Time Fourier Transform，STFT）的概念，用以测量声音信号的频率定位[]。

     给定一信号，其STFT定义为

     式中                                                 

及                 ，

并且窗函数应取对称函数。STFT的含义可解释如下：

    在时域用窗函数去截（注：将，的时间变量换成），对截下来的局部信号作傅立叶变换，即得在时刻得该段信号得傅立叶变换。不断地移动，也即不断地移动窗函数的中心位置，即可得到不同时刻的傅立叶变换。这些傅立叶变换的集合，即是，如图1所示。显然，是变量的二维函数。    

设语音波形分帧处理后得到的第n帧语音信号为 Xn(m)，则Xn(m)满足下式：

    其中，n＝0，1T，2T，…，并且N为帧长，T为帧移长度。某一帧的短时傅里叶变换的定义如下：

    式中w(n-m)是窗函数。不同的窗函数，可得到不同的傅里叶变换的结果。可以看出短时傅里叶变换有两个变量，即离散时间n及连续频率w。若令，则可得到离散的短时傅里叶变换如下：

它实际上就是  频率抽样。将上述某一帧语音信号的傅里叶变换写为 

    可以看出时变傅里叶变换是时间标号n的函数，当n变化时，窗函数w(n-m)沿着x(m)滑动。

                   图1 窗函数w(n-m)沿着x(m)滑动

可以得出结论：短时傅里叶变换实际就是窗选语音信号的标准傅里叶变换。这里，窗w(n-m)是一个“滑动的”窗口，它随n的变化而沿着序列X(n)滑动。由于窗口是有限长度的，满足绝对可和条件，所以这个变换是存在的。当然窗口函数不同，博里叶变换的结果也将不同。

对于w(n-m)窗来说，它除了具有选出x(m)序列中被分析部分作用外，其形状对时变傅里叶变换的特性也具有重要作用，从标准傅里叶变换可以方便的解释这种作用。如果被看成是w(n-m)x(m)序列的标准傅里叶变换，同时假设x(m)及w(m)的标准傅里叶变换存在，即：

    当n固定时，序列w(n-m)的傅里叶变换为 

    根据卷积定理，有：

因为上式右边两个卷积项均为关于角频率w的以2π为周期的连续函数，所以也可将其写成以下的卷积积分形式：

假设x(m)的DTFT是，且的DTFT是，那么是和的周期卷积。

    根据信号的时宽带宽的积为一常数这一基本性质，可知主瓣宽度与窗口宽度成反比，N越大，的主瓣越窄。为了使忠实再现的特性，相对于来说必须是—个冲激函数。所以为了使，需；但是N值太大时，信号的分帧又失去了意义。尤其是N大于语音的音素长度时，已不能反映该语音音素的频谱了。因此，应折衷选择窗的宽度N。另外，窗的形状也对短时博氏频谱有影响，如矩形窗，虽然频率分辨率很高(即主辩狭窄尖锐)，但由于第一旁瓣的衰减很小，有较大的上下冲，采用矩形窗时求得的与的偏差较大，这就是Gibbs效应，所以不适合用于频谱成分很宽的语音分析中。而汉明窗在频率范围中的分辨率较高，而且旁瓣的衰减大，具有频谱泄漏少的优点。所以在求短时频谱时一船采用具有较小上下冲的汉明窗。

    短时傅里叶变换有下列性质：  

    (1) 时移性

设  则

    (2) 频移性

设 ，则

短时傅里叶变换具有频移不变性，不具有时移不变性.但是，它在某一调制范围内即相差一相位因子的范围内保持时移不变性.

 满足“完全重构条件”：

则可由逆变换完全重构，即

语音信号的时域分析

语音信号的时域分析就是分析和提取语音信号的时域参数。进行语音分析时，最先接触到并且也是最直观的是它的时域波形。语音信号本身就是时域信号，因而时域分析是最早使用，也是应用最广泛的一种分析方法，这种方法直接利用语音信号的时域波形。时域分析通常用于最基本的参数分析及应用，如语音的分割、预处理、大分类等。

这种分析方法的特点是：表示语音信号比较直观、物理意义明确。实现起来比较简单、运算量少。可以得到语音的一些重要的参数。④只使用示波器等通用设备，使用较为简单等。

MATLAB数据采集箱中提供的进行语音信号分析的函数命令如下：

wavread :wavread 用于读取扩展名为“.wav”的声音文件。其调用形式为: y = wavread (‘filename’) 。其作用是读取wave 文件，将读取的采样数据送到y 中。    

sound:音频信号是以向量的形式表示声音采样的。sound 函数用于将向量转换为声音，其调用形式为:sound (y ,fs) ，作用是向扬声器送出向量y 中的音频信号(采样频率为fs) 。

通过wavread和plot()函数即可显示语音信号的时域波形。如图2所示。

图2  语音信号的时域波形  

贯穿于语音分析全过程的是“短时分析技术”。因为，语音信号从整体来看其特性及表征其本质特征的的参数均是随时间而变化的，所以它是一个非平稳过程，不能用处理平稳信号的数字信号处理技术对其进行分析处理。但是，由于不同的语音是由人的口腔肌肉运动产生的，相对于语音频率来说是非常缓慢的，所以从另一方面来看，虽然语音信号具有时变特性，但是在一个短时间范围内（一般可认为在10~30ms的短时间内），其特性基本保持不变即相对稳定，因而可以将其看作是一个准稳态过程，即语音信号具有短时平稳性。所以任何语音信号的分析和处理必须建立在“短时”的基础上，即进行“短时分析”，将语音信号分为一段一段来分析其特征参数，其中每一段称为一“帧”，帧长一般取为10~30ms。这样，对于整体的语音信号来讲，分析出的是由每一帧特征参数组成的特征参数时间序列。

一般每秒的帧数约为33~100帧，视实际情况而定。分帧虽然可以采用连续分段的方法，但一般要采用交叠分段的方法，这是为了使帧与帧之间平滑过渡，保持其连续性。前一帧和后一帧的交叠部分称为帧移。帧移与帧长的比值一般取为0~1/2。分帧是用可移动的有限长度窗口进行加权的方法来实现的，这就是用一定的窗函数w(n)来乘s(n)，从而形成加窗语音信号Sw(n)=s(n)w(n)。

在语音信号数字处理中常用的窗函数是矩形窗和汉明窗等，它们的表达式如下（其中N为帧长）：

矩形窗：

汉明窗：

这两种窗函数都有低通特性，通过分析这两种窗的频率响应幅度特性可以发现：矩形窗的主瓣宽度小，具有较高的频率分辨率，旁瓣峰值大，会导致泄漏现象；汉明窗的主瓣宽8*pi/N，旁瓣峰值低，可以有效的克服泄漏现象，具有更平滑的低通特性。因此在语音频谱分析时常使用汉明窗，在计算短时能量和平均幅度时通常用矩形窗。

语音信号的时域参数有短时能量、短时过零率、短时自相关函数和短时平均幅度差函数等，这是语音信号的一组最基本的短时参数，在各种语音信号数字处理技术中都要应用。

短时平均过零率分析是时域分析中最为简单的方法，即单位时间内时域波形穿过坐标轴的次数，用来区别浊音还是清音。语音中浊音过零率较低。清音的过零率较高。但由于这种高低只是相对而言，难以细化，因此用过零率来描述语言信号频率就比较粗略。

短时平均能量分析是另一种常用的时域分析方法。它也可以区分浊音和清音，清音的能量较小，浊音的能量较大。在高信噪比条件下可以区分有无语音。

短时自相关分析也是视频分析中比较重要的分析方法，常用语对基音周期的估计。在决定基音周期时，利用了短时自相关函数基音周期各整数倍的点上具有较高的峰起值。只要找到第一大峰值点的位置并计算它与原点的间隔，便能估计出基音周期。

本作业只对短时能量进行简单分析。由于语音信号的能量随时间变化，清音和浊音之间的能量差别相当显著。因此对语音的短时能量进行分析，可以描述语音的这种特征变化情况。

定义n时刻某语音信号的短时平均能量En为：

其中N为窗长,可见短时能量为一帧语音样值点的加权平方和。特殊地，当窗函数为矩形窗时，可简化为：

本次语音信号在矩形窗长N=400时的短时能量如下图3所示：

图3  N=400时的短时能量

短时平均能量的主要用途如下：

① 可以作为区分清音和浊音的特征参数。实验结果表明浊音的能量明显高于清音。通过设置一个能量门限值，可以大致判定浊音变为清音或者清音变为浊音的时刻，同时可以大致划分浊音区间和清音区间。

②在信噪比较高的情况下，短时能量还可以作为区分有声和无声的依据。

③可以作为辅助的特征参数用于语音识别中。

语音信号的频域分析

语音信号的频域分析就是分析语音信号的频域持征。从广义上讲，语音信号的频域分析包括语音信号的频谱、功率谱、倒频谱、频谱包络分析等，而常用的频域分析方法有带通滤波器组法、傅里叶变换法、线件预测法等几种。本文介绍的是语音信号的傅里叶分析法。因为语音波是一个非平稳过程，因此适用于周期、瞬变或平稳随机信号的标准傅里叶变换不能用来直接表示语音信号，而应该用短时傅里叶变换对语音信号的频谱进行分析，相应的频谱称为“短时谱 ”。

在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅立叶变换和逆变换。函数FFT用于序列快速傅立叶变换，其调用格式为y=fft(x)，其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT且和x相同长度；若x为一矩阵，则y是对矩阵的每一列向量进行FFT。函数FFT的另一种调用格式为y=fft(x,N)，式中，x，y意义同前，N为正整数。函数执行N点的FFT，若x为向量且长度小于N，则函数将x补零至长度N；若向量x的长度大于N，则函数截短x使之长度为N；若x 为矩阵，按相同方法对x进行处理。

利用上述函数即可画出语音信号的频谱图，如下图4所示。

图4 语音信号的频谱图

从图中可以确定语音参数，例如共振峰频率及基频。条纹的起点相当于声门脉冲的起点，条纹之间的距离表示基音周期。条纹越密表示基音频率越高，从图4中可以看出，当发浊音时，条纹相对密集些，其基音频率可以达到最大，说明浊音比清音的基音频率高。共振蜂表示基音脉冲的某些频率成分被加强，这在频谱图上呈现为条纹区更宽更黑。

参考文献

[1] 韩纪庆，张磊，郑铁然.语音信号处理[M].清华大学出版社，2004

[2] 易克初，田斌，付强.语音信号处理[M].国防工业出版社，2000

[3] 周辉，董正宏.数字信号处理基础及MATLAB实现[M].北京希望电子出版社，2006

[4] 邹理和.语音信号处理[M]. 北京：国防工业出版社，1985

[5] 丛玉良,王宏志.数字信号处理原理及其MATLAB实现[M].北京：电子工业出版社，2005

[6] 何强，何英.MATLAB扩展编程[M] .清华大学出版社，2002

[7] 王世一.数字信号处理[M].北京：北京理工大学出版社，2005

附录 程序

close all;        

clear all;

[x1]=wavread('xinhao.wav');  %读取语音信号的数据，赋给变量x1

figure(1);

plot(x1);                    %做原始语音信号的时域图形

title('语音信号时域波形');

xlabel('采样点');

ylabel('幅值');

n=length(x1);

y1=fft(x1,n);                %对信号做FFT变换

figure(2);

plot(abs(y1));               %原始语音信号的频域图形

axis([0 15000 0 250 ])

title('语音信号频谱图');

xlabel('Hz');

ylabel('幅值');

用MATLAB计算语音信号的短时平均能量

x=wavread('xinhao.wav');

s=fra(400,400,x);           %对输入信号进行分帧 

s2=s.^2;                    %一帧内各样点的能量 

enery=sum(s2,2);            %求一帧能量

plot(enery);                %画出N=50时的语音能量图

xlabel('帧数');             %横坐标

ylabel('短时能量E');        %纵坐标

legend('N=400');            %曲线标识

axis([0,190,5,15])

其中fra()为分帧函数，其MATLAB程序如下：

%fra.m

function f=fra(len,inc,x)    %分帧，len为帧长，inc为帧重叠样点

fh=fix(((size(x,1)-len)/inc)+1)%计算帧数

f=zeros(fh,len);            %设一个零矩阵

i=1;n=1;

while i<=fh %帧间循环

     j=1;

while j<=len %帧内循环

         f(i,j)=x(n);

         j=j+1;n=n+1;

     end

     n=n-len+inc;           %下一帧开始位置

     i=i+1;

end