现代信号处理仿真作业一(3.18谐波恢复).docx
- 文档编号:235360
- 上传时间:2023-04-28
- 格式:DOCX
- 页数:16
- 大小:25.42KB
现代信号处理仿真作业一(3.18谐波恢复).docx
《现代信号处理仿真作业一(3.18谐波恢复).docx》由会员分享,可在线阅读,更多相关《现代信号处理仿真作业一(3.18谐波恢复).docx(16页珍藏版)》请在冰点文库上搜索。
现代信号处理仿真作业一(3.18谐波恢复)
班级:
自研42 姓名:
李琳琳 学号:
2004211068
一一谐波恢复的基本理论与方法:
1.Pisarenko谐波分解理论
谐波过程可用差分方程描述,首先利用Pisarenko谐波分解理论推导谐波过程所对应的差分方程。
对单个正弦波x(n)=sin(2pfn+q)利用三角函数恒等式,有:
x(n)-2cos(2pf)x(n-1)+x(n-2)=0
对上式作z变换,得:
[1-2cos(2pf)z-1+z-2]X(z)=0
得到特征多项式:
1-2cos(2pf)z-1+z-2=0。
由此可见,正弦波的频率可以
由相应特征方程的一对共轭根来决定:
fi=|arctan[Im(zi)/Re(zi)]|/2p
将单个正弦波推广到多个正弦波的情形,得:
如果p个实的正弦波信号没有重复频率的话,则这p个频率应该由特征多项式
1 2p-1
1+az-1+K+a z-(2p-1)+z-2p=0
(1)
的根决定。
由此可得到p个实正弦波所组成的谐波过程可以用以下的差分方程进行描述:
å
2p
x(n)+ aix(n-i)=0
i=1
这是一个无激励的AR过程。
2.谐波恢复的ARMA建模法
å
2p
在无激励的AR模型差分方程x(n)+ aix(n-i)=0两边同乘x(n-k),
i=1
并取数学期望,则有:
å
2p
Rx(k)+ aiRx(k-i)=0,"k
(2)
i=1
正弦波过程一般是在加性白噪声中被观测的,设加性白噪声为w(n),
即观测过程为:
p
y(n)=x(n)+w(n)=åAisin(2pfin+qi)+w(n),其中
i=1
w(n)为0均值的高斯白噪声。
由于谐波过程与加性白噪声统计独立,所以有:
R(k)=R(k)+R(k)=R(k)+s2d(k) (3)
y x w x w
将(3)式代入
(2)式得:
å
2p
Ry(k)+ aiRy(k-i)=0,k>2p (4)
i=1
这就是加性白噪声中观测到的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=
Columns1through12
0 0 0.1200 0 0.1709
0.1207 0 0.1793 0 0.1963 0
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
Columns13through20
0.1333
0
0
0
0
0.1905
0
0
0.1997
0.1998
0.1997
0.1997
0.1995
0.1996 0.1996 0.2000
fv_mean=
0.0556
0.1997
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=
Columns1through12
-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
Columns13through20
-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
ar_mean=
-0.2845
0.8792
0.2932
0.0423
0.0394
-0.0400
ar_var=
0.0266
0.0208
0.0508
0.0451
0.0344
0.0239
fvMatrix=
Columns1through12
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.1997
0.1575
0.1527
0.1813
0.1998
0.1997
0.1997
0.1850
0.1691
0.1996
0.2488
0.1995
0.2285
0.1996
0.1996
0.1999
0.2137
0.1998
0.2345
0.1997
0.1995
Columns13through20
0.0864
0.0704
0
0
0.1165
0
0
0
0.1416
0.1838
0.1996
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
fv_mean=
0.0484
0.1806
0.2090
fv_var=
0.0031
0.0006
0.0003
其中,arMatrix为6*20的矩阵,用来存储20次独立运行的AR参数估计值,每一列为对应该次运行时所得的AR参数估计值;ar_mean和ar_var为AR参数估计值的统计结果,ar_mean为运行20次后AR参数的平均值,ar_var为20次AR参数的方差。
fvMatrix为3*20的矩阵,用来存储20次独立运行的正弦波频率估计值,每一列为对应该次运行时所得的正弦波频率估计值;fv_mean和fv_var为正弦波频率估计值的统计结果,fv_mean为运行20次后正弦波频率的平均值,fv_var为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=
Columns1through12
-1.0733
-1.0733
-1.0733 -1.0733 -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
2.2742
2.2742
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
-1.0684 -1.0684 -1.0684 -
1.0684
0.9931
0.9931
0.9931 0.9931 0.9931
0.9931
0.9931
0.9931
0.9931 0.9931 0.9931
0.9931
Columns13through20
-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
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
0.9931 0.9931 0.9931
ar_mean=
-1.0733
2.2742
-1.0684
0.9931
ar_var=1.0e-030*
0.2076
0.2076
0.2076
0.0130
fvMatrix=
Columns1through12
0.1999
0.1999
0.1999
0.1999
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
0.2134
0.2134
0.2134
0.2134
Columns13through20
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
fv_mean=
0.1999
0.2134
fv_var=
1.0e-032*
0.3244
0
其中,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.)分析实验结果
首先,在适当的阈值下利用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:
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:
1:
M
forj=1:
1:
pe+1
Re(i,j)=Rxx(pe+i+1-j+128);
end
end
kk=pe;
%最小二乘法估计ARMA模型的参数
b=-1*Re(:
[1]);
A=Re(:
[2:
1:
kk+1]);
%ar_lsline=(A'*A)\(A'*b);ar_lsline=inv(A'*A)*A'*b;arMatrix(:
loop)=ar_lsline;
%求谐波频率forj=2:
1:
kk+1
fc(j)=ar_lsline(j-1);
end
fc
(1)=1;
fz=roots(fc);count=1;
forj=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';
end
%分别求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:
1:
loopswhilekk~=4
Re=zeros(M,pe+1);
%w=randn(2000,1);w=wgn(1,2000,0);
forn=1:
1:
128
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,'unbiased');fori=1:
1:
M
forj=1:
1:
pe+1Re(i,j)=Rxx(pe+i+1-j+128);end
end
%求Re的有效秩kk[U,S,V]=svd(Re);
bottom=0;
forhh=1:
1:
pe+1
bottom=bottom+S(hh,hh)^2;
endoverhead=0;
v=overhead/bottom;kk=0;
whilev<0.998
kk=kk+1;
overhead=overhead+S(kk,kk)^2;v=overhead/bottom;
end
endkk;
%求Sp
Sp=zeros(kk+1);forj=1:
1:
kk
fori=1:
1:
pe+1-kkVj=V(:
j);
Sp=Sp+S(j,j)^2*(Vj([i:
1:
i+kk]))*(Vj([i:
1:
i+kk]))';
end
endinvSp=inv(Sp);
fc=invSp(:
1)/invSp(1,1);
arMatrix(:
loop)=fc(2:
5);fz=roots(fc);
count=1;
forj=1:
2:
kk
fw(count)=atan(imag(fz(j))/real((fz(j
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 信号 处理 仿真 作业 3.18 谐波 恢复
![提示](https://static.bingdoc.com/images/bang_tan.gif)