function fangzhen1_tls()
%生成增广矩阵Re
M=95;
pe=30;
loops=20;
fvMatrix=zeros(2,loops);
arMatrix=zeros(4,loops);
kk=0;
for loop=1:1:loops
while kk~=4
Re=zeros(M,pe+1);
%w=randn(2000,1);
w=wgn(1,2000,0);
for n=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);
%x(n)=(20^(1/2))*sin(2*pi*0.2*n);
end
Rxx=xcorr(x,'unbiased');
for i=1:1:M
for j=1:1:pe+1
Re(i,j)=Rxx(pe+i+1-j+128);
end
end
%求Re的有效秩kk
[U,S,V]=svd(Re);
bottom = 0;
for hh=1:1:pe+1
bottom = bottom + S(hh,hh)^2;
end
overhead=0;
v = overhead/bottom;
kk=0;
while v<0.998
kk = kk + 1;
overhead = overhead + S(kk,kk)^2;
v = overhead/bottom;
end
end
kk;
%求Sp
Sp=zeros(kk+1);
for j=1:1:kk
for i=1:1:pe+1-kk
Vj=V(:,j);
Sp=Sp+S(j,j)^2*(Vj([i:1:i+kk]))*(Vj([i:1:i+kk]))';
end
end
invSp=inv(Sp);
fc=invSp(:,1)/invSp(1,1);
arMatrix(:,loop)=fc(2:5);
fz=roots(fc);
count=1;
for j=1:2:kk
fw(count)=atan(imag(fz(j))/real((fz(j))));
fv(count)=abs(fw(count)/(2*pi));
count = count+1;
end
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')'
评论0