2??sin12?(M)?? ?Pxx(?)?1???? (2-11)
L?????Msin???(?)?0。因?(?)是下降的，当L??时，VarP由此可见，随着L的增加VarPxxxx????此Bartlett估计是一致估计。
?(?)与式(2-1)的E?I(?)?可见在二种情况的估计量的期望值比较式(2-10)的EPNxx??都是真值Pxx(?)与窗口函数WB(ej?)的卷积形式，但后者将前者WB中N改为M，
M?N/L?N。因而使WB(ej?)主瓣的宽度增大。由于主瓣的宽度愈窄愈接近?函数，
偏倚愈小。今式(3-10)中WB(ej?)的主瓣宽度大于后式中的主瓣宽度，因而
?(?)?Bias?I(?)?，而主瓣愈宽分辨率就愈差。因此Bias可用来说明谱的分辨BiasPxxN??率，Bias愈大说明谱分辨率愈差。一个固定的记录长度N，周期图分段的数目L愈大将使方差愈小，但M也愈小，因而使Bias愈大，谱分辨率变得愈差。因此Bartlett方法中Bias或谱分辨率和估计量的方差间是有互换关系的。M和N的选择一般是由对所研究的信号的预先了解来指导的。例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M足够大。又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度N=(LM)。
由此可见Bartlett法使谱估计的方差减小是用增加Bias以及降低谱分辨率的代价换来的。实际上，当N一定时，Bias与Var的互换性是谱估计的一个固有特性。
[例如] 为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为
Pxx(z)?1 a?1?0.1
(1?az?1)(1?az)这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其
H(z)?1/(1?az?1)的滤波器形成的[6]。为了简单设选用N?8的矩形窗函数。说明N?8的
周期图可以得到的Bias的情况。图2表示M=8分4段与16段二种情况平均后的周期图。显然L=4的方差比L=16的大。M=16，L=2及8的周期图表示在图3。图3中L不同造成的影响也是明显的。但是这二个曲线的起伏都很大，因此有理由认为
6
误差主要起因于方差。比较M?8与M=16的周期图可见，M愈大，将使周期图的起
伏愈增快的结论，在这里也同样成立。
比较图2与图3发现在这个例子中最好的选择是应用L=16，M=8的估计而不是L=8、M=16的估计，即宁可减小方差，牺牲Bias。在实际中，当然功率谱密度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预先知道的。通过改变M和L以及利用预先知道的情况，通常可以很好地进行选择。
平均周期图的方法特别适合于应用FFT算法。因此在FFT出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT出现以前主要是用窗口函数处理法平滑周期图。
图1 Px(?)与E[I8(?)]的特性[6] 图2 平滑后的周期图[6]（每段取8个数据）
图3 平均后的周期图[6] (每段取16个数据)
2.4 Welch法
Welch提出了对Bartlett法的修正使更适合于用FFT进行计算。他主要提出二方面的修正，其一是选择适当的窗函数?(n)，并在周期图计算前直接加进法，这样得到的每一段的周期图为
1(i)IM(?)?MU?x(n)?(n)ein?0N?12?j?n (2-12)
7
这里U?1M?Wn?0N?12(n)为归一化因子，而Bartlett法每段的周期图为
2i?j?n1(i)IM(?)?M?x(n)en?0N?1 (2-13)
这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段时，可使各段之间有重迭，这样将会使方差减小(当N与M一定时)，重迭可以达到50％。
3 经典谱估计方法的仿真和分析
要对谱估计方法进行分析和比较，需要用到功能强大的matlab编程软件[7]对这些方法编程仿真。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连 接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 3.1周期图法仿真
周期图法又称直接法，它是把随机序列x(n)的N个观测数据视为一个能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： clear;
Fs=1000; %采样频率 n=0:1/Fs:1;
%产生含有噪声的序列改变n的取值范围观察图形的变换 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); window=boxcar(length(xn)); %矩形窗 nfft=1024;
8
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx)); xlabel('频率') ylabel('功率/DB')
20100功率/DB-10-20-30-40-500510152025频率3035404550
图4 周期图当N=100功率谱
100-10-20功率/DB-30-40-50-60-70050100150200250频率300350400450500
图5 周期图法中N=1000功率谱
9
-10-20-30功率/DB-40-50-60-70-8000.511.522.5频率33.544.5x 1054
图6 周期图中N=100000功率谱
通过matlab程序仿真从图4至图6可以看出随着采样点数的增加，该估计是渐进无偏的。从图5和图6可以看出，采用周期突发估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办法也不能吃周期图变得更加平滑，这是周期图法的缺点。周期图法得出的估计谱方差特性不好：当数据长度N太大时，扑线的起伏加剧；N太小时谱的分辨率又不好。对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段xN(n)再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但偏差加大，分辨率下降。平滑
jw是用一个适当的窗函数W(e)与算出的功率谱SX(e)进行卷积，使谱线平滑。这
jw?种方法得出的谱估计是无偏的，方差也小，但分辨率下降。 3.2 相关法谱估计法仿真
间接法是先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计的方法。功率谱图如图7所示。 Matlab代码示例： clear;
Fs=1000; %采样频率 n=0:1/Fs:1;
10