function fangzhen1_ls()
pe=6;
%生成增广矩阵Re
M=100;
loops=20;
%pe=25;
%Re=zeros(M,pe+1);
fvMatrix=zeros(pe/2,loops);
arMatrix=zeros(pe,loops);
for loop=1:1:loops
% w = wgn(1,2000,0);
w = randn(2000,1);
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+50*loop);
%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
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;
%求谐波频率
for j=2:1:kk+1
fc(j) = ar_lsline(j-1);
end
fc(1) = 1;
fz=roots(fc);
count=1;
for j=1:2:kk
fw(count)=atan(imag(fz(j))/real((fz(j))));
fv(count)=fw(count)/(2*pi);
count=count+1;
end
fv=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')'
