p3181fai.m

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

%一般最小二乘方法（引入随机相位）
clear all;

N = 128;
Pi = 3.1415926;
n = [1:1:128];
p=6;   %可调AR阶数
q=16;


Times=100; %独立仿真次数
f=zeros(Times,2);
for k=1:Times
    
w = randn(1,N);    
fai1=rand*2*Pi-Pi;
fai2=rand*2*Pi-Pi;
x = sqrt(20)*sin(2*Pi*0.2*n+fai1)+sqrt(2)*sin(2*Pi*0.213*n+fai2)+w;
%r=R(N,x,T)       %N is the data number,T is the time delay

Re=zeros(p,p);
for i=1:p
    for j=1:p
        T=q+1-p+(i-1)+(j-1);
        Re(i,j)=R(N,x,T);
    end
end

for i=1:p
    r(i,1)=-R(N,x,q+i);
end
A=Re;
a1=((A.'*A)^(-1))*(A.')*r;
for i=1:p
    a(k,i)=a1(p+1-i);
end
z1=roots([1,a(k,:)]);
z=0;
num=1;
for i=1:p
    if 0.95<abs(z1(i))&abs(z1(i))<1.05
        z(num)=z1(i);
        num=num+1;
    end
end


num=1;
for i=1:length(z)  
    temp=atan(imag(z(i))/real(z(i)))/(2*Pi);
    if temp>0
        f(k,num)=temp;
        num=num+1;
    end
end
end

%让f从大到小排序
f=sort(f,2,'descend');

ua=sum(a)/Times;
da=0;
for i=1:Times
    da=da+((a(i,:)-ua).^2);
end
da=da/Times;

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

ua
da
uf

realf=[0.2 0.213]
temp2=size(f);
realff=zeros(1,temp2(2));
realff(1:2)=[0.2 0.213];
df
offset=abs(realff-uf)