现代信号处理仿真作业一（3.18谐波恢复）Word下载.docx
《现代信号处理仿真作业一（3.18谐波恢复）Word下载.docx》由会员分享，可在线阅读，更多相关《现代信号处理仿真作业一（3.18谐波恢复）Word下载.docx（16页珍藏版）》请在冰点文库上搜索。
y x w x w
将（3）式代入
（2）式得：
Ry（k）+ aiRy（k-i）=0,k>
2p （4）
这就是加性白噪声中观测到的p个正弦信号所组成的谐波信号的ARMA
过程的所服从的法方程。
3.谐波恢复方法
利用以上两种理论所推导的结果，可以通过以下方法进行本题的仿真实验（程序见附页）：
（1.） 根据题目要求给出仿真数据
x（n）=
20sin（2p0.2n）+
2sin（2p0.213n）+w（n）
其中w（n）为0均值，单位方差的高斯白噪声。
（2.） 求x（n）的自相关矩阵Rxx
（3.） 利用上面推导出来的（4）式给出该谐波过程的ARMA过程所服从的法方程。
其中公式中的y对应题目中的x。
（4.）若用最小二乘方法，则根据题目要求，可直接取AR阶数2p分别为4和6，然后用最小二乘法求解法方程，得出AR参数的估计值a1,a2,K,ap；
若用总体最小二乘方法，则取较大的M和Pe，列出法方程的增广矩阵，求出该矩阵的有效秩（即AR阶数的估计值）后，再用课本66页的总体最小二乘方法求得AR参数的估计值a1,a2,K,ap。
（5.）将所求得的AR参数a1,a2,K,ap代入
（1）中，得该谐波过程的特征方
程，求解该特征方程得其模为1的共轭复根。
（6.）由式fi=|arctan[Im（zi）/Re（zi）]|/2p即可求得正弦波频率的估计值。
一一 AR参数和正弦波频率估计值的统计结果
根据以上理论及方法分别编写最小二乘和总体最小二乘的matlab
程序，各独立运行20次，其结果及统计结果如下：
1.最小二乘
（1．）AR阶数为4arMatrix=
Columns1through12
-0.4598
-0.3507
-0.1278
-0.0943
-0.3434
-
0.0737
-0.2146
-0.2919
-0.3176
-0.2331
-0.7039
0.3275
0.8806
0.7543
0.8032
0.6490 0.9086
0.8011
0.6412
0.9384
0.6871
1.0946 1.0373
0.7350
0.1693
0.3133
0.4179
0.5340 0.2238
0.4562
0.4686
0.2366
0.3776
0.1731 -0.0778
0.3412
-0.0202
-0.0758
0.1129
-0.0202 0.0842
0.1405
-0.1023
0.1457
-0.1219
0.3371 -0.0110 -
0.0785
Columns13through20
-0.1600 -0.2792 -0.2207 -0.1843 -0.7587 -
0.2687 -0.2143 -0.4826
0.8284 0.7197 0.7327 0.6541 0.8180
1.0069 0.7344 0.8866
0.3846 0.3781 0.4042 0.4778 0.0283
0.2056 0.4112 0.1485
0.1174 -0.0642 -0.0133 -0.0708 -0.2617
0.2302 -0.0082 -0.0244
ar_mean=
-0.3053
0.8156
0.3036
0.0148
ar_var=
0.0326
0.0174
0.0257
0.0186
fvMatrix=
0 0 0.1200 0 0.1709
0.1207 0 0.1793 0 0.1963 0
0.2000 0.1999 0.1997 0.1998 0.1997
0.1999 0.1997 0.1997 0.1997 0.1998 0.1998
0.1997
0.1333
0.1905
0.1998
0.1995
0.1996 0.1996 0.2000
fv_mean=
0.0556
fv_var=
0.0064
0.0000
其中，arMatrix为4*20的矩阵，用来存储20次独立运行的AR参数估计值，每一列为对应该次运行时所得的AR参数估计值；
ar_mean和ar_var为AR参数估计值的统计结果，ar_mean为运行20次后AR参数的平均值，ar_var为20次AR参数的方差。
fvMatrix为2*20的矩阵，用来存储20次独立运行的正弦波频率估计值，每一列为对应该次运行时所得的正弦波频率估计值；
fv_mean和fv_var为正弦波频率估计值的统计结果，fv_mean为运行20次后正弦波频率的平均值，fv_var为20次正弦波频率的方差。
（2．）AR阶数为6arMatrix=
-0.4147
-0.4369
-0.1121
-0.3637
-0.3285
0.0468 -0.0936
-0.4142
0.0050
-0.3699
-0.0977
0.3209
1.0172
0.9038
0.6183
1.0506
0.9355
0.7342
0.7958
0.9177
1.0123
1.0642
0.8917
1.0613
0.3353
0.2528
0.3925
0.1859
0.2256
0.5002
0.2227
-0.1215
0.9125
0.3570
0.2892
0.1089
0.0778 -0.1941 0.1282 0.1866 0.1771
0.2078 0.0980 -0.1177 0.2724 -0.0708 0.1956
0.3256
0.1812 0.1873 -0.2136 0.0527 -0.0309 -
0.0761 -0.1245 -0.0807 0.4104 0.3344 -0.0451 -
0.0671
0.0662 -0.1642 0.0924 0.0129 0.0618
0.1047 -0.1642 -0.3487 0.2038 -0.1395 -0.0823
0.0513
-0.0709 -0.2426 -0.5091 -0.2368 -0.3134 -
0.4907 -0.4922 -0.3410
0.6266 0.8746 0.8122 0.7509 1.0563
0.9303 0.8535 0.6779
0.5671 0.4125 0.2955 0.2549 0.2981 -
0.1295 0.3570 0.1475
0.1128 0.0871 -0.2398 -0.0526 0.2745
0.0993 -0.5182 -0.2037
-0.0834 0.0801 0.1475 -0.0549 0.0759 -
0.2053 0.3977 -0.0969
0.1406 0.0408 -0.0559 -0.1252 0.1133 -
0.0711 -0.3375 -0.1985
-0.2845
0.8792
0.2932
0.0423
0.0394
-0.0400
0.0266
0.0208
0.0508
0.0451
0.0344
0.0239
0.0522 0 0.1030 0.1256 0.1242
0.1006 0 0 0.0338 0 0
0.1550
0.1789
0.1996
0.1033
0.1575
0.1527
0.1813
0.1850
0.1691
0.2488
0.2285
0.1999
0.2137
0.2345
0.0864
0.0704
0.1165
0.1416
0.1838
0.1984
0.1806
0.1820 0.1996 0.1995
0.1996 0.1997 0.2172 0.1998 0.1994
0.1996 0.2396 0.2024
0.0484
0.2090
0.0031
0.0006
0.0003
其中，arMatrix为6*20的矩阵，用来存储20次独立运行的AR参数估计值，每一列为对应该次运行时所得的AR参数估计值；
fvMatrix为3*20的矩阵，用来存储20次独立运行的正弦波频率估计值，每一列为对应该次运行时所得的正弦波频率估计值；
（3．）实验结果分析
分析以上两组实验结果，当最小二乘方法取AR阶数为4时，所估计出的频率有一个非常接近0.2，但是另一个几乎为0，即用AR参数构成的关于z的特征方程只有一对共轭复根，另一对非常接近于实轴。
也就是说此时虽然AR阶数与实际相吻合，但对频率的估计却是不准确的；
当最小二乘方法取AR除数为6时，所估计出的三个频率中有一个几近
于0，其余两个分别接近于实际频率0.2和0.213，可以用来作为对本题谐波恢复的频率估计。
2.总体最小二乘
反复运行总体最小二乘程序，可以看到绝大多数情况下AR阶数的估计值为4，由此可以推测观测数据对应的谐波信号个数为2。
因此，为了方便求统计结果，下面的20组数据均取自AR阶数估计值为4的那些运行结果。
（1．）总体最小二乘实验结果
arMatrix=
-1.0733
-1.0733 -1.0733 -1.0733 -
1.0733
2.2742
2.2742 2.2742 2.2742
-1.0684
-1.0684 -1.0684 -1.0684 -
1.0684
0.9931
0.9931 0.9931 0.9931
-1.0733 -1.0733 -1.0733 -1.0733 -1.0733 -
1.0733 -1.0733 -1.0733
2.2742 2.2742 2.2742 2.2742 2.2742
-1.0684 -1.0684 -1.0684 -1.0684 -1.0684 -
1.0684 -1.0684 -1.0684
0.9931 0.9931 0.9931 0.9931 0.9931
ar_var=1.0e-030*
0.2076
0.0130
0.2134
0.1999 0.1999 0.1999 0.1999 0.1999
0.1999 0.1999 0.1999
0.2134 0.2134 0.2134 0.2134 0.2134
0.2134 0.2134 0.2134
1.0e-032*
0.3244
（2．）分析实验结果
首先，在适当的阈值下利用SVD方法确定自相关矩阵的有效秩效果很好。
本实验中当取阈值为0.998时，所估计出的有效秩绝大多数情况下均为4，只有个别时候才会出现3或者5的情况。
其次，20次运行的统计结果表明用总体最小二乘方法所估计出的AR参数和谐波频率与实际值非常吻合。
其中，谐波频率的估计值分别为0.1999和0.2134，十分接近于实际频率0.2和0.213，而AR参数估计值和频率估计值的方差也都远小于最小二乘法中对应的方差，相差约三个数量级。
由此可见，谐波恢复时由维数较大的扩展阶相关矩阵估计有效秩并利用总体最小二乘方法进行AR参数和正弦波频率的估计，其精确程度远高于最小二乘方法。
一一 使用SVD确定样本自相关矩阵有效秩的优点和注意事项
1.优点
奇异值分解（SVD）方法是一种典型的线性代数定阶法，它主要用于
求解线性方程组，与该方程组相关联的矩阵不仅表征所期望的解的特征，
而且还常常传达动态性能的信息，在定理的支持下可以简便直观的确定样本自相关的有效秩。
2.注意事项
选择合适的阈值。
若阈值取的太低，则所估计的有效秩小于实际阶数，因此，阈值应取的尽量大一些。
但若阈值取的过高，也会造成有效秩大于实际阶数。
在本实验中，阈值取为0.998。
一一 附页：
程序代码
1.最小二乘方法
functionfangzhen1_ls（）pe=6;
%4
%生成增广矩阵Re
M=100;
loops=20;
%pe=25;
%Re=zeros（M,pe+1）;
fvMatrix=zeros（pe/2,loops）;
arMatrix=zeros（pe,loops）;
forloop=1:
1:
loops
% w=wgn（1,2000,0）;
w=randn（2000,1）;
forn=1:
128
x（n）=（20^（1/2））*sin（2*pi*0.2*n）+（2^（1/2））*sin（2*pi*0.213*n）+w（n+50*loop）;
endRxx=xcorr（x,'
unbiased'
）;
fori=1:
M
forj=1:
pe+1
Re（i,j）=Rxx（pe+i+1-j+128）;
end
kk=pe;
%最小二乘法估计ARMA模型的参数
b=-1*Re（:
[1]）;
A=Re（:
[2:
kk+1]）;
%ar_lsline=（A'
*A）\（A'
*b）;
ar_lsline=inv（A'
*A）*A'
*b;
arMatrix（:
loop）=ar_lsline;
%求谐波频率forj=2:
kk+1
fc（j）=ar_lsline（j-1）;
fc
（1）=1;
fz=roots（fc）;
count=1;
2:
kk
fw（count）=atan（imag（fz（j））/real（（fz（j））））;
fv（count）=fw（count）/（2*pi）;
count=count+1;
endfv=abs（fv）;
fv=sort（fv）;
fvMatrix（:
loop）=fv'
;
%分别求AR参数估计和频率估计的均值和方差
arMatrix
ar_mean=mean（arMatrix'
）'
ar_var=var（arMatrix'
fvMatrix
fv_mean=mean（fvMatrix'
fv_var=var（fvMatrix'
2.总体最小二乘方法
functionfangzhen1_tls（）
%生成增广矩阵ReM=95;
pe=30;
loops=20;
fvMatrix=zeros（2,loops）;
arMatrix=zeros（4,loops）;
kk=0;
forloop=1:
loopswhilekk~=4
Re=zeros（M,pe+1）;
%w=randn（2000,1）;
w=wgn（1,2000,0）;
forn=1:
x（n）=（20^（1/2））*sin（2*pi*0.2*n）+（2^（1/2））*sin（2*pi*0.213*n）+w（n+loop*50）;
end
Rxx=xcorr（x,'
pe+1Re（i,j）=Rxx（pe+i+1-j+128）;
%求Re的有效秩kk[U,S,V]=svd（Re）;
bottom=0;
forhh=1:
bottom=bottom+S（hh,hh）^2;
endoverhead=0;
v=overhead/bottom;
whilev<
0.998
kk=kk+1;
overhead=overhead+S（kk,kk）^2;
v=overhead/bottom;
endkk;
%求Sp
Sp=zeros（kk+1）;
forj=1:
fori=1:
pe+1-kkVj=V（:
j）;
Sp=Sp+S（j,j）^2*（Vj（[i:
i+kk]））*（Vj（[i:
i+kk]））'
endinvSp=inv（Sp）;
fc=invSp（:
1）/invSp（1,1）;
arMatrix（:
loop）=fc（2:
5）;
fz=roots（fc）;
count=1;
fw（count）=atan（imag（fz（j））/real（（fz（j
升级会员 