p3201.m

利用一般最小二乘法和SVD-TLS方法 进行AR参数估计和正弦波频率估计

%基本ESPRIT算法1(未加入随机相位)
clear all;

N = 128;
Pi = 3.1415926;
n = [1:1:128];
m= 10;

Times=100; %独立仿真次数
%f=zeros(Times,2);
for k=1:Times
    
w = randn(1,N);    
x = sqrt(20)*sin(2*Pi*0.2*n)+sqrt(2)*sin(2*Pi*0.213*n)+w;
%r=R(N,x,T)       %N is the data number,T is the time delay
Rxx = zeros(m,m);
for i=1:m
    for j=1:m
        if (i-j)>=0
            Rxx(i,j)=R(N,x,(i-j));
        else
            Rxx(i,j)=(R(N,x,(j-i))').';
        end
    end
end
Rxy = zeros(m,m);
for i=1:m
    for j=1:m
        if (i-j)>0
            Rxy(i,j)=R(N,x,(i-j-1));
        else
            Rxy(i,j)=(R(N,x,(j-i+1))').';
        end
    end
end
[V,D] = eig(Rxx);
for i=1:m
    d(i)=D(i,i);
end
%让d从大到小排序
d=sort(d,'descend');
for i=1:m
    if d(i)/d(1)<0.05
        p=i-1;
        break;
    end
end
%sigma=min(d);
sigma=sum(d(p+1:m))/(m-p);

I=diag(ones(1,m));
Z=zeros(m,m);
for i=1:m
    for j=1:m
        if (i-j)==1
            Z(i,j)=1;
        end
    end
end
Cxx=Rxx-sigma*I;
Cxy=Rxy-sigma*Z;
%[V,D] = EIG(A,B) 
[VV,DD] = eig(Cxx,Cxy);
for i=1:m
    dd(i)=DD(i,i);
end

ddd = abs(abs(dd)-1);
[ddd,II] = sort(ddd);

W=0;
num=1;
for i=1:p
    ii=II(i);
    if imag(dd(ii))>=0       
        temp = atan(imag(dd(ii))/real(dd(ii)));
        if real(dd(ii))>=0
            W(num)=temp;
            num=num+1;
        else
            W(num)=temp+Pi;
            num=num+1;
        end
    end
end
l= length(W);
f(k,1:l)=W/(2*Pi);
end

uf = sum(f)/Times;
df=0;
for i=1:Times
    df=df+((f(i,:)-uf).^2);
end
df=df/Times;

uf
df
f;