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

1、现代谱估计计算机仿真实验报告胡敏在许多工程应用中，利用观测到的一组样本数据估计并分析一个平稳随机信号的功率谱密度是十分重要的。例如，在雷达信号处理中，由回波信号的功率谱密度、谱峰的宽度、高度和位置，可以确定目标的位距离和运动速度；在阵列信号处理中，空间功率谱描述了信号功率随空间角度的分布情况。在许多信号处理应用中，谐波过程经常会遇到，它对应的功率谱为线谱，谐波过程的功率谱估计就是要确定谐波的个数，频率和功率（合称谐波恢复）。为了更好的学习现代信号处理中该部分的内容，我们做了相应的计算机仿真实验。1 实验目的1、深入理解现代谱估计和谐波恢复的基本理论，包括ARMA模型，ARMA谱估计，ARMA模

2、型识别，Pisarenko谐波分解，信号子空间和噪声子空间，旋转不变技术(ESPRIT)；2、熟悉与上述谱估计和谐波恢复理论相关的数学方法以及各自的特点，包括最小二乘估计（LS），奇异值分解（SVD），总体最小二乘估计（TLS），特征值分解和广义特征值分解； 3、体会ARMA功率谱估计中的Cadzow谱估计子和Kaveh谱估计子，ARMA模型的识别方法，Pisarenko谐波恢复方法，ARMA建模谐波恢复方法，MUSIC方法进行谐波恢复，两种ESPRIT方法（LS-ESPRIT和TLS-ESPRIT进行谐波恢复；2 实验原理2.1 ARMA谱估计相当多的平稳随机过程都可以通过用白噪声激励线性时

3、不变系统来产生，而线性系统又可以用线性差分方程进行描述，这种差分模型就是自回归滑动平均（ARMA）模型。而且，任何一个有理式的功率谱密度都可以用一个ARMA随机过程的功率谱密度精确逼近。ARMA随机过程定义为服足下列线性差分方程的离散随机过程： （1）式中是一离散白噪声；式（1）所示的差分方程称为ARMA模型，系统和分别称为自回归（AR）参数和滑动平均（MA）参数，而和分别叫做AR阶数和MA阶数。ARMA过程可以简记为ARMA，使用移位算子可以把它写作如下形式： （2）其中： 若，则平稳ARMA过程的功率谱密度为： （3）用（3）式进行谱估计必须事先辨识ARMA 模型和激励噪声的方差，而MA参

4、数的估计需要求解非线性方程。为了避开非线性运算，Cadzow和Kaveh分别提出了一种线性谱估计子：1、Cadzow谱估计子 （4） (5) (6) (7)其中为的协方差函数。因此用Cadzow谱估计子只需要确定AR阶数和AR参数就能进行ARMA谱估计。2、Kaveh谱估计子 （8） （9）Kaveh谱估计子需要确定AR阶数，AR参数和MA阶数来进行ARMA谱估计。3、AR阶数和参数的确定对于一个ARMA过程，可以推出其相关函数满足如下方程： (10)其中为系统的冲激响应，根据其定义可以得到： （11）由（10）式和（11）式就能推导得到著名的修正Yule-Walker方程，简称MYW方程：，

5、 （12）若已知AR阶数，就能通过求解个MYW方程来求解AR参数： （13）其中：该方程可以用最小二乘法估计的值： （14）而实际问题中，AR阶数往往是未知的，此时可用奇异值分解法确定AR阶数，总体最小二乘估计AR参数，合称SVD-TLS算法。考虑超定方程： （15）其中：，。若，就可以通过对奇异值分解： （16）中包含个奇异值，将其归一化：， （17）选择一个接近于零的数作为阈值，把大于此值得最大整数作为有效秩，它就是AR阶数。根据总体最小二乘方法可以得到矩阵： （18）其中：是酉矩阵第列的一个加窗段，定义为： （19）由的可以估计AR参数为：， （20）4、MA阶数和参数的确定在AR定阶和

6、参数估计的SVD-TLS算法中，取，令，构造超定矩阵：。计算其SVD，计算比值： （21）式中是对应的第个奇异值，若大于某个给定的阈值，则接受。在推导Kaveh估计子的过程中可以得到一组非线性方程组，使用Newton-Raphson算法求解该方程组可以得到MA参数，其过程如下：（1）、令MA参数的初始值为：，（2）、计算迭代的拟合误差函数：， （22）式中可由（9）式求解。（3）计算矩阵：（4）更新MA参数估计向量： （23）（5）判断MA参数估计向量是否收敛，若收敛，则迭代停止，若发散，令返回（2）计算，直到MA参数估计收敛为止。2.2 Pisarenko谐波分解理论考虑由个无重复频率的正弦

7、波组成的过程： （24）是在内均匀分布的随机数，则其个频率由特征多项式：， （25）的对共轭复根决定。一般在加性白噪声中观测该过程，所以观测过程是一个特殊的ARMA过程，且AR参数和MA参数完全相等。在使用和统计独立的假设下，可以得到一个重要的法方程： （26）其中：这表明是观测过程的自相关矩阵的特征值，其对应的特征向量正好是特征多项式（25）的系数，Pisarenko谐波分解法的思想就是找出最小的特征值并将其对应的特征向量代入（25）式以求得对共轭复根，再由下式确定频率： （27）2.3 ARMA建模法谐波恢复2.2中分析的ARMA过程不满足MYW方程的条件，但可以推导出其服从的法方程和MY

8、W方程的形式是一致的，所以谐波恢复的ARMA建模算法如下：（1）利用观测数据计算样本自相关矩阵：其中：，；（2）用SVD-TLS算法确定AR阶数和系数向量的总体最小二乘估计；（3）计算特征多项式（25）的共轭根对（，有（27）式确定频率。2.4 MUSIC方法考虑白噪声中的个谐波信号， （28），引入以下向量：则有： （29）对进行特征值分解得到： （30）式中是无加性噪声时信号自相关矩阵的特征值。所以的特征值由个信号特征值和个噪声特征值组成，与此对应把特征向量矩阵的列向量分为两部分，即分别称为和分别由信号特征向量和噪声特征向量组成；由和分别张成的空间分别叫做信号子空间和噪声子空间。可以推导出

9、： （31）所以用MUSIC方法进行谐波恢复的思想是：以很小的步长对进行搜索，寻找 （32）或 （33）的个极大值点，其对应的频率就是所求的谐波频率。2.5 旋转不变技术ESPRIT2.4中描述的问题，引入向量：于是有： （34） （35）酉矩阵称为旋转算符，在上面描述的过程是平移的结果，可以看作是最简单的旋转。向量和的互相关矩阵为： （36）其中：因为平移保持了和信号子空间的不变性，所以可以构造矩阵束：其中是的最小特征值。对矩阵束进行广义特征值分解，其特征值矩阵与矩阵有如下关系： （37）基本ESPRIT算法（LS-ESPRIT算法）的思想就是用矩阵束的广义特征值矩阵的前个特征值来估计谐波频

10、率，计算公式使用（27）式。TLS-ESPRIT算法的思想是：对进行奇异值分解，确定有效秩，并存储与个主奇异值对应的，和；求的广义特征值分解，得到个广义特征值，它们给出了个谐波频率。3 实验内容仿真的观测数据由下式给出：其中是一列高斯白噪声，。进行如下各项实验：1、 取AR阶数分别为4和6，用最小二乘法估计AR参数，然后使用Cadzow谱估计子进行功率谱估计，并试根据该谱确定谐波频率；2、 假定AR阶数未知，用SVD-TLS方法确定AR阶数和参数，然后使用Cadzow谱估计子进行功率谱估计，并试根据该谱确定谐波频率；3、 用ZHANG方法确定MA阶数，用Kaveh谱估计子进行功率谱估计，用Ne

11、wton-Raphson方法估计MA参数，结合2确定ARMA模型，计算ARMA功率谱，并试根据该谱确定谐波频率；4、 用Pisarenko谐波恢复方法确定谐波频率；5、 用ARMA建模谐波恢复方法确定谐波频率；6、 用MUSIC方法确定谐波频率；7、 用LS-ESPRIT方法确定谐波频率；8、 用TLS-ESPRIT方法确定谐波频率。4 实验过程1、编写基于最小二乘法和Cadzow谱估计子的计算机仿真程序lsestimate.m(代码见附录1)，独立运行程序5次，表1给出了频率的估计数据。从表中数据可以看出第三次程序运行效果比较理想，图1所示的是这次仿真得到的功率谱估计结果。表1 LS法+Ca

12、dzow估计子得到的频率估计数据AR阶数MA阶数第1次第2次第3次第4次第5次450.20310.28910.20310.20310.203160.20310.31250.20310.20310.203170.20310.19530.20310.20310.265680.20310.19530.20310.20310.1953670.20310.20310.20310.25780.203180.19530.19530.20310.20310.148490.20310.19530.20310.20310.2031100.20310.20310.20310.19530.2031取出奇点0.3125

13、，0.2891，0.2578，0.2656，0.1484后，计算数据计算方差和均值得到：均值：m=0.2015方差：s=0.00322、编写用SVD-TLS方法确定AR阶数和参数，ZHANG方法确定MA阶数，Newton-Raphson方法估计MA参数，用Cadzow谱估计子，Kaveh谱估计子，ARMA模型三种方法估计功率谱的程序svdtls.m(代码见附录2)。表2给出了独立运行计算机仿真程序20次的结果。图2给出了程序第二次运行和第十八次运行得到的用三种方法估计的功率谱。去除奇点计算数据均值和方差：均值：m_Cadzow=0.2008, m_ARMA=0.2013, m_Kaveh=0.

14、2013方差：s_Cadzow=0.0032, s_ARMA=0, s_Kaveh=0图1 AR阶数为4和6时，LS+Cadzow估计子得到的功率谱图2 TLS+Cadzow估计子+Kaveh估计子+ARMA得到的功率谱图3 MUSIC方法得到的谱3、编写基于Pisarenko谐波分解，ARMA建模方法、MUSIC方法、ESPRIT方法估计谐波频率的计算机仿真程序sinrecover.m(代码见附录3)，独立运行程序20次，表一给出了每次的频率估计结果。图3两次仿真中用MUSIC方法得到的谱。去除奇点计算数据均值和方差：均值：m_Pisarenko=0.2016,m_ARMA=0.2005,m

15、_MUSIC=0.2031,m_ESPRIT_LS=0.2014, m_ESPRIT_TLS=0.2010方差：s_Pisarenko=0.0021,s_ARMA=0.0005,s_MUSIC=0,s_ESPRIT_LS=0.0006,s_ESPRIT_TLS=0.0005表2、20次运行svdtls.m得到数据次数AR阶数MA阶数AR参数MA参数fCadzowfKavehfARMA164-1.3632, 0.5268, -3.2369-1.4935, -1.4795, -2.365250.8529, 7.4480, 20.3066-0.0332, 1.05640.20130.20130.20

16、132421.0569,0.8716,1.0702,0.882227.4962, 20.0525, 11.70630.20130.20130.20133641.9283, 3.5027, 3.16922.2823, 2.9045, 0.274472.2148, 56.5195, 37.026618.0968, 0.78440.20130.20130.2013458-0.0488, 0.5170, .02400.2315, -0.569616.9689, 0.0247, 2.9020, -3.6511, 3.4873, -4.9357-0.3594, -3.0908, -1.22660.2013

17、0.20130.2013550-0.3929, 0.5429, -0.10880.0126, -0.477817.15980.20130.20130.2013636-1.4370, 1.4764, -0.812229.5453, -23.2459, 13.0521-5.6810,2.7917,-3.2368,3.30010.00780.20130.2013734-1.5510, 1.5472, -0.924528.5816, -22.6882, 13.5667-7.2835, 3.92350.00780.20130.2013847-1.0295, 1.6533-0.6717, 0.409124

18、.4942, -18.1749, 12.5259-4.4965, 1.0124, -0.0920 -0.8637, -0.17410.20130.20130.2013952-1.1455, 1.3884, -0.48590.0420, 0.095124.7588, -17.2112, 11.37930.20130.20130.20131042-0.4675,0.6453,0.2868,-0.228315.6898, -1.0243, 2.59920.20130.20130.201311471.0061, 0.48791.2715, 0.474323.6177,14.2715,8.1385,7.

19、56193.1856,1.2324,1.7278,0.88940.20130.20130.201312620.3367, -0.3751, 0.7865-1.3191, -0.0176, -0.909629.6632, -6.2046, 7.26810.20130.20130.20131357-0.5971, 1.0239, 0.8547-0.4687, 0.859321.7835,-4.8307,6.0957,8.3019-2.8428,4.8565,-2.0569,1.56070.17190.17190.17191442-0.9832,1.6281,-0.6277,0.412926.280

20、0, -19.6938, 16.54430.21090.20130.20131552-0.4934, 1.4335, -0.33450.6037, -0.122221.9715, -10.6891, 13.28520.31250.20130.20131635-0.9209, 1.1604, -0.305921.1532, -12.5666, 4.76712.4555, -3.6615, 1.41380.19530.20130.20131737-1.3890, 1.4491, -0.759225.9700,-18.5269,7.2345,0.3776-1.8058,-0.9267,2.4637,

21、-2.62590.19530.20130.201318651.0014, 2.2738, 1.07552.9498, 0.1103, 1.227661.7779, 32.6979, 48.681523.5966, 27.1928, 13.64960.20130.20130.201319320.5695, 0.2604, 1.144119.4426, 4.0816, 2.06290.20130.20130.201320522.1808, 2.7018, 2.48412.1877, 1.801561.4016, 49.5070, 27.04970.19530.20130.2013表3、20次运行程

22、序sinrecover.m得到的频率估计数据次数fPisa1fPisa2fARMA1fARMA1fMUSICfE_LSfE_TLS10.19650.23000.201100.20310.20140.200520.20100.11520.20030.21380.20310.20240.201930.20480.13920.20090.06280.20310.20150.201040.20310.11310.20130.10170.20310.20130.201150.20100.12690.20050.12680.20310.20250.202160.20080.05230.20040.0840

23、0.20310.20090.200970.200400.200500.20310.20030.200380.19670.23000.20080.13200.20310.20140.200690.200800.200300.20310.20070.2007100.20310.13680.20040.23340.20310.20180.2016110.20100.14040.19980.17090.20310.20120.2011120.20320.14310.20050.13090.20310.20110.2008130.20180.09770.200600.20310.20100.200914

24、0.20070.11050.20050.14280.20310.20200.2015150.20250.15170.20010.14200.20310.20090.2006160.20320.14560.19890.22000.20310.20100.2006170.20370.14580.20060.11350.20310.20170.2014180.20070.16650.20070.15370.20310.20250.2019190.20300.13210.20040.19480.20310.20060.2002200.20360.13040.201000.20310.20080.200

25、35 实验讨论5.1 一般最小二乘估计和总体最小二乘估计的比较实验过程中发现的第一个问题是，LS估计的数值稳定性不如TLS。LS方法估计的结果出现歧点的频率比较高，而TLS方法后频率估计只出现一次歧点。在实验程序中增加三种试验：1、减少噪声的影响（将噪声电平减少十倍），发现这样对LS法的数值稳定性改变很小；2、把样本数量增加一倍后，运行结果出现歧点的次数明显减少。分析原因，是因为样本数量的增加提高自相关矩阵估计的精确度，因为LS法中使用的矩阵理论上是满秩的，但自相关函数计算所使用的估计子是渐进无偏估计，加上噪声的影响，可能导致矩阵亏缺，致使最终结果的不稳定，所以大样本可提高估计准确度；3、增大

26、，估计错误率的提高，使用MYW方程时，选择（即方程组选择的起点）理论上只要大于实际的值就能进行准确估计，但实验结果并非如此，究其原因还是跟自相关函数估计子有关，因为随着数据间隔的增加，实际参加计算的样本数在减少，而自相关函数估计子中使用的样本数是不变的，所以当自相关函数偏离零点越远是，估计值就越小于实际值，而且受噪声的影响就越大。TLS法的合理的原因就在于同时考虑了MYW方程中样本相关矩阵的误差和方程右边样本相关向量的误差（LS法中只认为方程右边样本相关向量含有误差），也就是考虑了总体的误差，实验结果也证明了这么做的合理性。在后面ARMA建模或TLS-ESPRIT谐波恢复方法的应用中也证明了T

27、LS这种合理性，可以看出TLS-ESPRIT估计的频率偏差和方差明显优于LS-ESPRIT得到的结果。5.2 三种功率谱估计方法的比较比较几幅图像，首先看到是Cadzow估计子有明显的能量泄漏（或者称频率泄漏）现象，即频谱上没有出现正弦波应有理想线谱，而在整个频域内都有能量分布，第一个原因当然有噪声的影响，但是从后面的Kaveh估计子和ARMA模型估计的频谱中可以看出，在如此高的信噪比条件（10dB）下噪声的影响是很小的，所以主要原因肯定是由算法导致的：对于Cadzow估计子，从（4）可以看出它为了避开MA参数估计的非线性计算，把噪声方差和MA参数都整合到另外个参数（）中，这样做的缺陷是忽略了

28、MA阶数的影响，在的情况下可以把这个过程理解为有理式的部分分式分解，但如果，这么做显然是不合理的，至少对增加的个数来使（4）式成立。Kaveh估计子在这方面做了改进，至少考虑了MA阶数的影响，从（8）式可以看出Kaveh估计子把噪声方差和MA参数都整合到另外个参数（）中。考虑到实验中的噪声方差是1，这部分的影响也可以忽略，所以其效果可以达到与ARMA模型（考虑AR阶数和系数，MA阶数和系数）相比拟的程度。在实验中改变噪声方差后，就能发现Kaveh 估计子不如ARMA模型理想。5.3 四种谐波恢复方法的比较首先注意到实验中的谐波信号是确定性信号，不能满足平稳随机过程的条件，而且是相关的。所以用平

29、稳随机过程和不相关信号的处理方法来估计信号功率是行不通的，所以实验中的谐波恢复只是估计了谐波的频率。Pisarenko谐波分解法遇到问题是在矩阵的所有特征值中如何确定哪个是噪声方差对应的特征值，实验程序中的方法是选最小的那个特征值，这样做显然是不合理的，从信号子空间和噪声子空间的理论分析中可知，对于实验中的信号，谐波信号对应的特征值只有四个，如果事先不知道谐波个数的话，就必须构造超定的矩阵来进行处理，其对应有多个最小的特征值，所以就必须逐步进行降维处理，直到只有一个最小的特征值，这个过程的实现是相当困难的，因为只有一个最小的判定条件根本就无法量化，但实验程序中假设谐波个数已知的条件下，发现其估

30、计准确率还是很高的。ARMA建模法可以取得较好的效果，得益于SVD较好估计了谐波的个数，TLS能较好估计特征多项式的系数向量。MUSIC方法和ESPRIT的方法能准确地估计谐波的频率，但是MUSIC方法对噪声处理能力不如ESPRIT方法，这点在5.5中讨论。5.4 关于频率分辨率提高第二个谐波的幅度跟第一个谐波一致，进行实验，发现从功率谱上还是不能发现第二个谐波，所以谱估计对频率的分辨率是有限的，这个限度一方面跟谱估计子有关，一方面跟计算谱值的点数有关，点数增加可以提高分辨率。将第二频率提高到0.25以上，就能实现较好的分离，图四给出了提高频率间隔前后的谱估计图像。相比之下，几种谐波恢复方法（

31、除了Pisarenko方法）的频率分辨率都比较高，特别是在增加样本数据点的情况下，除了Pisarenko方法经常出现歧点（我觉得原因开自方法本身的局限性）外，其它几种方法都能准确地进行频率估计。图3 提高频率间隔前后的功率谱估计图5.5 关于信噪比提高另一个谐波的频率到0.313，观察信噪比对上述方法的影响。实际上，模拟观测信号中，第二个谐波与噪声的信噪比是0dB，但发现三种谱估计方法有较强的去噪声能力，图4中给出了三种谱估计方法在信噪比0dB的条件下发现信号的情景。相比之下，四种谐波恢复方法噪声背景的处理能力比较弱。特别是MUSIC方法，只有当信噪比提高到5dB以上才能进行正确的频率分离。图

32、5给出了提高信噪比到10dB后，MUSIC方法得到的谱。图4 谱估计方法在信噪比0dB下发现信号 图5 MUSIC方法在提高信噪比后发现信号5.6 关于几个阈值的选取程序编写中发现阈值的取法并没有固定的标准，教材上给出阈值的例子在实验中计算效果并不很理想，所以阈值应该针对特点的问题来选取，本实验中确定阈值的方法是：先运行几次程序，观察实验数据，选择一个可行的数作为阈值，然后经过多次尝试调整之后取定的。附录1、程序lsestimate.mclc,clear;t=1:128;%数据时间向量N=length(t);%数据个数f=t/N;%频率向量w=2*pi*f;%角频率向量z=exp(j*w);r

33、andn('state',sum(100*clock);%每次计算给随机数产生设置不同的起点wn=randn(size(t);%功率为1的高斯白噪声x=sqrt(20)*sin(2*pi*0.2*t)+sqrt(2)*sin(2*pi*0.213*t)+wn;%观测数据%估计自相关函数R=xcorr(x);%估计协方差函数COR=xcov(x);CORP=COR;CORP(N)=COR(N)/2;%用一般最小二乘法估计功率谱和谐波频率pls=4;%AR模型的阶数取4for qels=5:8;%选取修正Yule-Walker方程的起点 %以下求相关函数估计值组成的Hankel矩阵

34、 RLS4=R(qels:(qels+pls-1); for k=1:(pls-1) RLS4=R(qels-k):(qels-k+pls-1);RLS4; end bls4=-R(qels+1):(qels+pls);%Yule-Walker方程右边的值 als4(:,(qels-pls)=(RLS4'*RLS4)(-1)*RLS4'*(bls4)'%最小二乘估计公式enddisplay('AR阶数为4，qe取5,6,7,8时，用一般最小二乘法估计的AR参数:')als4=als4(pls:-1:1,:)%采用Cadzow谱估计子进行功率谱估计als4

35、=ones(1,4);als4;als4=als4'CCadzowp4=CORP(N:N+pls);for i=1:pls CCadzowp4=CCadzowp4;CORP(N-i:N-i+pls);endnCadzowp4=als4*CCadzowp4;Nzp4=zeros(4,N);Azp4=zeros(4,N);Nnzp4=zeros(4,N);Anzp4=zeros(4,N);for k=0:pls Nzp4=Nzp4+nCadzowp4(:,k+1)*(z.(-k); Azp4=Azp4+als4(:,k+1)*(z.(-k); Nnzp4=Nnzp4+nCadzowp4(:

36、,k+1)*(z.k); Anzp4=Anzp4+als4(:,k+1)*(z.k);endPWxCadzowp4=Nzp4./Azp4+Nnzp4./Anzp4;PWxCadzowp4=abs(PWxCadzowp4);figure(1);for i=1:4 subplot(4,1,i); stem(f,PWxCadzowp4(i,:),'filled'); flsp4(i)=find(PWxCadzowp4(i,1:(N/2-1)=max(PWxCadzowp4(i,1:(N/2-1)/N;enddisplay('AR阶数为6，qe取7,6,7,8时，用Cadzow

37、谱估计子估计频率:')flsp4pls=6;%AR模型的阶数取6for qels=7:10;%选取修正Yule-Walker方程的起点 %以下求相关函数估计值组成的Hankel矩阵 RLS6=R(qels:(qels+pls-1); for k=1:(pls-1) RLS6=R(qels-k):(qels-k+pls-1);RLS6; end bls6=-R(qels+1):(qels+pls);%Yule-Walker方程右边的值 als6(:,(qels-pls)=(RLS6'*RLS6)(-1)*RLS6'*(bls6)'%最小二乘估计公式enddispl

38、ay('AR阶数为6，qe取7,8，9，10时，用一般最小二乘法估计的AR参数:')als6=als6(pls:-1:1,:)%采用Cadzow谱估计子进行功率谱估计als6=ones(1,4);als6;als6=als6'CCadzowp6=CORP(N:N+pls);for i=1:pls CCadzowp6=CCadzowp6;CORP(N-i:N-i+pls);endnCadzowp6=als6*CCadzowp6;Nzp6=zeros(4,N);Azp6=zeros(4,N);Nnzp6=zeros(4,N);Anzp6=zeros(4,N);for k=0

39、:pls Nzp6=Nzp6+nCadzowp6(:,k+1)*(z.(-k); Azp6=Azp6+als6(:,k+1)*(z.(-k); Nnzp6=Nnzp6+nCadzowp6(:,k+1)*(z.k); Anzp6=Anzp6+als6(:,k+1)*(z.k);endPWxCadzowp6=Nzp6./Azp6+Nnzp6./Anzp6;PWxCadzowp6=abs(PWxCadzowp6);figure(2);for i=1:4 subplot(4,1,i); stem(f,PWxCadzowp6(i,:),'filled'); flsp6(i)=find(P

40、WxCadzowp6(i,1:(N/2-1)=max(PWxCadzowp6(i,1:(N/2-1)/N;enddisplay('AR阶数为6，qe取7,8，9，10时，用Cadzow谱估计子估计频率:')flsp62、程序svdtls.mclc,clear;t=1:128;%数据时间向量N=length(t);%数据个数f=t/N;%频率向量w=2*pi*f;%角频率向量z=exp(j*w);randn('state',sum(100*clock);%每次计算给随机数产生设置不同的起点wn=randn(size(t);%功率为1的高斯白噪声x=sqrt(20)

41、*sin(2*pi*0.2*t)+sqrt(2)*sin(2*pi*0.213*t)+wn;%观测数据%估计自相关函数R=xcorr(x);%估计协方差函数COR=xcov(x);CORP=COR;CORP(N)=COR(N)/2;%用SVD-TLS法估计功率谱和谐波频率peTLS=6;qeTLS=10;%选取修正Yule-Walker方程的起点%以下求增广矩阵RTLS=R(N+qeTLS+1):-1:(N+qeTLS+1-peTLS);for k=2:(peTLS+2) RTLS=RTLS;R(N+qeTLS+k):-1:(N+qeTLS+k-peTLS);end%以下对RTLS做奇异值分解

42、SVDUTLS STLS VTLS=svd(RTLS);%奇异值归一化，并确定有效秩pSTLS1=STLS/STLS(1,1);pTLS=0;k=peTLS;while pTLS=0 if STLS1(k,k)>=0.005 pTLS=k; else k=k-1; endenddisplay('SVD-TLS法确定的AR阶数为：')pTLS%计算SpSpTLS=zeros(pTLS+1,pTLS+1);for i=1:pTLS for k=1:(peTLS+1-pTLS) SpTLS=SpTLS+STLS(i,i)*VTLS(k:k+pTLS,i)*(VTLS(k:k+p

43、TLS,i)' endend%AR参数的估计值为：SpnTLS=SpTLS(-1);aTLS=SpnTLS(:,1)/SpnTLS(1,1);display('SVD-TLS法确定的AR参数为：')aTLS%采用Cadzow谱估计子进行功率谱估计CCadzow=CORP(N:N+pTLS);for i=1:pTLS CCadzow=CCadzow;CORP(N-i:N-i+pTLS);endnCadzow=aTLS'*CCadzow;Nz=zeros(1,N);Az=zeros(1,N);Nnz=zeros(1,N);Anz=zeros(1,N);for k=0

44、:pTLS Nz=Nz+nCadzow(k+1)*(z.(-k); Az=Az+aTLS(k+1)*(z.(-k); Nnz=Nnz+nCadzow(k+1)*(z.k); Anz=Anz+aTLS(k+1)*(z.k);endPWxCadzow=Nz./Az+Nnz./Anz;PWxCadzow=abs(PWxCadzow);subplot(3,1,1);stem(f,PWxCadzow,'filled');title('Cadzow谱估计子估计的功率谱');%通过功率谱估计频率display('Cadzow谱估计子估计的频率：')fCadzo

45、w=find(PWxCadzow(1:(N/2-1)=max(PWxCadzow(1:(N/2-1)/N%MA阶数和参数的辨识Q=qeTLS;%Q=qe>q%以下求增广矩阵RTLStemp=R(N+Q):-1:(N+Q-pTLS);for i=1:(pTLS+5) % 使RTLStemp为超定矩阵 RTLStemp=RTLStemp;R(N+Q+i):-1:(N+Q+i-pTLS);endUTLStemp STLStemp VTLStemp=svd(RTLStemp);OQ1=STLStemp(pTLS+1,pTLS+1);flag=0;while flag<=0.3 %以下求增广

46、矩阵 Q=Q-1; RTLStemp=R(N+Q):-1:(N+Q-pTLS); for i=1:(pTLS+5) % 使RTLStemp为超定矩阵 RTLStemp=RTLStemp;R(N+Q+i):-1:(N+Q+i-pTLS); end UTLStemp STLStemp VTLStemp=svd(RTLStemp); OQ2=STLStemp(pTLS+1,pTLS+1); flag=abs(OQ1-OQ2)/OQ2); qTLS=Q;enddisplay('zhang方法确定的MA阶数为：')qTLS%Kaveh估计子系数和Newton-Raphson算法估计MA系

47、数ck=zeros(qTLS+1,1);for k=1:(qTLS+1) for m=1:pTLS+1 for n=1:pTLS+1 ck(k)=ck(k)+aTLS(m)*conj(aTLS(n)*COR(N+k-1+n-m); end endend%Newton-Raphson算法bNewton=zeros(qTLS+1,1);%以下是MA系数的初始值biNewton=ck(1)(1/2);zeros(qTLS,1);while min(abs(bNewton)=0%以下是第i次的拟合误差函数I=find(abs(bNewton)=0);biNewton(I)=0;biNewton=biN

48、ewton+bNewton;fki=-ck;for k=1:qTLS+1 for m=1:(qTLS+2-k) fki(k)=fki(k)+biNewton(m)*biNewton(m+k-1); endendFi=hankel(biNewton)+rot90(hankel(biNewton(qTLS+1:-1:1),3);%以下是第i+1次跌代的MA系数bi1=biNewton-Fi(-1)*fki;bNewton=bi1;%判断第i+1次跌代的MA系数的收敛性for k=1:qTLS+1 if abs(bi1(k)-biNewton(k)/bi1(k)<=0.05 bNewton(k

49、)=bi1(k); endendbiNewton=bi1;enddisplay('Newton-Raphson算法确定的MA参数为：')bNewton%用ARMA模型和Kaveh估计子进行功率谱估计Bz=zeros(1,N);Cz=zeros(1,N);Bnz=zeros(1,N);Cnz=zeros(1,N);for k=0:qTLS Bz=Bz+bNewton(k+1)*(z.(-k); Cz=Cz+ck(k+1)*(z.(-k); Bnz=Bnz+bNewton(k+1)*(z.k); Cnz=Cnz+ck(k+1)*(z.k);endCz(1)=Cz(1)/2;Cnz(

50、1)=Cnz(1)/2;%ARMA模型功率谱估计PWxzARMA=(Bz.*Bnz)./(Az.*Anz);PWx=abs(PWxzARMA);subplot(3,1,2);stem(f,PWx,'filled');title('ARMA谱估计的功率谱');%通过功率谱估计频率display('ARMA谱估计的频率：')fARMAP=find(PWx(1:(N/2-1)=max(PWx(1:(N/2-1)/N%Kaveh估计子功率谱估计PWxzKaveh=(Cz+Cnz)./(Az.2);PWx=abs(PWxzKaveh);subplot(3,1,3);stem(f,PWx,'filled');title('Kaveh谱估计子估计的功率谱');%通过功率谱估计频率display(&#