基于MATLAB-的语谱图显示与分析

基于MATLAB的语谱图显示与分析李富强1万红1,2黄俊杰1（1.郑州大学；2.郑州解放军信息工程大学）摘要：文章给出了一种基于MATLAB的语谱图显示原理及实现方法。该法可对语音信号(或其它类型的似平稳信号)进行频谱图的伪彩色映射及显示。语谱图的类型有宽带和窄带两种，可按需求设置伪彩色显示的映射。文中介绍了语谱图显示的原理，MATLAB的相关功能函数，并给出了一个程序实例及其实验结果。在MATLAB中进行了仿真，结果表明宽带及窄带语谱图对时间分辨率和频率分辨率的不同反映，并且该法对设备要求低，编程简单，可以在普通通用计算机上实现。关键字：语谱图；宽带语谱图；窄带语谱图；伪彩色显示；MATLABTheDisplayandAnalysisofSonogramBasedonMATLABLiFuqiangWanHong（SchoolofElectronicEngineering,ZhengzhouUniversity,Zhengzhou450002China）Abstract:ThisarticleintroducesthedisplayprincipleandrealizationofspectrogrambasedonMATLAB.Thespectrogramofspeechsignalsandothersimilarquasi-stationaryonescanbedonebythismethod.Therearetwokindsofspectrogram:wide-bandandnarrow-bandspectrograms.Inaddition,differentpseudo-color-mappingdisplaysareoptional.Itcoversthedisplayprincipleofspectrogram,thecorrelativefunctionsofMATLABandthesourceprogram,aswellastheexperimentalresults.ThesimulationinMATLABturnsoutthatwide-bandandnarrow-bandspectrogramsrepresentdifferentdegreesoftimeandfrequencyresolution,respectively.ThismethodcanbeimplementedeasilywithoutcomplicatedprogrammingbasedonuniversalPC.Keywords:spectrogram；wide-bandspectrogram；narrow-bandspectrogram；pseudo-colordisplay；MATLAB1语谱图语音的发音过程中，声道通常都是处于运动状态的，因此它的共振峰特性也是时变的。不过这个时变过程比起振动过程来说要缓慢得多，因此一般可以假定它是短时平稳的，每一时刻我们都可以用这时刻附近的一短段（例如15ms）语音信号分析得到一种频谱。对语音信号连续地进行频谱分析就可以得到一种二维图谱，其横坐标表示时间，纵坐标表示频率，而每像素的灰度值大小反映相应时刻和相应频率的信号能量密度。这种时频图称为语谱图（Sonogram或Spectrogam）[1]，这种反映语音信号动态频谱特性的时频图在语音分析中有重要实用价值，被视为可视语言。从语谱图上不仅能看出任一时刻发音器官的共振峰特征，而且可以看出语音的基音频率，是否清音、爆破音等。语谱分析在语音识别、合成及编码中很有意义。1.1语谱图的产生基理语音信号是一种典型的非平稳信号，但是其非平稳性是由发音器官的物理运动过程而产生的，此过程与声波振动的速度相比较缓慢，可以假定在10～30ms这样的短时间内是平稳的。傅立叶分析是分析线性系统和平稳信号稳态特性的强有力的手段，而短时傅里叶分析，也叫时间依赖傅立叶变换，就是在短时平稳的假定下，用稳态分析方法处理非平稳信号的一种方法。设离散时域采样信号为()xn，0,1,,1nN=−L，其中n为时域采样点序号，N是信号长度。然后对信号进行分帧处理，则()xn表示为()nxm，0,1,,1mN=−L，其中n是帧序号，m是帧同步的时间序号，N为帧长（一帧内的采样点数）。信号{}()xn短时傅里叶变换为[2]：()()()jjmnmXexmwnmeωω∞−=−∞=−∑gg(1)其中{}()wn为窗序列，则信号()xn的离散时域傅里叶变换（DTFT）[3]为：10(,)()NjjmnmXnexmeωω−−==∑(2)采用离散傅立叶变换（DFT）得：210(,)()kmNjNnmXnkxmeπ−−==∑(3)其中01kN≤≤−，则(,)Xnk就是()xn的短时幅度谱估计，而时间m处频谱能量密度函数（或功率谱函数）(,)Pnk为：2(,)(,)((,))(((,)))PnkXnkXnkconjXnk==×(4)则(,)Pnk是二维的非负实值函数，并且不难证明它是信号()xn的短时自相关函数的傅里叶变换。用时间n作为横坐标，k作纵坐标，将(,)Pnk的值表示为灰度级所构成的二维图像就是语谱图。通过变换1010log((,))Pnk能得到语谱图的dB表示。1.2语谱图的伪彩色映射对(,)Pnk实施伪彩色映射，可以得到具有较高分辨率和较好视觉效果的伪彩色语谱图。把(,)Pnk的最大值max(,)Pnk映射为归一化1电平，而把最小值min(,)Pnk映射为归一化0电平。再将(,)Pnk线性映射为0～1间的电平iL，计算机监视器再跟据iL的值把语谱图以伪彩色模式显示出来。为了得到更好的显示效果，可以选择适当的基准值Base，把小于Base值限定在此基准电平上，把大于Base的值线性映射到0～1的归一化彩色值。则彩色值矩阵{}(,)Llnk=的数学表示如下：(,)(,)((,)(,))BnkBaselnkMaxBnknkBase−=∀−（5）且(,),(,)(,),(,)PnkPnkBaseBnkBasePnkBase⎧=⎨≤⎩（6）2相关MATLAB功能函数1967年美国MathWorks公司推出了MATLAB(Matrixlaborotary的简称)。MATLAB自推出之后不断更新和扩充，目前已经推出了7.x版本。它的交互式集成界面能够帮助用户快速地完成数值分析、矩阵运算、数字信号处理、仿真建模、系统控制和优化等功能。由于MATLAB的这些特性，它已经成为科研工作和工程仿真中的高效助手[4]。2.1用于短时频谱分析的功能函数1)分帧加窗函数[5]：下面程序中采用的是hanning窗。根据要求还可选用hamming(N)，blackman(N)及bartlett(N)，其中N为窗长（帧长）。2)快速傅里叶变换函数：fft(s)。s为加窗后的一帧信号。因为实时信号FFT的频域样值关于中点(即采样频率的1/2)对称，所以fft(s)给出的矩阵数据只有前半部分有用。3)复数取共扼函数：conj(z)。2.2伪彩色映射和语谱图显示的功能函数[6]1)伪彩色映射函数：colormap(MAP)。其中MAP是所采用的伪彩色映射矩阵，默认值为JET，可以通过MAP=colormap获得当前的伪彩色映射矩阵，它可以是一个任意行的矩阵，但其必须有且只有三列，并分别表示红色、绿色和蓝色的饱和度。2)频谱图显示函数：imagesc(t,f,L)。当数据比例映射使用整个色域时使用该函数，其它情况下与image()函数作用相同。其中t是时间坐标，f是频率坐标，L则是从功率谱值经伪彩色映射后的彩色电平值。3流程图及实现程序3.1流程示意图图1语谱图产生流程图3.2程序functionsogram(Winsiz,Shift,Base,Mode,Gray);[Sg,Fs]=wavread(Wavnam);n=floor((length(Sg)-Winsiz)/Shift)+1;A=zeros(Winsiz/2+1,n);fori=1:nn1=(i-1)*Shift+1;n2=n1+(Winsiz-1);s=Sg(n1:n2);s=s.*hanning(Winsiz);z=fft(s);z=z(1:(Winsiz/2)+1);z=z.*conj(z);z=10*log10(z);A(:,i)=z;endL0=(ABase);L1=(ABase);B=A.*L0+Base*L1;L=(B-Base)./(max(max(B))-Base);y=[0:Winsiz/2]*Fs/Winsiz;x=[0:n-1]*Shift;ifMode==1colormap(‘default’);elsemymode=gray;mymode=mymode(Gray:-1:1,:);colormap(mymode);endimagesc(x,y,L);axisxy;3.3程序中参数说明Wavnam：定义待分析信号的路径及文件名。语音信号由计算机声卡录制，以.wav格式存储。程序中wavread(Wavnam)读入的语音信号样值赋给矩阵Sg，也即待分析信号x(n)，采样频率赋给Fs。Winsiz：定义帧长，一般应取2的幂次，其目的是适合FFT的要求。通过对Winsiz的选值可实现宽带频谱或窄带频谱的选择。Shift：定义帧移值。一般此值小于或等于Winsiz。Shift值越小，时域分辨率越高。Base：基准电平值。注意，此值的设定需根据实际经验，可以通过在多次运行此程序中给出不同的Base值，观察所获得的频谱图的视觉和分辨率效果，选择一个合适的Base值，如果没有特别要求，可取默认值Base=0。Mode：定义显示模式。1伪彩色映射，0为灰度映射。上述程序中，伪彩色映射采用默认值’default’也即jet，可选的其它映射还有：bone，cool，copper，flag，hot，hsv，pink，prism。Gray：灰度显示层数，当Mode=0时有效，且Gray只能在1～64间取值，为了获得较好灰度显示效果，一般取值64。3语谱图显示与分析语谱图中的花纹有横杠、乱纹和竖直条等。横杠是与时间轴平行的几条深黑色带纹，它们是共振峰。从横杠对应的频率和宽度可以确定相应的共振峰频率和带宽。在一个语音段的语谱图中，有没有横杠出现是判断它是否是浊音的重要标志。竖直条（又叫冲直条）是语谱图中出现与时间轴垂直的一条窄黑条。每个竖直条相当于一个基音，条纹的起点相当于声门脉冲的起点，条纹之间的距离表示基音周期。条纹越密表示基音频率越高。图2是“你好”的灰度语谱图显示。Winsiz=256，Shift=32，Base=0，Mode=0，Gray=64。此图为灰度映射图，且窗长较短，可以获得较高的时间分辨率。在此宽带语谱图上，功率谱的谐波结构在频率轴方向被平滑掉了，谐波间隙被填平了，因此在频域看不出基音频率和频谱的谐波形式结构，但可以明显地看出共振峰结构和语谱包络，特别是时间分辨率高，可以清楚地看到浊音共振峰时变过程。图2“你好”宽带语谱图灰度显示图图3“你好”窄带语谱图灰度显示图图4“你好”宽带语谱图伪彩色显示图图3参数为：Winsiz=2048，shift=128，Base=0，Mode=0，Gray=64。此图仍为灰度映射图，但窗长较长，可以获得较高的频率分辨率。在此窄带语谱图上，滤波器的带宽能在浊音区分辨出每个谐波。语谱的谐波形式的结构在频率方向清晰可见，并能看到基音频率的时变过程。图4是“你好”的伪彩色宽带语谱图。参数为：Winsiz=256，shift=32，Base=0，Mode=1。Gray值任意。5结论文章阐述了语谱图的产生，以及如何在MATLAB中编程绘出它的伪彩色图及灰度图的方法。对相关的MATLAB函数作了简要介绍，通过对录制语音信号的分析，分别得到它的宽带语谱图及窄带语谱图，并分析了窗长对时间分辨率和频率分辨率的影响，并对宽带语谱图分别进行了灰度与伪彩色显示。通过计算机声卡，可以很方便的完成语音录制，又由于MATLAB具有很