基于MATLAB的语音频域特征分析

1、经典功率谱估计与现代功率谱估计的对比 摘要本文主要介绍了在MATLAB环境下，从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点以及实现方法。关键词功率谱；功率谱估计；经典功率谱估计；现代功率谱估计；语谱图；共振峰信号的频谱分析是研究信号特性的重要手段之一，通常是求其功率谱来进行频谱分析。功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计。 功率谱估计cleara=wavread

2、(鸟语花香.wav);subplot(2,1,1),plot(a);title(original signal); grid;N=256;h=hamming(N);for m=1:N b(m)=a(m)*h(m);end y=20*log(abs(fft(b);subplot(2,1,2);plot(y);title(短时谱);xlabel(频率(Hz);ylabel(功率谱(db);grid; %相关法：相关法是利用维纳-辛钦定理该方法先由序列x(n)估计出自相关函数R（n），然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 xn,Fs,bits=wavread(鸟语花香.wav)

3、;n=0:1/Fs:1;nfft=512; cxn=xcorr(xn,unbiased); CXk=fft(cxn,nfft);Pxx=abs(CXk); index=0:round(nfft/2-1);k=index*Fs/nfft; a=log(10); b=log(Pxx(index+1);c=b/a;plot_Pxx=10*c; plot(k,plot_Pxx); xlabel(frequeney（hz）相关法);ylabel(power spectraldensity);title(recorrelation psd estimate);周期图法周期图法是直接将信号的采样数据x(n)

4、进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系：式中，N为随机信号序列x(n)的长度。在离散的频率点,有：其中，FFTx(n)为对序列x(n)的Fourier变换，由于FFTx(n)的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱估计。N=512;Nfft=256; %数据长度和FFT所用的数据长度xn,Fs,bits=wavread(鸟语花香.wav);%fs是采样频率n=0:N-1;t=n/Fs; %采用的时间序列Pxx=10*log10(abs

5、(fft(xn,Nfft).2)/N);%Fourier振幅谱平方的平均值，并转换为dBf=(0:length(Pxx)-1)*Fs/length(Pxx); %给出频率序列plot(f,Pxx); %绘制功率谱曲线xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=256);grid on;分段平均周期图法（Bartlett法）将信号序列x(n)，n=0,1,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分

6、段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。%Samp9_6%运用信号不重叠分段估计功率谱 xn,Fs,bits=wavread(鸟语花香.wav);N=1024;Nsec=256;n=0:N-1;t=n/Fs; %数据点数，分段间隔，时间序列pxx1=abs(fft(xn(1:256),Nsec).2)/Nsec; %第一段功率谱pxx2=abs(fft(xn(257:512),Nsec).2)/Nsec; %第二段功率谱pxx3=abs(fft(xn(515:7

7、68),Nsec).2)/Nsec; %第三段功率谱pxx4=abs(fft(xn(769:1024),Nsec).2)/Nsec; %第四段功率谱Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4); %平均得到整个序列功率谱f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出功率谱对应的频率subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2); %绘制功率谱曲线xlabel(频率/Hz);ylabel(功率谱 /dB);title(平均周期图(无重叠) N=4*256 采用语音北风);grid on%运用信

8、号重叠分段估计功率谱 pxx1=abs(fft(xn(1:256),Nsec).2)/Nsec; %第一段功率谱pxx2=abs(fft(xn(129:384),Nsec).2)/Nsec; %第二段功率谱pxx3=abs(fft(xn(257:512),Nsec).2)/Nsec; %第三段功率谱pxx4=abs(fft(xn(385:640),Nsec).2)/Nsec; %第四段功率谱pxx5=abs(fft(xn(513:768),Nsec).2)/Nsec; %第五段功率谱pxx6=abs(fft(xn(641:896),Nsec).2)/Nsec; %第六段功率谱pxx7=abs(

9、fft(xn(769:1024),Nsec).2)/Nsec; %第七段功率谱Pxx=10*log10(pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2); %绘制功率谱曲线xlabel(频率/Hz); ylabel(功率谱/dB);title(平均周期图(重叠一半) N=1024);grid on通过实验仿真可以直观地看出以下特性：(1)功率谱估计中的相关函数法和周期图

10、法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。(2)平均周期图法收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs现象而造成的。在经典谱估计进行功率谱估计时，由于将所有在窗口外的数据都视为0，这便使得谱估计的质量下降。现代谱估计是利用待研究信号的先验知识，对信号在窗口外的数据作出某种较合理的假设，以达到提高谱估计质量的目的。现代功率谱估计主要是利用白噪声输入参数模型之后得到输出序列，当改变系统参数时得到的序列也不同，这样当改变参数使得输出与已知有限序列相同或者近似时就可以利用下面的公

11、式求得其功率谱。 白噪声功率谱：Pe（）=2 输出序列谱估计：Py（）=|H（ej）|2Pe（）=2|H（ej）|2 可以说现代谱估计实际上是对模型参数的估计。现代谱估计的参数模型有自回归滑动平均（ARMA）模型、自回归（AR）模型、滑动平均（MA）模型，Wold分解定理阐明了三者之间的关系：任何有限方差的ARMA或MA模型的平稳随机过程可以用无限阶的MA模型表示。但是由于只有AR模型参数估计是一组线性方程，而实际的物理系统往往是全极点系统，因而AR应用最广。 最大熵法（Maxmum entropy method， MEM法）周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的

12、功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),rxx(p)为基础，外推自相关序列rxx(p+1),rxx(p+2),，保证信息熵最大。最大熵功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与AR自回归（全极点滤波器）模型谱是等价的。MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem，其调用格式为：Pxx,f,a=pmem(x,p,Nfft,Fs,corr)式中，x为输入信号序列或输入相关矩阵；p为全极点滤波器阶

13、次；a为全极点滤波器模型系数向量；xcorr是把x认为是相关矩阵。%最大熵法xn,Fs,bits=wavread(鸟语花香.wav);N=1024;Nfft=256;n=0:N-1;t=n/Fs; %数据长度、分段长度和时间序列window=hanning(256); Pxx1,f=pmem(xn,14,Nfft,Fs); %采用最大熵法，采用滤波器阶数14，估计功率谱subplot(2,1,1),plot(f,10*log10(Pxx1); %绘制功率谱xlabel(频率/Hz);ylabel(功率谱/dB);title(最大熵法(MTM) Order=14);grid on%采用Welch

14、方法估计功率谱noverlap=128; %重叠数据dflag=none;subplot(2,1,2)psd(xn,Nfft,Fs,window,noverlap,dflag); %采用Welch方法估计功率谱xlabel(频率/Hz);ylabel(功率谱/dB)title(Welch方法估计的功率谱);grid on经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N点有限长序列x(n)，虽然其估计出的自相关函数也是有限长的，但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，不像经典谱估计那样受窗函数的影响。因而

15、现代谱的分别率比较高，而且现代谱线要平滑得多，从上图可以清楚看出。小结：1.功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。2.现代谱估计方法曲线明显比经典谱估计方法光滑，说明其处理结果的方差比经典谱估计方法处理的结果小，证明了现代谱估计方法利用参数模型对窗口外的数据的假设是比较合理有效的，能达到提高谱估计质量的目的，这也是现代谱估计优于经典谱估计的主要原因。除此之外，我们还可以采用多窗口法、多信号分类法、加窗平均周期图法、Welch法估计等来功率谱估计。语谱图语谱图反映了语音信号的动态频率特性，在语音分析中具有重要的实用价值。水平

16、方向是时间轴，垂直方向是频率轴，图上的灰度条纹代表各个时刻的语音短时谱。x,fs,nbits=wavread(鸟语花香.wav); %fs是一个标量，指定采样频率specgram(x,512,fs,300); % specgram(a,nfft,fs,window), window is a periodic Hann (Hanning) window of length nfft.xlabel(时间(s);ylabel(频率(Hz);title(语谱图);a女生“鸟语花香”的窄带语谱图b女生“鸟语花香”的宽带语谱图图a画出了一段女生“鸟语花香”的窄带语谱图，因为帧长的样点数N=采样率fs*窗时间长图中使用长为N/fs=N/fs=25/8000hz=3.1ms的汉明窗。从图中可以看出，此时浊音区的谐波非常明显，同时从图中也能分辨出共振峰的大致位置。图b画出了一段女生“鸟语花香”的宽带语谱图，图中使用长为300/8000hz=37ms的汉明窗。由于频谱的谐波结构在频率轴方向被平滑掉了，各个共振峰自成一片显得更加突出，从图中