现代谱估计计算机仿真实验报告.doc
- 文档编号:566688
- 上传时间:2023-04-29
- 格式:DOC
- 页数:21
- 大小:899.50KB
现代谱估计计算机仿真实验报告.doc
《现代谱估计计算机仿真实验报告.doc》由会员分享,可在线阅读,更多相关《现代谱估计计算机仿真实验报告.doc(21页珍藏版)》请在冰点文库上搜索。
现代谱估计计算机仿真实验报告
胡敏
在许多工程应用中,利用观测到的一组样本数据估计并分析一个平稳随机信号的功率谱密度是十分重要的。
例如,在雷达信号处理中,由回波信号的功率谱密度、谱峰的宽度、高度和位置,可以确定目标的位距离和运动速度;在阵列信号处理中,空间功率谱描述了信号功率随空间角度的分布情况。
在许多信号处理应用中,谐波过程经常会遇到,它对应的功率谱为线谱,谐波过程的功率谱估计就是要确定谐波的个数,频率和功率(合称谐波恢复)。
为了更好的学习现代信号处理中该部分的内容,我们做了相应的计算机仿真实验。
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随机过程定义为服足下列线性差分方程的离散随机过程:
(1)
式中是一离散白噪声;式
(1)所示的差分方程称为ARMA模型,系统和分别称为自回归(AR)参数和滑动平均(MA)参数,而和分别叫做AR阶数和MA阶数。
ARMA过程可以简记为ARMA,使用移位算子可以把它写作如下形式:
(2)
其中:
若,则平稳ARMA过程的功率谱密度为:
(3)
用(3)式进行谱估计必须事先辨识ARMA模型和激励噪声的方差,而MA参数的估计需要求解非线性方程。
为了避开非线性运算,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方程:
,(12)
若已知AR阶数,就能通过求解个MYW方程来求解AR参数:
(13)
其中:
该方程可以用最小二乘法估计的值:
(14)
而实际问题中,AR阶数往往是未知的,此时可用奇异值分解法确定AR阶数,总体最小二乘估计AR参数,合称SVD-TLS算法。
考虑超定方程:
(15)
其中:
,,。
若,就可以通过对奇异值分解:
(16)
中包含个奇异值,将其归一化:
,(17)
选择一个接近于零的数作为阈值,把大于此值得最大整数作为有效秩,它就是AR阶数。
根据总体最小二乘方法可以得到矩阵:
(18)
其中:
是酉矩阵第列的一个加窗段,定义为:
(19)
由的可以估计AR参数为:
,(20)
4、MA阶数和参数的确定
在AR定阶和参数估计的SVD-TLS算法中,取,令,构造超定矩阵:
。
计算其SVD,计算比值:
(21)
式中是对应的第个奇异值,若大于某个给定的阈值,则接受。
在推导Kaveh估计子的过程中可以得到一组非线性方程组,使用Newton-Raphson算法求解该方程组可以得到MA参数,其过程如下:
(1)、令MA参数的初始值为:
,,,
(2)、计算迭代的拟合误差函数:
,(22)
式中可由(9)式求解。
(3)计算矩阵:
(4)更新MA参数估计向量:
(23)
(5)判断MA参数估计向量是否收敛,若收敛,则迭代停止,若发散,令返回
(2)计算,直到MA参数估计收敛为止。
2.2Pisarenko谐波分解理论
考虑由个无重复频率的正弦波组成的过程:
(24)
是在内均匀分布的随机数,则其个频率由特征多项式:
,,(25)
的对共轭复根决定。
一般在加性白噪声中观测该过程,所以观测过程是一个特殊的ARMA过程,且AR参数和MA参数完全相等。
在使用和统计独立的假设下,可以得到一个重要的法方程:
(26)
其中:
这表明是观测过程的自相关矩阵的特征值,其对应的特征向量正好是特征多项式(25)的系数,Pisarenko谐波分解法的思想就是找出最小的特征值并将其对应的特征向量代入(25)式以求得对共轭复根,再由下式确定频率:
(27)
2.3ARMA建模法谐波恢复
2.2中分析的ARMA过程不满足MYW方程的条件,但可以推导出其服从的法方程和MYW方程的形式是一致的,所以谐波恢复的ARMA建模算法如下:
(1)利用观测数据计算样本自相关矩阵:
其中:
,;
(2)用SVD-TLS算法确定AR阶数和系数向量的总体最小二乘估计;
(3)计算特征多项式(25)的共轭根对(,有(27)式确定频率。
2.4MUSIC方法
考虑白噪声中的个谐波信号
,(28)
,引入以下向量:
则有:
(29)
对进行特征值分解得到:
(30)
式中是无加性噪声时信号自相关矩阵的特征值。
所以的特征值由个信号特征值和个噪声特征值组成,与此对应把特征向量矩阵的列向量分为两部分,即
分别称为和分别由信号特征向量和噪声特征向量组成;由和分别张成的空间分别叫做信号子空间和噪声子空间。
可以推导出:
(31)
所以用MUSIC方法进行谐波恢复的思想是:
以很小的步长对进行搜索,寻找
(32)
或(33)
的个极大值点,其对应的频率就是所求的谐波频率。
2.5旋转不变技术ESPRIT
2.4中描述的问题,引入向量:
于是有:
(34)
(35)
酉矩阵称为旋转算符,在上面描述的过程是平移的结果,可以看作是最简单的旋转。
向量和的互相关矩阵为:
(36)
其中:
因为平移保持了和信号子空间的不变性,所以可以构造矩阵束:
其中是的最小特征值。
对矩阵束进行广义特征值分解,其特征值矩阵与矩阵有如下关系:
(37)
基本ESPRIT算法(LS-ESPRIT算法)的思想就是用矩阵束的广义特征值矩阵的前个特征值来估计谐波频率,计算公式使用(27)式。
TLS-ESPRIT算法的思想是:
对进行奇异值分解,确定有效秩,并存储与个主奇异值对应的,和;求的广义特征值分解,得到个广义特征值,它们给出了个谐波频率。
3实验内容
仿真的观测数据由下式给出:
其中是一列高斯白噪声,。
进行如下各项实验:
1、取AR阶数分别为4和6,用最小二乘法估计AR参数,然后使用Cadzow谱估计子进行功率谱估计,并试根据该谱确定谐波频率;
2、假定AR阶数未知,用SVD-TLS方法确定AR阶数和参数,然后使用Cadzow谱估计子进行功率谱估计,并试根据该谱确定谐波频率;
3、用ZHANG方法确定MA阶数,用Kaveh谱估计子进行功率谱估计,用Newton-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所示的是这次仿真得到的功率谱估计结果。
表1LS法+Cadzow估计子得到的频率估计数据
AR阶数
MA阶数
第1次
第2次
第3次
第4次
第5次
4
5
0.2031
0.2891
0.2031
0.2031
0.2031
6
0.2031
0.3125
0.2031
0.2031
0.2031
7
0.2031
0.1953
0.2031
0.2031
0.2656
8
0.2031
0.1953
0.2031
0.2031
0.1953
6
7
0.2031
0.2031
0.2031
0.2578
0.2031
8
0.1953
0.1953
0.2031
0.2031
0.1484
9
0.2031
0.1953
0.2031
0.2031
0.2031
10
0.2031
0.2031
0.2031
0.1953
0.2031
取出奇点0.3125,0.2891,0.2578,0.2656,0.1484后,计算数据计算方差和均值得到:
均值:
m=0.2015
方差:
s=0.0032
2、编写用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.2013
方差:
s_Cadzow=0.0032,s_ARMA=0,s_Kaveh=0
图1AR阶数为4和6时,LS+Cadzow估计子得到的功率谱
图2TLS+Cadzow估计子+Kaveh估计子+ARMA得到的功率谱
图3MUSIC方法得到的谱
3、编写基于Pisarenko谐波分解,ARMA建模方法、MUSIC方法、ESPRIT方法估计谐波频率的计算机仿真程序sinrecover.m(代码见附录3),独立运行程序20次,表一给出了每次的频率估计结果。
图3两次仿真中用MUSIC方法得到的谱。
去除奇点计算数据均值和方差:
均值:
m_Pisarenko=0.2016,m_ARMA=0.2005,m_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参数
f
Cadzow
f
Kaveh
f
ARMA
1
6
4
-1.3632,0.5268,-3.2369
-1.4935,-1.4795,-2.3652
50.8529,7.4480,20.3066
-0.0332,1.0564
0.2013
0.2013
0.2013
2
4
2
1.0569,0.8716,1.0702,0.8822
27.4962,20.0525,11.7063
0.2013
0.2013
0.2013
3
6
4
1.9283,3.5027,3.1692
2.2823,2.9045,0.2744
72.2148,56.5195,37.0266
18.0968,0.7844
0.2013
0.2013
0.2013
4
5
8
-0.0488,0.5170,.0240
0.2315,-0.5696
16.9689,0.0247,2.9020,-3.6511,3.4873,-4.9357
-0.3594,-3.0908,-1.2266
0.2013
0.2013
0.2013
5
5
0
-0.3929,0.5429,-0.1088
0.0126,-0.4778
17.1598
0.2013
0.2013
0.2013
6
3
6
-1.4370,1.4764,-0.8122
29.5453,-23.2459,13.0521
-5.6810,2.7917,-3.2368,3.3001
0.0078
0.2013
0.2013
7
3
4
-1.5510,1.5472,-0.9245
28.5816,-22.6882,13.5667
-7.2835,3.9235
0.0078
0.2013
0.2013
8
4
7
-1.0295,1.6533
-0.6717,0.4091
24.4942,-18.1749,12.5259
-4.4965,1.0124,-0.0920
-0.8637,-0.1741
0.2013
0.2013
0.2013
9
5
2
-1.1455,1.3884,-0.4859
0.0420,0.0951
24.7588,-17.2112,11.3793
0.2013
0.2013
0.2013
10
4
2
-0.4675,0.6453,0.2868,-0.2283
15.6898,-1.0243,2.5992
0.2013
0.2013
0.2013
11
4
7
1.0061,0.4879
1.2715,0.4743
23.6177,14.2715,8.1385,7.5619
3.1856,1.2324,1.7278,0.8894
0.2013
0.2013
0.2013
12
6
2
0.3367,-0.3751,0.7865
-1.3191,-0.0176,-0.9096
29.6632,-6.2046,7.2681
0.2013
0.2013
0.2013
13
5
7
-0.5971,1.0239,0.8547
-0.4687,0.8593
21.7835,-4.8307,6.0957,8.3019
-2.8428,4.8565,-2.0569,1.5607
0.1719
0.1719
0.1719
14
4
2
-0.9832,1.6281,-0.6277,0.4129
26.2800,-19.6938,16.5443
0.2109
0.2013
0.2013
15
5
2
-0.4934,1.4335,-0.3345
0.6037,-0.1222
21.9715,-10.6891,13.2852
0.3125
0.2013
0.2013
16
3
5
-0.9209,1.1604,-0.3059
21.1532,-12.5666,4.7671
2.4555,-3.6615,1.4138
0.1953
0.2013
0.2013
17
3
7
-1.3890,1.4491,-0.7592
25.9700,-18.5269,7.2345,0.3776
-1.8058,-0.9267,2.4637,-2.6259
0.1953
0.2013
0.2013
18
6
5
1.0014,2.2738,1.0755
2.9498,0.1103,1.2276
61.7779,32.6979,48.6815
23.5966,27.1928,13.6496
0.2013
0.2013
0.2013
19
3
2
0.5695,0.2604,1.1441
19.4426,4.0816,2.0629
0.2013
0.2013
0.2013
20
5
2
2.1808,2.7018,2.4841
2.1877,1.8015
61.4016,49.5070,27.0497
0.1953
0.2013
0.2013
表3、20次运行程序sinrecover.m得到的频率估计数据
次数
fPisa1
fPisa2
fARMA1
fARMA1
fMUSIC
fE_LS
fE_TLS
1
0.1965
0.2300
0.2011
0
0.2031
0.2014
0.2005
2
0.2010
0.1152
0.2003
0.2138
0.2031
0.2024
0.2019
3
0.2048
0.1392
0.2009
0.0628
0.2031
0.2015
0.2010
4
0.2031
0.1131
0.2013
0.1017
0.2031
0.2013
0.2011
5
0.2010
0.1269
0.2005
0.1268
0.2031
0.2025
0.2021
6
0.2008
0.0523
0.2004
0.0840
0.2031
0.2009
0.2009
7
0.2004
0
0.2005
0
0.2031
0.2003
0.2003
8
0.1967
0.2300
0.2008
0.1320
0.2031
0.2014
0.2006
9
0.2008
0
0.2003
0
0.2031
0.2007
0.2007
10
0.2031
0.1368
0.2004
0.2334
0.2031
0.2018
0.2016
11
0.2010
0.1404
0.1998
0.1709
0.2031
0.2012
0.2011
12
0.2032
0.1431
0.2005
0.1309
0.2031
0.2011
0.2008
13
0.2018
0.0977
0.2006
0
0.2031
0.2010
0.2009
14
0.2007
0.1105
0.2005
0.1428
0.2031
0.2020
0.2015
15
0.2025
0.1517
0.2001
0.1420
0.2031
0.2009
0.2006
16
0.2032
0.1456
0.1989
0.2200
0.2031
0.2010
0.2006
17
0.2037
0.1458
0.2006
0.1135
0.2031
0.2017
0.2014
18
0.2007
0.1665
0.2007
0.1537
0.2031
0.2025
0.2019
19
0.2030
0.1321
0.2004
0.1948
0.2031
0.2006
0.2002
20
0.2036
0.1304
0.2010
0
0.2031
0.2008
0.2003
5实验讨论
5.1一般最小二乘估计和总体最小二乘估计的比较
实验过程中发现的第一个问题是,LS估计的数值稳定性不如TLS。
LS方法估计的结果出现歧点的频率比较高,而TLS方法后频率估计只出现一次歧点。
在实验程序中增加三种试验:
1、减少噪声的影响(将噪声电平减少十倍),发现这样对LS法的数值稳定性改变很小;2、把样本数量增加一倍后,运行结果出现歧点的次数明显减少。
分析原因,是因为样本数量的增加提高自相关矩阵估计的精确度,因为LS法中使用的矩阵理论上是满秩的,但自相关函数计算所使用的估计子是渐进无偏估计,加上噪声的影响,可能导致矩阵亏缺,致使最终结果的不稳定,所以大样本可提高估计准确度;3、增大,估计错误率的提高,使用MYW方程时,选择(即方程组选择的起点)理论上只要大于实际的值就能进行准确估计,但实验结果并非如此,究其原因还是跟自相关函数估计子有关,因为随着数据间隔的增加,实际参加计算的样本数在减少,而自相关函数估计子中使用的样本数是不变的,所以当自相关函数偏离零点越远是,估计值就越小于实际值,而且受噪声的影响就越大。
TLS法的合理的原因就在于同时考虑了MYW方程中样本相关矩阵的误差和方程右边样本相关向量的误差(LS法中只认为方程右边样本相关向量含有误差),也就是考虑了总体的误差,实验结果也证明了这么做的合理性。
在后面ARMA建模或TLS-ESPRIT谐波恢复方法
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 估计 计算机仿真 实验 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)