现代谱估计计算机仿真实验报告

,.现代谱估计计算机仿真实验报告胡敏在许多工程应用中，利用观测到的一组样本数据估计并分析一个平稳随机信号的功率谱密度是十分重要的。例如，在雷达信号处理中，由回波信号的功率谱密度、谱峰的宽度、高度和位置，可以确定目标的位距离和运动速度；在阵列信号处理中，空间功率谱描述了信号功率随空间角度的分布情况。在许多信号处理应用中，谐波过程经常会遇到，它对应的功率谱为线谱，谐波过程的功率谱估计就是要确定谐波的个数，频率和功率（合称谐波恢复）。为了更好的学习现代信号处理中该部分的内容，我们做了相应的计算机仿真实验。1实验目的1、深入理解现代谱估计和谐波恢复的基本理论，包括ARMA模型，ARMA谱估计，ARMA模型识别，Pisarenko谐波分解，信号子空间和噪声子空间，旋转不变技术(ESPRIT)；2、熟悉与上述谱估计和谐波恢复理论相关的数学方法以及各自的特点，包括最小二乘估计（LS），奇异值分解（SVD），总体最小二乘估计（TLS），特征值分解和广义特征值分解；3、体会ARMA功率谱估计中的Cadzow谱估计子和Kaveh谱估计子，ARMA模型的识别方法，Pisarenko谐波恢复方法，ARMA建模谐波恢复方法，MUSIC方法进行谐波恢复，两种ESPRIT方法（LS-ESPRIT和TLS-ESPRIT进行谐波恢复；2实验原理2.1ARMA谱估计相当多的平稳随机过程都可以通过用白噪声激励线性时不变系统来产生，而线性系统又可以用线性差分方程进行描述，这种差分模型就是自回归—滑动平均（ARMA）模型。而且，任何一个有理式的功率谱密度都可以用一个ARMA随机过程的功率谱密度精确逼近。ARMA随机过程定义为服足下列线性差分方程的离散随机过程)(nx：qjjpiijnebneinxanx11)()()()(（1）式中)(ne是一离散白噪声；式（1）所示的差分方程称为ARMA模型，系统paa,1和qbb,,1分别称为自回归（AR）参数和滑动平均（MA）参数，而p和q分别叫做AR阶数和MA阶数。ARMA过程)(nx可以简记为ARMA),(qp，使用移位算子1z可以把它写作如下形式：)()()()(nezBnxzA（2）其中：ppzazazA111)(,.qqzbzbzB111)(若),0(~)(2Nne，则平稳ARMA),(qp过程的功率谱密度为：222222)()()()()(jjezxeAeBzAzBPj（3）用（3）式进行谱估计必须事先辨识ARMA模型和激励噪声的方差2，而MA参数的估计需要求解非线性方程。为了避开非线性运算，Cadzow和Kaveh分别提出了一种线性谱估计子：1、Cadzow谱估计子jjezezxzAzNzAzNzAzAzBzBP)()()()()()()()()(11112（4）piiiznzN0)((5)piikikan0)(pk,1,0(6)0)(02/)()(kkCkkCkxx(7)其中)(kCx为)(nx的协方差函数。因此用Cadzow谱估计子只需要确定AR阶数和AR参数就能进行ARMA谱估计。2、Kaveh谱估计子jjezqqkkkezxzczAzAzAzAzBzBP)()(1)()()()()(1112（8）pjxjipikjikCaac0)(qk,1,0（9）Kaveh谱估计子需要确定AR阶数，AR参数和MA阶数来进行ARMA谱估计。3、AR阶数和参数的确定对于一个ARMA),(qp过程，可以推出其相关函数满足如下方程：02)()()(ixihihR(10)其中)(nh为系统的冲激响应，根据其定义可以得到：,.nqkkpiibknbinha00)()(（11）由（10）式和（11）式就能推导得到著名的修正Yule-Walker方程，简称MYW方程：0)()(1pixixilRalR，ql（12）若已知AR阶数p，就能通过求解p个MYW方程来求解AR参数：rRa（13）其中：)()2()1()2()()1()1()1()(qRpqRpqRpqRqRqRpqRqRqRxxxxxxxxxRTpaaa,,,21aTxxxpqRqRqR)(,),2(),1(r该方程可以用最小二乘法估计a的值：rRRRaTT1^)(（14）而实际问题中，AR阶数往往是未知的，此时可用奇异值分解法确定AR阶数，总体最小二乘估计AR参数，合称SVD-TLS算法。考虑超定方程：0aRee（15）其中：)()1()()2()1()2()1()()1(eexexexeexexexeexexexepMqRMqRMqRpqRqRqRpqRqRqRRTpppeeaaaaa,,,,,,,1121aepM，ppe，qqe。若pqpqee，就可以通过对eR奇异值分解：HeΣRVU（16）中包含1ep个奇异值，将其归一化：11/kkdefkk，11epk（17）,.选择一个接近于零的数作为阈值，把kk大于此值得最大整数k作为有效秩p，它就是AR阶数。根据总体最小二乘方法可以得到矩阵：Hijpjpniijjjpvv)(1112)(S（18）其中：vij是酉矩阵V第j列的一个加窗段，定义为：Tijjpkvjivjivv)],(,),,1(),,([（19）由)(pS的)(pS可以估计AR参数为：)1,1(/)1,1()()(ppiiaSS，pi,1,0（20）4、MA阶数和参数的确定在AR定阶和参数估计的SVD-TLS算法中，取qqQe，令1QQ，构造超定矩阵：)()1()()1()()1()()1()(2pMQRMQRMQRpQRQRQRpQRQRQRxxxxxxxxxeRpM。计算其SVD，计算比值：)1(1,1)1(1,1)(1,1QppQppQpp（21）式中)(1,1Qpp是e2R对应Q的第1p个奇异值，若大于某个给定的阈值，则接受qQ。在推导Kaveh估计子的过程中可以得到一组非线性方程组，使用Newton-Raphson算法求解该方程组可以得到MA参数，其过程如下：（1）、令MA参数的初始值为：0n，0)0(0cb，0)0(kb，qk,,1,0（2）、计算n迭代的拟合误差函数：kqinkininkcbbf0)()()(，qk,,1,0（22）式中kc,.)(0)(1)(0)()(1)(0)()()(1)()(1)(0)(00nnqnnqnnnqnqnnqnnnbbbbbbbbbbbbF（)()()()1(iinnfFbb1nn2.2Pisarenko谐波分解理论考虑由p个无重复频率的正弦波组成的过程：piiiinfAnx1)2sin()(（24）i是在],[内均匀分布的随机数，则其p个频率由特征多项式：012)12(1211pppzzaza，ipiaa2，pi,,1,0（25）的p对共轭复根决定。一般在加性白噪声),0(~)(2wNnw中观测该过程，所以观测过程)(ny是一个特殊的ARMA)2,2(pp过程，且AR参数和MA参数完全相等。在使用)(nx和)(nw统计独立的假设下，可以得到一个重要的法方程：aaR2wy（26）其中：)0()12()2()12()0()1()2()1()0(yyyyyyyyyyRpRpRpRRRpRRRRTpaaa221,,,,1a这表明2w是观测过程)(ny的自相关矩阵yR的特征值，其对应的特征向量正好是特征多项式（25）的系数，Pisarenko谐波分解法的思想就是找出yR最小的特征值并将其对应的特征向量代入（25）式以求得p对共轭复根，再由下式确定频率：2/)]Re(/)arctan[Im(iiizzf（27）,.2.3ARMA建模法谐波恢复2.2中分析的ARMA)2,2(pp过程不满足MYW方程的条件，但可以推导出其服从的法方程和MYW方程的形式是一致的，所以谐波恢复的ARMA建模算法如下：（1）利用观测数据计算样本自相关矩阵：)()1()()2()1()2()1()()1(MRMpRMpRRpRpRRpRpRyeyeyyeyeyyeyeyeR其中：ppe2，pM2；（2）用SVD-TLS算法确定AR阶数p2和系数向量a的总体最小二乘估计；（3）计算特征多项式（25）的共轭根对（*),iizz，有（27）式确定频率。2.4MUSIC方法考虑白噪声中的p个谐波信号pijnwinwesnxi1)()(，1,,1,0Nn（28）),0(~)(2Nnw，引入以下向量：TNxxxn)]1(,),1(),0([)(xTNjjiieen],,,1[)()1(aTpnnnn)](,),(),([)(ssss21Tpaaa)](,),(),([21A则有：IAPAxxR2)()(HHxxnnE（29）HnnE)()(ssP对xxR进行特征值分解得到：IIUAPAUURU22212)0,0,,,(pHHxxHaadiag（30）式中221,,paa是无加性噪声时信号自相关矩阵的特征值。所以xxR的特征值由p个信号特征值和pN,.个噪声特征值组成，与此对应把特征向量矩阵的列向量分为两部分，即][GSU分别称为S和G分别由信号特征向量和噪声特征向量组成；由S和G分别张成的空间分别叫做信号子空间和噪声子空间。可以推导出：pTpTH,,,,,,)(212100Ga（31）所以用MUSIC方法进行谐波恢复的思想是：以很小的步长对进行搜索，寻找)()(1)(aGGaHHP（32）或)())((1)(aSSIaHHP（33）的p个极大值点，其对应的频率就是所求的谐波频率。2.5旋转不变技术ESPRIT2.4中描述的问题，引入向量：Tmnxnxnxn)]1(,),1(),([)(xTmnwnwnwn)]1(,),1(),([)(wTmnynynyn)]1(,),1(),([)(yTmnxnxnx)](,),2(),1([TNjjiieen],,,1[)()1(aTpnnnn)](,),(),([)(ssss21Tpaaa)](,),(),([21A),,,(21pjjjeeediag于是有：,.)()()(nnnwAsx（34）)1()()(nnnwsAy（35）酉矩阵称为旋转算符，在上面描述的过程)(ny是)(nx平移的结果，可以看作是最简单的旋转。向量)(ny和)(nx的互相关矩阵为：ZAAPyxR2)()(HHHxynnE（36）其中：0100100Z因为平移保持了)(ny和)(nx信号子空间的不变性，所以可以构造矩阵束xyxxCC,：HxxxxxxCAPAIRIR2minHxyxyxyCAAPZRZR2min其中min是xxR的最小特征值。对矩阵束xyxxCC,进行广义特征值分解，其特征值矩阵与矩阵有如下关系：0