1、 1 经典功率谱估计与现代功率谱估计的对比 摘要 本文主要介绍了在 MATLAB 环境下， 从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点以及实现方法。 关键词 功率谱；功率谱估计；经典功率谱估计；现代功率谱估计 ； 语谱图 ；共振峰 信号的频谱分析是研究信号特性的重要手段之一，通常是求其功率谱来进行频谱分析 。 功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。然而，实际应用中的平稳随机信号通常是有限长的 ，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计。 功率谱
2、估计 clear a=wavread(鸟语花香 .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; 2 相关法：相关法是利用维纳 -辛钦定理该方法先由序列 x(n)估计出自相关函数 R（ n），然后对R(n)进行傅立叶变换，便得到
3、x(n)的功率谱估计。 xn,Fs,bits=wavread(鸟语花香 .wav); 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(相关法 ); ylabel(power spectraldensity); title(recorrelation
4、psd estimate); 3 周期图法 周期图法是直接将信号的采样数据 x(n)进行 Fourier 变换求取功率谱密度估计的方法。假定有限长随机信号序列为 x(n)。它的 Fourier变换和功率谱密度估计 )(fSX 存在下面的关系： 2|)(|1)( fXNfS X 式中， N 为随机信号序列 x(n)的长度。在离散的频率点 fKf ,有： 1, . . . ,1,0,|)(|1|)(|1)( 22 NknXFFTNkXNkS X 其中， FFTx(n)为对序列 x(n)的 Fourier 变换，由于 FFTx(n)的周期为 N，求得的功率谱估计以 N 为周期，因此这种方法称为周期
5、图法。下面用例子说明如何采用这种方法进行功率谱估计 。 N=512;Nfft=256; %数据长度和 FFT 所用的数据长度 xn,Fs,bits=wavread(鸟语花香 .wav);%fs 是采样频率 n=0:N-1;t=n/Fs; %采用的时间序列 Pxx=10*log10(abs(fft(xn,Nfft).2)/N);%Fourier 振幅谱平方的平均值，并转换为 dB f=(0:length(Pxx)-1)*Fs/length(Pxx); %给出频率序列 plot(f,Pxx); %绘制功率谱曲线 xlabel(频率 /Hz);ylabel(功率谱 /dB); title(周期图 N
6、=256);grid on; 4 分段平均周期图法（ Bartlett 法） 将信号序列 x(n)， n=0,1, ,N-1，分成互不重叠的 P 个小段，每小段由 m 个采样值，则 P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列 x(n)的功率谱估 计。 平均周期图法还可以对信号 x(n)进行重叠分段，如按 2:1 重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。 这两种方法都称为平均周期图法，一般后者比前者好。 %Samp9_6 %运用信号不重叠分段估计功率谱 xn,Fs,bits=w
7、avread(鸟语花香 .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:768),Nsec).2)/Nsec; %第三段功率谱 pxx4=abs(fft(xn(769:1024),Nsec).2)/Nsec; %第四段功率谱 Pxx=10*log10(pxx1+pxx2+pxx3+pxx4)/4); %平均得
8、到整个序列功率谱 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 %运用信号重叠分段估计功率谱 pxx1=abs(fft(xn(1:256),Nsec).2)/Nsec; %第一段功率谱 pxx2=abs(fft(xn(129:384),Nsec).2)/Nsec; %第二段功率谱 p
9、xx3=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(fft(xn(769:1024),Nsec).2)/Nsec; %第七段功率谱 Pxx=10*log10(pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并
10、转化为 dB f=(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 5 通过实验仿真可以直观地看出以下特性： (1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。 (2)平均周期图法收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽分辨率低。这
11、是由于对随机6 序列的分段处理引起了长度有限所带来的 Gibbs现象而造成的。 在经典谱估计进行功率谱估计时，由于将所有在窗口外的数据都视为 0，这便使得谱估计的质量下降。现代谱估计是利用待研究信号的先验知识，对信号在窗口外的数据作出某种较合理的假设，以达 到提高谱估计质量的目的。现代功率谱估计主要是利用白噪声输入参数模型之后得到输出序列，当改变系统参数时得到的序列也不同，这样当改变参数使得输出与已知有限序列相同或者近似时就可以利用下面的公式求得其功率谱。 白噪声功率谱： Pe（ ） = 2 输出序列谱估计： Py（ ） =|H（ ej ） |2Pe（ ） = 2|H（ ej ） |2 可以说
12、现代谱估计实际上是对模型参数的估计。现代谱估计的参数模型有自回归滑动平均（ ARMA）模型、自回归（ AR）模型、滑动平均（ MA）模型， Wold分解定理阐明了三者之间的关系：任何有限方差 的 ARMA或 MA模型的平稳随机过程可以用无限阶的 MA模型表示。但是由于只有 AR模型参数估计是一组线性方程，而实际的物理系统往往是全极点系统，因而 AR应用最广。 最大熵法（ Maxmum entropy method， MEM法 ） 周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。最大熵功率谱估计的目的是最大限度地保留截断
13、后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列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 为全极点滤波器阶次； a 为全极点滤波器模型系数向量； xcorr是把 x 认为是相
14、关矩阵。 %最大熵法 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 方法估计功率谱 noverl
15、ap=128; %重叠数据 dflag=none; subplot(2,1,2) psd(xn,Nfft,Fs,window,noverlap,dflag); %采用 Welch 方法估计功率谱 xlabel(频率 /Hz);ylabel(功率谱 /dB) title(Welch 方法估计的功率谱 );grid on 7 经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不 受此限制。这是因为对于给定的 N点有限长序列 x(n)，虽然其估计出的自相关函数也是有限长的，但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度， 不像 经典谱估计那样受窗函数
16、的影响。因而现代谱的分别率比较高，而且现代谱线要平滑得多，从上图可以清楚看出。 小结： 1.功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。 2.现代谱估计方法曲线明显比经典谱估计方法光滑，说明其处理结果的方差比经典谱估计方法处理的结果小，证明了现代谱估计方 法利用参数模型对窗口外的数据的假设是比较合理有效的，能达到提高谱估计质量的目的，这也是现代谱估计优于经典谱估计的主要原因。 除此之外，我们还可以采用 多窗口法、多信号分类法 、 加窗平均周期图法 、 Welch法估计 等来功率谱估计。 语谱图 语谱图 反映了语音信号的动态频率特
17、性，在语音分析中具有重要的实用价值。水平方向是时间轴，垂直方向是频率轴，图上的灰度条纹代表各个时刻的语音短时谱。 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(语谱图 ); 8 a 女生“鸟语花香”的窄带语谱图 b 女生“鸟语花香”的宽带语谱图 9 图
18、 a 画出了一段女生“鸟语花香 ”的窄带语谱图， 因为 帧长的样点数 N=采样率 fs*窗时间长 图中使用长为 N/fs=N/fs=25/8000hz=3.1ms的汉明窗。从图中可以看出，此时浊音区的谐波非常明显，同时从图中也能分辨出共振峰的大致位置。 图 b 画出了一段女生“鸟语花香”的宽带语谱图，图中使用长为 300/8000hz=37ms的汉明窗。由于频谱的谐波结构在频率轴方向被平滑掉了，各个共振峰自成一片显得更加突出，从图中可以清楚观察到共振峰的位置以及共振峰随时间变化。 语谱图的时间分辨率和频率分辨率是由窗函数的特性决定的。时间分辨率高，可以看 出时间波形的每个周期及共振峰随时间的变化，但频率分辨率低，不足以分辨由于激励所形成的细微结构，称为宽带语谱图；而窄带语谱图正好与之相反。 宽带语谱图可以获得较高的时间分辨率，反映频谱的快速时变过程；窄带语谱图可以获得较高的频率分辨率，反映频谱的精细结构。两者相结合，可以提供与语音特性相关的信息。语谱图上因其不同的灰度，形成不同的纹路，称之为 “ 声纹 ” 。声纹因人而异，因此可以在司法、安全等场合得到应用。
