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

研究生课程：现代信号处理高阶统计量分析 课程编号：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 0NE rr，即估计的数学期望等于其真值；一致估计：2lim()0NE rr，即估计的方差随数据增长而趋于零的估计；实际中，遇到的情况是利用有限长度数据估计一过程的高阶累积量谱。一般有两种主要方法：(1)常规（Fourier型）方法(2)参数方法：AR、MA、ARMA 或Volterra模型 本章讨论常规方法以及他们的统计特性和计算复杂性：常规方法可分成下列三类：(1)间接法（Indirect）：是定义（2.24）式或(3.45)式的近似；(2)直接法(Direct)：是定义（2.27）式或(3.47)式的近似(3)复解调制法。（Complex demodulation）虽然常规法直接且其实现可利用FFT，但估计的统计方差和频率分辨率的限制对其应用能力有严格的限制。授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 594.2 间接法 4.2 间接法 该方法是先用有限长数据估计高阶统计量，然后用多维窗函数产生高阶谱。因此，类似于用Blackman-Tuckey法估计功率谱。4.2.1高阶统计量估计 设是给定的数据集，它可表示一个严格平稳随机过程或确定性序列的实现，则有以下估计方法：)(),2(),1(Nxxx步骤1.将数据分 K 段，每段 M 个样点，即 N=KM 然而，如果数据样本对应一个确定性能量信号，则数据分段是不合适的。同样，如果过程是确定性周期的，则 M 应等于信号的周期或周期的整数倍；步骤2.每段数据去均值；步骤3.假设1,1,0),()(Mkkxi是每段数据),1(Ki，则高阶矩估计为)()()(1),(1)(1)()(11)(21niisskininkxkxkxMm （4.1）其 中，,3,2nKi,2,1,2,1,0k，),0max(111ns，)11,11min(2,1MMnsM，nkL。注意到确定估计的n阶矩函数的支撑区。nL步骤4.平均所有段，即 KininnxnmKm111)(11),(1),(，nkLn,3,2 （4.2）),(11nxnm(),111nnm认为是对的一般估计。如果信号是确定性的，则，即K=1，因此 N=M),(11nxnm),(11nxnm步骤5.对于随机信号可用累积量与矩的关系式（(2.5)式）生成累积量。如果每段去均值，则有),(11nxnc )()(1212xxmc),(),(213213xxmc 研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59)()()()()()(),(),(12332132222321232143214xxxxxxxxmmmmmmmc （4.3）其中，,1,2,3knL k Rosenblatt和Van Ness1965年指出这种估计有两个基本要求：(1)估计应是渐近无偏的；(2)当K 或1K，N时，估计的方差应趋于零，即一致估计。4.2.2 高阶谱估计 n阶矩谱估计：111111111(,)(,)(,)exp(nnnnnLLxxnnnnLLMmwj111111111 1(,)(,)(,)exp(nnnnnLLxxnnnnLLCcwj),1nun11)nn1n1n （4.4）而累积量谱估计为 1)nn （4.5）其中，是有界支撑域的连续窗函数，且(1uw是带宽，通常取。注意到，对于nnL/11K，可选择和nLM，以使当MLn,11(,)nw时，有。0)/M(2Ln此外，利用高阶谱的对称性可减少（4.4）和（4.5）式的计算复杂度。4.3.3 窗函数 与常规功率谱估计类似，为了获得平滑的估计，必须选择合适的窗。用于高阶谱估计的窗函数12211211(,)(,)(,)(wwww应当满足下列性质：1具有高阶矩或累积量的对称性，例如，在实信号双谱估计中，窗函数应满足：22,)2在估计的高阶统计量的支撑区外为零，即 11(,)0nw，,1,2,in1L in 研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 593在原点等于1，即1)0,0(w (归一化条件)4具有实非负傅立叶变换，即0),(11nW，i，1,2,1ni，窗函数也应具有有限能量。利用标准一维滞后窗函数可以容易地生成一类满足这些性质的窗函数，即有 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，定义为EV B，其中 12121221,2BWdd 2121221,2VWdd Sasaki 等1975证明了双谱 近似正比于指数V，双谱偏差上确界正比于指数J，21212121(),2JW2dd 表4.2列出了MSE最优窗和图4.2所示的乘积窗函数对应的J，V，B值。正如所期望的，MSE最优窗具有最小的E值。Sasaki窗具有最小的J值，这是因为他最初是由使双谱偏差上确界最小化而导出的。另一方面，Parzen窗具有最小的V值。例4.1（见教材）例4.2（见教材）研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 594.3 直接法 由（2.27）或（3.47）多谱的定义提供了一种估计多谱的方法。如果有有限长观察数据，并估计它的傅立叶变换系数，作三重（四重）乘积并平均即得双谱（三谱）的估计。类似于平均周期图法或Welch法。4.3.1 高阶周期图（Higher-Order Periodogram）设()x k，是零均值实平稳时间序列，其DFT为 1,1,0Nk 102()()expNxkFx kjNk ，0,1,1N (4.15)如果()x k没有去均值，则可置0)0(xF。高阶周期图定义为 111111(,)()()()xnnxxnxMFFFN1n (4.16)Brilliges和Rosenblatt指出，对于大的N，有 111(,)(,)xxnnnE MM1n (4.17)并有渐近方差 11112121211var Re(,)var Im(,)()()()xxnnnnxxxnnMMN MMM (4.18)其中2()xM是()x k的真实功率谱，11(,)xnnM是()x k的真实n阶矩谱。显然，高阶周期图估计是渐近无偏的，而其方差正比于功率谱的乘积。进而，当增大时，估计方差也增大。因此，高阶周期图估计是非一致估计（inconsistent）。N通常有两种减小估计方差的途径：(1)在邻域频率上平滑周期图；（频域平滑）(2)平均不相邻时间段上的周期图估计。（Welch法，时域平滑）4.3.2 矩谱的直接估计 设是平稳随机或确定序列的观测值。假设采样周期)()1(Nxx1T，nnN1是研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59高阶谱域沿水平或垂直方向频率样点间所需的间隔。最后，假设高阶谱是在频率00.5（归一化频率）之间估计。直接法如下：步骤1.将数据分为K段，每段含M个样点，即N=K*M，且每段数据去均值。如果为确定性信号，则仅用一个记录(N=M)。如果为了获得FFT算法适合的长度，可对每段数据补零。步骤2.假设(),0,1,.,1ixkkM是第i段数据，计算DFT系数 1()()02()()exp)MiixkFxkjMk，0,1,.,1,1,.,MiK (4.19)步骤3.通常和nNM之间的关系为nnNMM，其中是正奇整数，如。换言之，确定了高阶谱平滑过相邻频率的长度。由于nM21nnMJnMM为偶，而为奇，我们可以折衷的值（最近的整数）。在大小为nMnnN1nnnMMM 个1n的维超矩形窗上进行频域平均估计阶矩谱：n)11(nixF(11)(1)(1(*ixkinnnkFM)1nk()n)1kk)1n,(1)i1111nJJJJkxFnnnnn Ki,2,1 (4.20)当时，在频域不平均，有 1(0)nnMJ*()()()()111111(,)()()(),1,2,iiiinnxxnxnMFFFiK (4.21)步骤4.最后，在K段上平均可得给定数据的阶矩谱，即 nKininnxnMKM111)(11),(1),(4.22)其中，jnj)2(1,2,1nj0,1,1jM。4.3.3 累积量谱估计的Zurbenko方法 根据Zurbenko（1986）的研究，前几节讲述的常规方法（直接和间接类）有一主研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59要的缺点，即，不可能对自变量的所有值获得累积量谱估计。从（2.27）式可见，n阶累积量谱定义在对应于021n的所有点（n为隐含变量）上。然而，直p接法和间接法产生关于主族（其自变量满足关系,1,2,ikn ip p10,1kkn),(213xC和)(2xC先前的常规方法，是一线性非高斯过程的双谱和功率谱，可以获得渐近无偏一致的双谱估计）非 常 不 稳 定 的 估 计。例 如，假 设且。利用),()(32xxCC)，对应于所有,(2131xC和2除021外。从),()(32C获得的功率谱函数将是有偏的。这点对于试图从n阶累积量谱估计重构低阶累积量谱时均成立。对于本例，为了获得一个xCx渐近无偏和一致的Zurbenko方法允许从估计的高阶累积量谱重构渐近无偏和一致的低阶的累积量谱。功率谱估计，我们必须从功率谱估计的开始就重复用常规方法。本节讲述的设(1),(2),()xxx N是一平稳随机过程的采样值，MiKN)1(，Zurbenko步骤1.生成DFT系数 方法步骤如下：,12()()exp()L ML Mk LFx kjMk (4.24.1)和n阶周期图 nkkMisnMisFMI1,21,)(1),(4.24.2)，。1,1,0Ks0i注意到，对于210n，实际的n阶周期图可从),(21,nMisI获得。步骤2.构造n阶矩谱估计 112,1201(,)(,)Knis MsMIKn (4.24.3)步骤3.构造累积量谱估计)()()!1()(),(1121ppnMMpMC (4.24.4)研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59其中，求和是对所有的划分p21，0k，pk,2,1，np1。2.2.2节 给 出 了 当时 的 划 分。Zurbenko指 出，对 于 任 意 的4,3,2p1,021nn，只要NMi1，C12(,)n 是渐近无偏的和一致的估计。4.4 常规方法的复解调制法 Godfrey于1965年提出了基于复解调制的高阶矩谱直接估计方法的等效方法。该方法的基本思想是，一旦计算出复调制，就可容易地生成任意阶的矩谱。假设已经计算出第i段数据的DFT系数，)()(ixF1,1,0M，（4.19）式，则复调制方法步骤如下：步骤 1.加一窄带滤波器，然后将频率移到0点，即构成新序列：其他,0|),()()()(JkkFFixik （4.25.1）和 （4.25.2）其他,02|),()()()(JkkFFixik其中，i=1,2,n-1是（4.20）式中定义的指数。iJJ步骤 2.将其变换到时域生成复解调制。2/2/2expJJkikiSJskjFZ （4.26.2）和 /2/22expJiiSkkJskZFjJ （4.26.2）其中。2JJ步骤 3.用下式估计第i段数据的n阶矩谱 ()121 1()()()11110(,)1 ()()()innJiiissnsnsnMZZZJ 1n 研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59 1,2,.,iK （4.27）步骤 4.如（4.22）式中一样平均得到给定数据的矩谱()()12112111(,)(,)KiinnnniMMK 容易证明（4.27）式方法与（4.20）式的直接法是等效的。4.5 常规方法的统计特性 一般，间接法和直接法得到的估计不同，然而，在特定条件下可以相同。已经证明常规方法得到的是渐近无偏和一致估计。进而，每个高阶谱估计都是渐近复高斯的（例如，对大的M和N）且与在不同频率处的所有其它高阶谱估计不相关。4.5.1 双谱估计的统计特性 假设和分别是一严格平稳零均值随机过程的真实功率和双谱。设)(2xM2(,),(213xM31xM 是由直接法或间接法用长度为N的随机过程的单次实现计算的一致双谱估计。这些方法的主要结果是对于足够大的记录长度和总长度，直接和间接法提供近似无偏估计，即：312312(,)(,)xxE MM （4.28）有渐近方差 3122312312Re(,)1 Im(,)(,)2xxVarMVarM （4.29.1）（直接法）（间接法）)()()()()()(),(2122212232122212232123xxxxxxMMMMKNMMMMKVL，120 （4.29.2）其中，K 是分段数，M 是每段的样本数，V 是双谱窗的总能量，（对于矩形窗 V=1）。窗宽由（4.4）或（4.5）式定义，3L33/(21)NMJ,由（4.20）式定义。为了使当3J,MK 时方差趋于零，应当取 2J+1 为)(KO。从（4.29）式显见，如果用矩研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59形窗，间接法取 L3=N3，则两种常规方法给出近似相同的估计。注意到（4.29）式给出的双谱方差仅是位于主区中频率点),(2123),(21上，而没有给出在其边界频率点上。在主区域边界（如)0,0,12121,2）上双谱方差的确切表示可在 Van Ness1966以及 Lii 和 Rosenblatt1989的文献中找到。Brillinger 和 Rosenblatt1967表明，对于大的 M 和 N，误差双相干 312312312)(,1(nC(4)(,)(,(0,1)xxcMMN )1)n （4.30）是近似复高斯变量，其均值为零，方差为 1。4.5.2高阶谱估计的性能 高阶累积量谱估计的统计特性的推导是非常复杂的，因为一个严平稳零均值过程的四阶或高阶累积量不仅是矩，而有（2.5）或（4.3）给出的形式。尽管如此，Lii 和 Rosenblatt1989证明了用常规方法获得的误差累积量谱估计是具有零均值的渐近联合高斯的。因此，常规方法提供渐近无偏的高阶谱估计：即，n11(,),xxnnE C （4.36）n 阶累积量谱估计的渐近方差可在 Lii 和 Rosenblatt1989中找到。考虑间接三谱估计的渐近方差，并与双谱估计的方差比较，即：32322212243214,xxxxCCCCMKVL21，1230 （4.37）其中 V 是三维窗的总能量，是该过程的真实功率谱。比较（4.37）和（4.29）式可见，主要差别在于和。这表明 L4的值应选的小于 L3以便保证三谱估计合理的方差。实际上，如果选取 L2用于功率谱估计，则为了获得 n 阶累积量谱的平滑估计，Ln的选取应满足)(22xC34L23L112,n=3,4,.nnLL （4.38）从（4.29）和（4.37）式可见，减小常规方法高阶谱估计方差的途径如下：1.增加数据分段数 K；2.减小累积量域窗函数的支撑区域(Ln)的大小或增加频率平滑窗(Jn)的大小；3.增加每段数据的长度 M。研究生课程：现代信号处理高阶统计量分析 课程编号：0211007（博）0221023（硕）西安电子科技大学 授课教师：姬红兵教授 更新日期 2010 年 4 月 1 日 59然而增加数据分段数 K 取决于计算时间并且可能引入潜在的非平稳性。另一方面，在大小为(2Jn+1)的大的多面体上的频域平均或用较小Ln值的累积量窗会降低频率分辨率，且可能增加偏差。在短数据情况下，通过重叠数据段可以增加数据分段数 K。常规高阶谱估计方法具有以下优点：1.由于可以借助 FFT 算法，所以易于实现；2.当观测数据较长时，可以获得较好的估计性能。然而，由于傅立叶变换的不确定性原理(uncertainty principle)，常规方法在多谱域的谐波分辨能力受到限制。因此，在检测间距非常近的两个频率分量间的二次相位耦合时，它的效果较差。此外，对于参数过程，常规方法得到的双谱估计的逼真性也差。作业：作业：参考本章例 4.2 和例 4.2，仿真研究高阶谱估计间接法的性能：（1）研究各种窗函数对于信号双谱估计性能的影响；（2）研究不同数据长度 N、数据分段 K 以及 M 对于双谱估计性能的影响；（3）研究加入高斯白色噪声和有色噪声情况下的信号双谱估计性能；对比分析实验结果，给出结论。