1、1,第六章 功率谱估计的经典方法,功率谱就是指功率谱密度（power spectral density,PSD），简称谱。,2,研究二阶平稳随机过程特征功率谱密度揭示随机过程中所隐含的周期及相邻的谱峰等有用信息，要用有限长的N个样本数据去估计该平稳随机过程的功率谱密度谱估计的方法。此种估计是建立在时间平均的方法之上，并假定具有遍历性。,3,经典谱估计线性、非参数化方法：周期图法，相关图法等。采用经典的傅里叶变换及窗口截断。对长序列有良好估计。 现代谱估计非线性、参数化方法：最大似然估计，最大熵法，AR模型法，预测滤波器法，ARMA模型等。对短序列的估计精度高，与经典法相互补充。是融合经典变换理
2、论、统计估计理论、系统辨识、信息论、时间序列分析及计算方法等理论与技术新学科。应用广泛，发展迅速。,4,经典谱估计现代谱估计完整的谱估计理论,长数据序列必需应用经典谱估计性能良好 短数据序列现代谱估计估计精度分辨度大大高于经典法。,5,6.1 概述,6,功率谱密度当x(t)为广义平稳过程时，其能量通常是无限的，则需研究其功率在频域上的分布，即功率谱密度。对于平稳随机过程，谱分析是采用自相关函数：,自相关函数是滞后积 的数学期望。Wiener-Kinchine定理将自相关函数与功率谱密度联系起来。,7,功率谱的物理意义及定义(Wiener-Kinchine维纳-钦辛定理)：,是以0对称，周期为2
3、。反变换为：,自相关函数是下列随机过程的数学期望。,8,6.2 自相关函数和谱密度的估计,已知随机序列 的一段样本数据 , 利用这段数据估计自相关函数的方法有： ： 为遍历性零均值广义平稳随机过程，无偏估计(对实数序列）：,该估计的特点：尽管是无偏估计，但不是一致估计， 方差很大，不是一种好的估计方法,9,：有偏估计：,该估计的特点：尽管求平均时，只用N去除不合理，但其服从渐进一致估计的原则，比无偏估计法误差小，因此工程上常用该方法计算,对于实数序列,10,利用自相关函数计算功率谱的实质：周期图,间接法：相关图法,11,即加窗截断为有限长序列，则有：,12,功率谱的估计可写成：,为数据序列,的
4、离散傅氏变换DFT。,直接法：周期图法,13,估计方差为1，长度N=32的白噪声序列的自相关函数和周期图。 实现程序为：edit ransig.m 然后编辑ransig.m文件为 运行结果。,14,function m,ms,sq,cc,cs=ransig() N=32; x=randn(1,N); disp(均值); m=mean(x); disp(均方值); ms=x*x/N; disp(方差); xstd=std(x,1); sq=xstd.*xstd; cc=xcorr(x,unbiased); figure(1); plot(x); xlabel(n); ylabel(x(n);,t
5、itle(origal signal); figure(2); plot(cc); xlabel(n); ylabel(cor); title(auto-correlation); cs=10*log10(abs(fft(x).2)/(N+1); figure(3); f=(0:length(cs)-1)/length(cs); plot(f,cs); xlabel(Frequency); ylabel(Power Spectrum (dB); title(Periodogram); grid;,15,均值 m =0.2326均方值 ms =0.8507方差 sq =0.7966,16,17,
6、6.3 相关图法 (Correlogram Method),根据Wiener-Kinchine定理，先估计出有限长信号x(n) 的自相关函数，即：,易见自相关函数是偶函数，其长度为2N-1.实际计算时，由于x(n)只有N个观测值，则对于每个延迟|m|，,可用的数据只有N-1-|m|个，故可将上式改写成：,18,第二步，求,的DFT,得到x(n)的功率谱估计：,由于功率谱是经自相关函数间接求出的，故称间接法。 平稳随机过程的自相关函数应为：,19,现在，来讨论相关图法的性能，即：,估计偏差：,分析相关图法谱估计的性能：,周期图是随机变量集合的函数，须从统计观点讨论周期图的估计性能，即上述收敛问题
7、。,这是由于：,20,所以估计的偏差为：,故可得出结论：m0，为有偏估计。 当N，m有限，偏差0，故为渐进无偏估计。,注意到：当|m|接近于N时，滞后积序列ym(n)非常短，用来计算平均值的样值数目非常少，得到的 将远离 。,21,周期图不是R(m)的傅立叶变换，为有偏估计,22,估计方差：推导过程略。(l=n-k),显然，当N，N|m|+|l|,注意到：当|m|接近于N时， 的方差将变得很大，式 不再成立。,23,m越大，分辨率越高，但自相关的偏差及方差也相应增大，通常取m=N/10 N/5，较好。,自相关函数的无偏一致估计为：,但是这个估计的方差永远大于前面的有偏一致估计，而且用它计算出的
8、功率谱出现负值，与功率谱非负的性质不符。所以工程上常用前者估计。,24,6.4 周期图法 (Periodogram Method),定义：长度为N的实平稳随机信号序列,的周期图为：,式中,25,分析周期图法谱估计的性能：,下面分析偏差和方差。,26,不难看出： 是个三角函数（两个矩形函数的卷积）， 被称为Bartlett窗函数，用 表示。,27,自相关函数的真值,28,Bartlett窗函数及其傅立叶变换示意图,1,29,30,当 时,为渐进无偏估计。,31,方差的讨论（假定x是方差为 的白噪声随机过程）。,由于周期图的方差与随机过程的4阶矩有关，计算一般随机过程的方差非常困难。但可以计算高斯
9、白噪声的周期图的方差，其结果对一般随机信号有重要参考价值。,32,33,34,矩分解定理,35,36,37,38,39,40,41,说明： 、 方差不等于零，周期图的方差不趋近于零，不满足一致估计的条件，周期图不是功率谱的一致估计。、方差较大， 不是功率谱的好估计。使得取一次估计得到的周期图的一致性和稳定性差。,说明此时周期图的方差为常数，与记录数据长度无关，即周期图的方差不会随记录数据的长度的增加而下降。,42,高斯白噪声序列的周期图,3个周期图的平均值都近似等于1，但是其方差并不随数据记录长度增加而下降。,43,结论：周期图法不是功率谱的一致估计，而且周期图都在真实功率谱附近随机起伏。,包
10、括主瓣平滑作用和旁瓣泄漏作用，,但是，增加信号的数据点数。数据量增大，相当于滞后窗加宽，相应的傅立叶变换的主瓣变窄，即功率扩散的频率范围变窄。同时，方差不会降低，随机起伏现象仍然存在。且数据点越多，随机起伏越密集。,44,随机信号：取N=64和256，观测周期图以及周期图的平均。 讨论：点数增大，滞后窗加宽，其傅立叶变换的主瓣变窄，功率扩散的频率范围变窄。,45,function =sinvnperiod() clf; %case 1 :N=64 x(n)=sin(2*pi*f*n+p)+v(n) N=64or256; n=0:N-1; f=0.2; wn1=randn(1,N); wn2=r
11、andn(1,N); wn3=randn(1,N); wn4=randn(1,N); wn5=randn(1,N); xn1=sin(2*pi*f*n)+wn1; xn2=sin(2*pi*f*n)+wn2; xn3=sin(2*pi*f*n)+wn3; xn4=sin(2*pi*f*n)+wn4; xn5=sin(2*pi*f*n)+wn5;,px1=10*log10(abs(fft(xn1).2)/(N+1); px2=10*log10(abs(fft(xn2).2)/(N+1); px3=10*log10(abs(fft(xn3).2)/(N+1); px4=10*log10(abs(f
12、ft(xn4).2)/(N+1); px5=10*log10(abs(fft(xn5).2)/(N+1); figure(1); x=(0:length(px1)-1)/length(px1); plot(x,px1,x,px2,x,px3,x,px4,x,px5); pxmean=10*log10(px1+px2+px3+px4+px5)/5); figure(2); plot(x,pxmean);,edit sinvnperiod.m,46,47,48,49,50,经典谱估计的缺点,加窗造成主瓣平滑作用，使得频率分辨率低； 加窗效应，相关图法主观认为未观测数据都等于0，造成频谱能量的泄漏；
13、 周期图法不是功率谱的一致估计，而且周期图都在真实功率谱附近随机起伏； 周期图法假设数据是以N为周期的周期性延拓，把不真实的信息加于随机过程之上，限制了频率分辨率和谱估计的质量指标，对短时间序列误差太大。,51,一些改进方法的思路,将长度为N的序列分k段进行谱估计，分别求出每一段的功率谱，然后加以平均，甚至可以每段数据有部分重叠，得平均周期图。这时，如果各段数据相互独立，则所得估计的方差为原来不分段时得1/k，达到一致估计。缺点是点数减少，分辨率下降，这样会造成有偏的。 不用矩形窗，采用合适的窗函数来消除矩形窗旁瓣带来的谱失真，但同时频谱宽度增大。即改进的周期图法。 以上两者结合为平均改进周期
14、图法，即韦尔奇法。,经典谱估计法,1,2,3,52,例如 为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为,53,与 的特性,54,平均后的周期图(每段取16个数据),55,经典谱估计法总结：,实质：根据测得的数据，通过计算 FFT（傅立叶变换）获得相应的谱密度信息。 特点：计算简单，效率高； 缺点：由于将测得数据以外的数据均看作零（实质相当于乘上了一个矩形窗口），因此频率分辨率低。频谱能量向旁瓣泄漏，且不是功率谱的一致估计。即使进行改进，也不能从根本上改善周期图的性能，只适用于数据较长和频率分辨率要求不高的场合。,56,6.5 改进的周期图法,
15、分析：滞后窗的加窗效应源自无限长序列取样时加的数据窗。 周期图的期望值被主瓣平滑的程度和功率向旁瓣泄露的多少，取决于加在信号x上的数据窗w(n)的类型。 当w(n)为矩形窗，主瓣窄，但是旁瓣泄漏也最严重。 所以，可以选择合适的窗，比如海明窗，以减弱旁瓣泄漏，使弱窄带信号也显现出来，从而提高频率分辨率，但是旁瓣的下降是以主瓣的变宽为代价的。,57,修正的周期图加窗周期图法：,定义：N为数据w(n)的长度，U为数据窗的平均能量：,修正周期图是功率谱的渐进无偏估计，但不是一致估计。,58,汉宁窗函数：海明窗函数：,59,Bartlett法周期图的平均：,由上式分析知：如果得到周期图的期望的一致估计，
16、即为功率谱密度的一致估计。而一个随机变量的一组互不相关的观测数据的算术平均是该随机变量的均值的一致估计。那么，如果对一个随机信号的若干互不相关的取样序列的周期图进行平均，得到的周期图的平均将是该随机信号的功率谱的一致估计。,60,定义：xi(n)是x(n)的K个互不相关的取样序列，每个序列长度为L。则xi(n)的周期图为则这些周期图的平均为：,性能分析略。（将近无偏估计，且在K和L都趋于无穷大时为一致估计）,61,性能分析略。（将近无偏估计，且在K和L都趋于无穷大时为一致估计）,Bartlett法：对于长度为N的一次观测x(n)，将它分成长度为L，互不重叠的K段序列，N=KL，每个子序列为：则
17、Bartlett周期图为：,62,例子：为高斯白噪声。,用平均周期图法和Bartlett法求信号的功率谱密度估计。,63,N=1024; Ns=256; pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期图 无重叠分段平均法 pxx1=abs(fft(xn(1:256),Ns).2)/Ns; pxx2=abs(fft(xn(257:512),Ns).2)/Ns; pxx3=abs(fft(xn(513:768),Ns).2)/Ns; pxx4=abs(fft(xn(769:1024),Ns).2)/Ns; pxxb=10*log10(pxx1+pxx2+pxx3+px
18、x4)/4);,64,65,66,分析结果： 与周期图谱估计比较，无重叠分段平均周期图法Bartlett法估计曲线均较为平坦，但是波峰加宽，频率分辨率下降。,67,Welch法：,对Bartlett法进行两点修正：. 让子序列xi(n)有部分重叠；. 对每个子序列都加数据窗w(n)。设子序列为xi(n) 长度为L,相邻两个子序列有L-D点重叠，于是第i个子序列为：K个子序列总长度为N=L+D(K-1)。若D=L，则N=KL，即为Bartlett法。,性能分析略。（将近无偏估计，估计方差比Bartlett法大。）,68,pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期
19、图 无重叠分段加窗平均周期图 w=hanning(256); pxx1=abs(fft(w.*xn(1:256),Ns).2)/norm(w)2; pxx2=abs(fft(w.*xn(257:512),Ns).2)/norm(w)2; pxx3=abs(fft(w.*xn(513:768),Ns).2)/norm(w)2; pxx4=abs(fft(w.*xn(769:1024),Ns).2)/norm(w)2; pxxp=10*log10(pxx1+pxx2+pxx3+pxx4)/4); 重叠分段加窗平均周期图法，即Welch法 pxx1=abs(fft(w.*xn(1:256),Ns).
20、2)/Ns; pxx2=abs(fft(w.*xn(129:384),Ns).2)/norm(w)2; pxx3=abs(fft(w.*xn(257:512),Ns).2)/norm(w)2; pxx4=abs(fft(w.*xn(385:640),Ns).2)/norm(w)2; pxx5=abs(fft(w.*xn(513:768),Ns).2)/norm(w)2; pxx6=abs(fft(w.*xn(641:896),Ns).2)/norm(w)2; pxx7=abs(fft(w.*xn(769:1024),Ns).2)/norm(w)2; pxxw=10*log10(pxx1+pxx
21、2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);,69,70,71,分析： 从上述例子看到：与Bartlett法相比，Welch法与它的周期图谱估计频率分辨率近似，但是方差减小了。 对于给定数据量N，增加子序列重叠点数，将增加子序列之间相关性，加大运算量。一般情况下，Welch法以重叠50%70%为宜。,72,描述谱估计技术的性能参数有两个： 1、变异性： 2、品质因数： ，它是变异性与分辨率之积。经典谱估计方法的每种技术的品质因数都是近似的，都与记录数据长度成反比，但每种方法的分辨率和方差都不相同（见下页表格数据）。然而，它们的总的性能指标受到数据记录长度的限制。所以，经典谱估
22、计以傅立叶变换为基础，具有计算效率高的优点，但是由于将未观测数据认为0和数据加窗，而具有频率分辨率低、旁瓣泄漏等严重的缺陷。,73,74,复习题： 1、功率谱如何定义？谱估计的任务是什么？经典谱估计和现代谱估计的主要区别是什么？经典谱估计方法有哪些优缺点？ 2、如何估计随机信号的自相关序列？自相关序列的无偏估计和渐近无偏估计有何区别？ 3、什么是周期图？什么是计算周期图的直接方法和间接方法？ 4、滞后窗对周期图的期望有何影响？什么是谱估计的频率分辨率和旁瓣泄漏？,75,5、周期图的偏差和方差如何计算？周期图的分辨率如何计算？ 6、改善周期图质量有哪些方法？ 7、Welch周期图法和Bartlett法周期图法的主要区别是什么？ 8、比较各种周期图的性能。 9、在使用Welch周期图法时，如何考虑序列分段长度和相邻分段间的重叠部分的长度的选取？,
