高阶谱理论中的若干问题解析

1、高阶谱理论中的若干问题解析胡念青 作者介绍：胡念青，男，四川师范大学文理学院计科系系主任、副教授。 （四川师范大学文理学院，四川 成都 610110）摘要：本文详细阐述了在信号处理中引入高阶谱（PolySpectrum）理论的重要意义；给出了高阶谱的准确数学定义及相关性质，系统介绍了高阶谱估计的经典算法和现代算法；首次提出了现代算法中的模型阶次判别理论；并分析讨论了非线性系统X2（t）+aX(t)与线性系统h(t)在级联与并联情况下对线性部份的辨识问题；分析讨论了在非高斯假定下、最小相位条件不成立时，过程（模型）传递函数的相位谱估计问题。关健词：高阶谱；信号处理；模型辨识一 引论谱分析在随机过

2、程论、时间序列分析及信号处理中属于一个非常重要的理论问题，也是一个极为有用的处理工具。由维纳创建的功率谱理论，已成为谐波分析、参数估算、信号模型识别、系统辨识及预测控制等多种问题中不可缺少的应用工具。然而，必须看到，与过程的二阶矩相联系的功率谱无论在理论上，还是在应用中都存在着无法避免的局限性。特别是在非高斯及非线性问题的处理中更是如此。因为功率谱仅包含了过程与二阶矩相当的信息量，故只有在高斯情况下，它才能给过程以完整的统计描述。相反在非高斯及非线性情况下，它对过程所能提供的信息描述，便显得极为不够了。由此，前苏联著名的工程数学家Kolmogorov提出了将高阶(大于二阶)矩作付里叶变换这一思

3、想，并进一步发展由Shiryaev提出了高阶谱的概念。之后，由Brillinger初步建立了高阶谱理论，并逐步在海洋波、地震波分析、经济时间序列分析、流体力学及无线电信号处理中找到了广泛的应用。那么到底引用高阶谱理论对实际的工程问题的分析处理有何价值和意义呢？我们试图首先通过对无线电信号处理中的几个问题的分析来对此予以阐释。1. 非线性信号模型分析我们知道，任何具有连续功率谱的平稳随机过程Xt均能找到一个一般的非线性模型表示：Xt=，其中为一个不相关的过程。如果Xt的功率谱密度函数PSD满足这Poley-Winer条件，则上述表达式可以给出单边形式：Xt=显然，是一个不相关的过程( E

4、3;=0 )。这一点仅仅包含了过程的二阶信息。只有在高斯情况下，才能得到Xt的完整统计描述。相反，则可以看到这一信号模型之缺陷。例如在预测问题中，仅依赖这一线性模型，预测精度将得不到保证。于是便相应产生了非线性信号模型。这里我们考虑维纳模型：Volterra级数形式，Ut为输入，Xt为输出：Xt=+引入所谓的广义传递函数： 假定平稳输入过程：Ut=则 Xt=+其中，表示了Ut中w1、w2处的频率分量对Xt中w1+w2处的频率分量之贡献，余类推。当假定Ut=时， 则：Xt=+由此可以看到非线性因素的作用。输出中不再仅仅包含单一的频率分量w0了，而是出现了所谓频率增值现象。因此，我们已不可能仅仅利

5、用二阶特性来确定上述各级传递函数了。为此，必须引进高阶谱。假定Xt中的K阶中心自矩为：C(s1，s2，sk-1)=EXt·Xt+s1Xt+sk-1EXt·EXt+s1EXt+sk-1则可以定义K阶谱：hk(s1，s2，sk-1)= 以同样的方式还可引进两个过程Xt和Yt之间的高阶互谱。下面考虑特殊情况，即Xt的Volterra级数中只包含一个单一项：Xt=且Ut为高斯，则：Xt= Xt与Ut这三阶互矩为：EXt·Ut+s1·Ut+s2 = ·E由于是高斯的(Ut为高斯过程)，则积分号的那四项积之均值可以分解为三项两两相乘之均值的和。即：E

6、83;E+ E·E+ E·E而E=0 w1+w20； E=hu(w1)dw1 w1+w2=0；其中hu(w)为过程Ut之功率谱密度函数。EXt·Ut+s1·Ut+s2 =· + 上式中之第一项恰为EXt·EUt+s1·Ut+s2所以将积分变量w3、w4换为w1、w2，有：Cuux(s1，s2) EXt·Ut+s1·Ut+s2= EXt·EUt+s1·Ut+s2+2由高斯互谱的定义(即它与高阶互矩之关系)，有：huux(w1，w2)=2 于是：由此，我们可以看到引入高阶谱可以估计出单一的

7、各级级数，当Xt的Volterra级数展开中包含多个项式时，可以用维纳分析法将级数展开式重写为一系列正交项之和。其中引入Hermite多项式，这样由新的参数构成的传递函数便可应用上面推导的表达式予以估计。1. 谐波分析中二次相位对的判别问题在谐波分析中存在这样的一种现象，即由于两个谐波分量的相互作用，将可能在差频或和频处产生一个新的频率分量。这种由于某种非线性因素的存在而引起的相位关系，称为二次相位对，一个简单的例子便是幅度调制。显然，二次相位对仅可能出现在谐波相关的地方(即三个频率分量处，一个是另外两个之和或之差)。在某些实际应用中，需要判断功率谱中谐波相关处的峰值能量到底是否由二次相位对所

8、引起。由于功率谱完全忽略了相位关系，因此，它不可能对此问题给出解答。于是我们引入Bispectrum理论。例：过程X(n)= 其中1>2>0、4>5>0、3=1+2、6=4+5，1、2、3、4、5为独立的0，2上均匀分布之随机变量，6=4+5。显然(1、2、3)、 (4、5、6 )均为谐波相关。但只有6处的能谱分量是由二次相位对导致的。即4、5的能谱分量因某种非线性因素的存在而在6=4+5处产生了一个新的能谱分量。而3处的能谱分量却是独立的。该过程的功率谱中包含了6个脉冲分量，从功率谱的谱象中，是无法检验出二次相位对是否存在。考查Xn的三阶自矩R(k,l) EX(n)X

9、(n+k)X(n+l)R(k,l)=1/4cos(4l+5k)+ cos(4l+6k)+ cos(4k+5l)+ cos(6k-5l)+ cos(4k-6l)+ cos(5k-6l)非常重要的，我们发现：三阶自矩中仅仅包含了二次相位对所对应的频率分量。利用其付氏变换，得到Bispectrum在其谱象中，频率(4、5)处将会出现一个峰值，它直接标明了二次相位对出现的地方，而在(1、2)处都不存在这样一个峰值。因此，这个事实给我们提供了一个极为有用的方法，即以Bispectrum作为工具，来判别在谐波相关处是否存在二次相位对。综上所述，从我们对无线电信号处理中几个问题的分析中，我们看到了高阶谱理论

10、中对于功率谱理论所无法解决的问题。诸如非高斯、非线性情况下一些信号处理问题，提供了一个极有价值、极富意义的手段。二 高阶谱的定义及性质定义：假定(k)为一由一类离散或连续K重随机过程X1(t)Xk(t)构成的一个集合，它满足下列三个条件：A. 对于1jk和1h1,hkk， 存在。这里表示K阶互乘矩：EB. 对于1jk，=U为常数，连续情况下- <u< ，离散情况U=0，1，2， 该条件实为一联合平稳假定。C. 对于1jk，在勒贝格测度中，平面w1+w2+wj=0上存在绝对连续的增量(w1+w2+wj) dw1dw2dwj其中(w)为狄拉克delta函数，并满足=·(w1+

11、w2+wj) dw1dw2dwj其中为j阶中心矩。如果一个K重随机过程X1(t)Xk(t)属于上述定义的集合(k)，则 被定义为该过程的K-1阶高阶谱(polyspectrum)。上面给出了高阶谱的一般定义。它实质上为一种互谱，因它涉及了K个过程X1(t)Xk(t)。下面我们将上述定义局限于一个单一的随机过程X(t)，即假定 X(t)X(t)(K个X(t)，属于集合(k)，则有：=·(w1+w2+wk) dw1dw2dwk于是简称为随机过程X(t)的K-1阶谱。上述积分，其积分限为-：离散情况；-：连续情况如果假定一个过程 X(t)X(t)(K个X(t)属于(k)，且满足条件：dt1

12、dt2dtj-1<, 1jk；和： dw1dw2dwj-1<, 1jk；(在离散情况下，上式变为一个和式)则存在付里叶反变换：= dt1dtk或：从上述定义中可以直观地看到，一个过程的高阶谱可以认为是其高阶中心矩的付里叶变换，并在一定条件下，二者构成一变换对。因此高阶谱在许多方面具有付里叶变换谱之性质。具体地说，若Xt是一实过程则：=由Cramer定理知，随机过程Xt(二阶平稳)具有谱表达式：Xtx(t)= 其中心自矩：= E根据两边等式的变量情况，可以知道：w1+w2+wk+1=0由此可以得到高阶谱的谱分解表达式： ·dw1dw2dwk=(w1+w2+wk+1) E由于

13、在目前的实际应用中更多地是遇到二阶谱，即Bispectrum，所以这里再给Bispectrum一个详细的讨论。考虑一个离散过程X(k)，它是三阶平稳且均值为0。有三阶自矩：R(m，n) EX(k) ·X(k+m)·X(k+n)其二维付里叶变换即为所谓的Bispectrum：( w1 、w2 )B(w1，w2)= 可以直接地推导出R(m,n)的对称性：R(m,n)R(n,m)R(-n,m-n)R(n-m,-m)这一性质将在我们后面提出的阶次判别方法中得到应用。即只要确定于R(m,n)在n=0、m=0及m,n0，所界定的无穷三角形区域内的值，便可确定整个R(m,n)的全部值。由

14、R(m,n)的对称性，亦可得到B(w1，w2)之对称性： B(w1，w2)= B(w2，w1)= B*(-w2，-w1)= B(-w1-w2，w2)= B(w1，-w1-w2)B(w1，w2)于w1，w2是以2为周期变化的，由上述对称性知，只要知道了三角形区域w1w2，w20则可以完全确定Bispectrum。下面就几个特殊平稳过程的分析，给出其相应的Bispectrum谱表达式。1. 高斯过程：Xt= x(t)B(w1，w2)dw1dw2=E对于高斯过程，其不同的频率分量是独立的，所以其Bispectrum恒为零。依此类推，其各阶( 2)高阶谱均为零。这是一个非常重要的事实，从中我们可以将高

15、阶谱视为判别一个过程是否为高斯过程之一准则，亦可作为任一过程偏离高斯过程的一个量度。2. 假定Vt是一独立随机变量，EVt=0、EVt2=1、EVt3=又考虑一线性过程Xt=，其中h(t)为一线性装置的冲激响应，h(t)= ，过程的功率谱密度函数为f(w)= 其Bispectrum为：B(w1,w2)= ·H(w1) ·H(w2) ·H*(w1+w2)3. 假定Vt为一高斯过程EVt=0、EVt2=1 ，且其功率谱密度函数为g(w)。又Xt=(Vt) , (t)为一非线性函数。如果为一奇函数(-t) =-(t)，则Xt的Bispectrum将恒为零。现假定(t)中

16、包含了偶函数部份Xt=Vt+Vt2， (a为常数)，则Xt的功率谱密度函数为： f(w)=g(w)+2a2,而其Bispectrum为：B(w1,w2)=2ag(w1)g(w2)+g(w1)g(w1+w2)+g(w2)g(w1+w2)+O(a2) 4. 假定过程Xt=其中h(t)= ，而 T-1、T0、T1为Poission过程随机点出现的时刻ETk+1Tk=。其功率谱f(w)= Bispectrum为：B(w1,w2)= ·H(w1) ·H(w2) ·H*(w1+w2)三 Bispectrum估计中的经典算法现在，我们将注意力集中到谱的估计方法上来，这里我们仅讨

17、论最为常用的估计方法，先考虑所谓经典算法，亦称为“付里叶”型估算法。Xtx(t)= 则：B(1，2)d1d2=EdY(1)dY(2)dY(1+2)由此考虑一零均值过程X(t)(若非均值为0，可以将估计出来的均值从过程中移除)，现接收到Xt中一段数据，对该数据段作有限离散付里叶变换，得到其系数Y，Y可以视为(近视等价) ，继而可将得到的付氏系数进行三重乘积，并在整个数据段上予以平均以消除其中的统计不确定性，便可得到Bispectrum的估计。这一原理与利用周期图加窗平滑方法，对功率谱进行估计是一致的。具体算法如下：假定我们需要估计从0/2到0/2频率范围内的Bispectrum(实际上由对称性知

18、仅需要估计出00/2范围内的值均可)，要求谱线宽度(即谱点的距离为)0。这样要求在时间域上的抽样率至少为每单位时间抽0个样本，且要求不相重叠的每个样本记录的长度至少为N =0/0。现假定共收到K个数据记录，其中每个记录的长度为N=M*N0。这里N的选择应该对FFT运算有帮助，即应为2的幂项，同时为方便起见，我们又须将M选择为奇数，所以我们必须对N0作适当的调整(N0至少为0/0)，即可通过末端加零，或允许各记录重叠的方式来满足FFT所要求的次数，总的数据长度应该是Ntot=K* N0。由Bispectrum的对称性及三角特性021、1+20/2，亦即仅需估计其上的值即可。由此，可进行如下操作A

19、. 若需要，应将估计出来的非零均值从过程中减掉；B. 若需要，通过在末端加零的方式来满足一个有利于FFT运算的长度N，并将改进后的记录记为X0、X1.、XN-1；C. 对该数据记录，作快速付里叶变换FFT：Yq=，0qN/2D. 在频率点处估算功率谱密度函数f()，方法是对M=2L+1个连续的周期图作平均：注意到和q是通过关系式=q(0/N)联系上的。E. 估计Bispectrum：B（1，2）即B（q1·0/2，q2·0/2），这里我们可以加上一些平滑窗，如：或者加上hexagonal窗： F. 最后在k个数据段上进行平均：由此估计出bicoherence：综上所述，Bi

20、spectrum的经典算法可以归纳为：将得到的数据样本划分为一定数量的数据段，在每一数据段上作DFT变换，并对变换结果施以加窗平滑，然后在同一数据段上赋以不同值进行三重乘积并平均，最后再对各样本划分记录施以总体平均，便得到Bispectrum。由上面估算方法的具体步骤可以看出，运算量最大的地方出现在第五步，我们知道：N年复数点的DFT变换，大约需要2Nlog2N次运算，由此可以初步估算出该算法所需要的运算次数为：(M+1/16· N0+log2M+log2 N0)·Ntot 四 Bispectrum估计的现代算法上面我们讨论了Bispectrum估计的经典算法，可以看出其运

21、算量较大，运算过程较为烦琐，而且由于付里叶变换的“不确定性原理”，其在估计谱的分辨率有估真度方面亦受到影响。通过分析利用AR模型来估计功率谱的方法，我们将其沿用到高阶矩，便产生了即将讨论的现代算法。给定一组接收数据样本Xt，我们预告假定Xt是一线性过程，且其背景信号模型为非高斯白噪声(NGWN)激励下的一个AR模型。考虑离散情况：X(n)+ =W(n)其中W(n)为NGWN：EW(n)=0，EW2(n)=Q, EW3(n)=对于m<n，X(m)与W(n)是独立的。假定W(n)是三阶平稳过程，相应地X(n)亦满足三阶平稳性，又假定AR模型是稳定的，亦即AR滤波器的传递函数H(z)=1/A(

22、z)是稳定的，其中A(z)=1+设k、l0，在模型方程两边乘以X(n-k)、X(n-l)，并取均值，得到R(-k,-l)+ 这里R(m,n)为X(n)的三阶矩，(k，l)为二维冲激函数。现考虑m=n线上的2p+1个三阶矩值，它们构成方程R·a=b其中：a1，a1apT，b，00TR 上述矩阵是属于Toplitz型的，但不一定对称。由前面的讨论可知，在这种线性模型的背景下，过程的Bispectrum又可表为: B(w1，w2)= H(w1) ·H(w2) ·H*(w1+w2)，其中： w1 、w2 ，由此可见，欲求B(w1，w2)，仅需估计出H(w)，亦即先将模型识

23、别出来即可，这就需要估计出系数a1ap，以拟合出模型。那么这种拟合成系数的估计，需要在什么样的背景中进行呢？从方程R·a=b可以很方便地给出解答。这种拟合是在二阶矩等价的条件下进行的，即从接收的数据，可估计出 m=0,±1±p处的2p+1个三阶矩点，再利用方程R·a=b，用Toplitz结构算法，解出系数，这样便将模型识别出来了。具体的算法如下：假定接收到N=K·M个数据 X1，X2XN=K·M1. 首先估计三阶自矩序列：A. 将所得样本划分成k个数据段，每段长度为M；B. 在每一个记录段上求得ri(m,n)= i=1kC. 将所得之

24、ri(m,n)在所有记录上平均，得到R(m,n) 的估计； = 2. 利用估计模型阶次 (详见后面)； 3. 取在m=n=0,±1±p上的2+1个值，构成出方程R·a=b，其中 =,00T ，是估计出来的激励噪声的三阶矩，利用Toplitz方程的快速算法，可以解出系数a,即其估计值；4. 利用构造出(w) =1/1 w 可求得Bispectrum的估计值:(w1,w2)(w1) (w1) *(w1+w2)在这里仅需用B(w1，w2)的对称性，估计出其三角形区域W2W10 ,w1+w20 内的值即可。对上述算法，我们再作一些细节的说明：A. 如果接收数据确由AR模型

25、产生，可以证明该算法给出AR模型参数以一致估计；可以看到矩阵 在一般情况下既非对称又非正定，因而没有一个正交矢量空间可以利用；B. 即使模型本身是稳定的，而由估计参数拟合的估计模型不再一定稳定，但因我们最终的问题是估算Bispectrum，而非用来预测和控制，所以稳定性问题显得不重要了；C. 由于模型本身并非高斯型的，所以在最后求解的形式上，最大熵谱法与前面讨论的AR模型法并不等价；D. 我们已经强调，上述模型拟合是在三阶矩序列等价的背景下进行的，因而为Yule-Wolker方程求得的估计系数是不能用来估计Bispectrum。参考文献：1兰华,胡屏.相关与非相关噪声下正弦小信号恢复互高阶谱矩

26、估计方法J.电子测量技术.2001(4).2姚文冰.高阶谱分析在抗噪声语音编码中的应用J.华中科技大学学报(自然科学版)2001,29(9).3马彦,石要武.色噪声背景下正弦频率估计的互高阶谱LS方法.系统工程与电子技术 J. 2005(8).4王威,张宁.高阶谱理论及其在雷达信号处理中的应用J.电子器件1997,20(1).5兰华,石要武.混合噪声背景下正弦参数估计的互高阶谱Piasarenko方法J.计量学报2001,28(1).6M.Ysorer Raghuveer and Chrysostomosl Nikis, Bispectrum Estimation, A Paremetric

27、ApproachJ. IEEE Trans on Assp Vol 33 No.4 1985.7D.R Brllinger, An introduction to polyspetraJ. Ann Math Statist Vol36 pp1351-1374 1995.8K.S.Lii.M.Rosenblantt and C.Van Atta, Bispectral measurements in turbulenceJ. J.Fuid Mech Vol 77 pp45-62 2003.9, Simple effective MA and ARMA techniquesM. Proc ICASSP 2005 PP1426-1429.10M.Huzii, Estimation of coefficients of an autogressive process by using a higher order momentJ. J Time Ser. Analog Vol 2 pp87-93 2005.The Polyspectrum Theoryand Its Application in Singnal Processing Hu nian-qing（Arts and Science Co