现代信号处理第四章-高阶谱估计的常规方法

研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学第四章高阶谱估计的常规方法4.1引言功率谱估计回顾谱估计的现代史是从Tuckey于1949年的突破开始的。Blackman-Tuckey给出了用Wiener相关法从采样数据序列得到功率谱估计的实现方法——BT法。1965年FFT的出现产生了周期图法；BT法与周期图法的致命弱点是频率分辨率的限制。为了克服这一缺点，1967年Burg受到他本人在地震应用研究中线性预测方法的启发，导出了最大熵谱估计法；E.Parzen于1968年正式提出了AR谱估计方法；此后十几年来发展了许多高分辨的谱估计方法，称为现代谱估计方法，而BT法与周期图法称为传统谱估计法。无论哪种方法，每一种谱估计技术都可以认为是一种模型法。具体地说，就是根据对过程的先验知识建立一个近似实际过程的模型，其次，利用观察数据或自相关函数来估计假设的模型参数，最后作谱估计。常用的模型有：周期图与BT（正弦谐波总和模型）、AR、MA、ARMA、Prony、最大似然等。各种估计之间的性能变化，就是由于假设的模型与实际过程拟合的好坏不同引起的。虽然不同模型可以产生类似的结果，但是有些模型所需要的参量可能相对少一些。因此，从表征过程的角度来看，这些方法就更为合适。无偏估计：ˆlim[[]]0NErr，即估计的数学期望等于其真值；一致估计：2ˆlim[()]0NErr，即估计的方差随数据增长而趋于零的估计；实际中，遇到的情况是利用有限长度数据估计一过程的高阶累积量谱。一般有两种主要方法：(1)常规（Fourier型）方法(2)参数方法：AR、MA、ARMA或Volterra模型本章讨论常规方法以及他们的统计特性和计算复杂性：常规方法可分成下列三类：(1)间接法（Indirect）：是定义（2.24）式或(3.45)式的近似；(2)直接法(Direct)：是定义（2.27）式或(3.47)式的近似(3)复解调制法。（Complexdemodulation）虽然常规法直接且其实现可利用FFT，但估计的统计方差和频率分辨率的限制对其应用能力有严格的限制。授课教师：姬红兵教授[e-mail]更新日期2010年4月1日59研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学授课教师：姬红兵教授[e-mail]更新日期2010年4月1日594.2间接法该方法是先用有限长数据估计高阶统计量，然后用多维窗函数产生高阶谱。因此，类似于用Blackman-Tuckey法估计功率谱。4.2.1高阶统计量估计设是给定的数据集，它可表示一个严格平稳随机过程或确定性序列的实现，则有以下估计方法：)(),2(),1(Nxxx步骤1.将数据分K段，每段M个样点，即N=K·M然而，如果数据样本对应一个确定性能量信号，则数据分段是不合适的。同样，如果过程是确定性周期的，则M应等于信号的周期或周期的整数倍；步骤2.每段数据去均值；步骤3.假设1,,1,0),()(Mkkxi是每段数据),,1(Ki，则高阶矩估计为)()()(1),,(1)(1)()(11)(21niisskininkxkxkxMm（4.1）其中，，，,3,2nKi,2,1,2,1,0k，),,,0max(111ns，)11,,11min(2,1MMnsM，nkL。注意到确定估计的n阶矩函数的支撑区。nL步骤4.平均所有段，即KininnxnmKm111)(11),,(1),,(ˆ，nkLn,,3,2（4.2）),,(ˆ11nxnm(),,111nnm认为是对的一般估计。如果信号是确定性的，则，即K=1，因此N=M),,(11nxnm),,(ˆ11nxnm步骤5.对于随机信号可用累积量与矩的关系式（(2.5)式）生成累积量。如果每段去均值，则有),,(ˆ11nxnc)(ˆ)(ˆ1212xxmc),(ˆ),(ˆ213213xxmc研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学授课教师：姬红兵教授[e-mail]更新日期2010年4月1日59)(ˆ)(ˆ)(ˆ)(ˆ)(ˆ)(ˆ),,(ˆ),,(ˆ12332132222321232143214xxxxxxxxmmmmmmmc（4.3）其中，,1,2,3knLkRosenblatt和VanNess[1965年]指出这种估计有两个基本要求：(1)估计应是渐近无偏的；(2)当K或1K，N时，估计的方差应趋于零，即一致估计。4.2.2高阶谱估计n阶矩谱估计：111111111ˆˆ(,,)(,,)(,,)exp{(nnnnnLLxxnnnnLLMmwj1111111111ˆˆ(,,)(,,)(,,)exp{(nnnnnLLxxnnnnLLCcwj),1nun11)}nn1n1n（4.4）而累积量谱估计为1)}nn（4.5）其中，是有界支撑域的连续窗函数，且(1uw是带宽，通常取。注意到，对于nnL/11K，可选择和nLM，以使当MLn,11(,,)nw时，有。0)/M(2Ln此外，利用高阶谱的对称性可减少（4.4）和（4.5）式的计算复杂度。4.3.3窗函数与常规功率谱估计类似，为了获得平滑的估计，必须选择合适的窗。用于高阶谱估计的窗函数12211211(,)(,)(,)(应当满足下列性质：1．具有高阶矩或累积量的对称性，例如，在实信号双谱估计中，窗函数应满足：22,)2．在估计的高阶统计量的支撑区外为零，即11(,,)0nw，,1,2,,in1Lin研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学授课教师：姬红兵教授[e-mail]更新日期2010年4月1日593．在原点等于1，即1)0,,0(w(归一化条件)4．具有实非负傅立叶变换，即0),,(11nW，i，1,,2,1ni，窗函数也应具有有限能量。利用标准一维滞后窗函数可以容易地生成一类满足这些性质的窗函数，即有1112112(,,)()()()()nnwdddd1n（4.6）其中，()d是一维函数，具有性质：()()dd()0d，nL1)0(d()0D，（4.7）然而，并不是所有的一维窗函数都满足条件：对于任意，()0D。例如：Hanning窗在频域具有负旁瓣。图4.1示出一些常用的一维窗函数。图4.2示出（4.9）式最优MSE窗，以及用（4.6）式产生的二维窗函数。这些窗函数可以用双谱偏差上确界（J）、双谱方差（V）、以及估计双谱与真实值之间的均方误差（MSE）。Rao和Gabr[1984]证明了MSE直接正比于一个效率指数E，定义为EVB，其中12121221,2BWdd2121221,2VWddSasaki等[1975]证明了双谱近似正比于指数V，双谱偏差上确界正比于指数J，21212121(),2JW2dd表4.2列出了MSE最优窗和图4.2所示的乘积窗函数对应的J，V，B值。正如所期望的，MSE最优窗具有最小的E值。Sasaki窗具有最小的J值，这是因为他最初是由使双谱偏差上确界最小化而导出的。另一方面，Parzen窗具有最小的V值。例4.1（见教材）例4.2（见教材）研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学授课教师：姬红兵教授[e-mail]更新日期2010年4月1日594.3直接法由（2.27）或（3.47）多谱的定义提供了一种估计多谱的方法。如果有有限长观察数据，并估计它的傅立叶变换系数，作三重（四重）乘积并平均即得双谱（三谱）的估计。类似于平均周期图法或Welch法。4.3.1高阶周期图（Higher-OrderPeriodogram）设{()}xk，是零均值实平稳时间序列，其DFT为1,,1,0Nk102()()expNxkFxkjNk，0,1,,1N(4.15)如果()xk没有去均值，则可置0)0(xF。高阶周期图定义为111111(,,)()()()xnnxxnxMFFFN1n(4.16)Brilliges和Rosenblatt指出，对于大的N，有111ˆ(,,)(,,)xxnnnEMM1n(4.17)并有渐近方差11112121211ˆˆvarRe(,,)varIm(,,)()()()xxnnnnxxxnnMMNMMM(4.18)其中2()xM是()xk的真实功率谱，11(,,)xnnM是()xk的真实n阶矩谱。显然，高阶周期图估计是渐近无偏的，而其方差正比于功率谱的乘积。进而，当增大时，估计方差也增大。因此，高阶周期图估计是非一致估计（inconsistent）。N通常有两种减小估计方差的途径：(1)在邻域频率上平滑周期图；（频域平滑）(2)平均不相邻时间段上的周期图估计。（Welch法，时域平滑）4.3.2矩谱的直接估计设是平稳随机或确定序列的观测值。假设采样周期)}()1({Nxx1T，nnN1是研究生课程：现代信号处理－高阶统计量分析课程编号：0211007（博）0221023（硕）西安电子科技大学授课教师：姬红兵教授[e-mail]更新日期2010年4月1日59高阶谱域沿水平或垂直方向频率样点间所需的间隔。最后，假设高阶谱是在频率0~0.5（归一化频率）之间估计。直接法如下：步骤1.将数据分为K段，每段含M个样点，即N=K*M，且每段数据去均值。如果为确定性信号，则仅用一个记录(N=M)。如果为了获得FFT算法适合的长度，可对每段数据补零。步骤2.假设(),0,1,...,1ixkkM是第i段数据，计算DFT系数1()()02()()exp{)MiixkFxkjMk，0,1,...,1,1,...,MiK(4.19)步骤3.通常和nNM之间的关系为nnNMM，其中是正奇整数，如。换言之，确定了高阶谱平滑过相邻频率的长度。由于nM21nnMJnMM为偶，而为奇，我们可以折衷的值（最近的整数）。在大小为nMnnN1nnnMMM个1n的维超矩形窗上进行频域平均估计阶矩谱：n)11(nixF((1ˆ1)(1)(1(*ixkinnnkFM)1nk()n)1kk)1n,,(1)i1111nJJJJkxFnnnnnKi,,2,1(4.20)当时，在频域不平均，有1(0)nnMJ*()()()()111111ˆ(,,)()()(),1,2,,iiiinnxxnxnMFFFiK(4.21)步骤4.最后，在K段上平均可得给定数据的阶矩谱，即nKininnxnMKM111)(11)