p319.m

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

clear all;

N = 128;
m =80;
Pi = 3.1415926;
n = [0:1:N-1];
fsin=0;

Times=100; %独立仿真次数
for k=1:Times 
    
x = zeros(m,N); 
for i=1:m     
    rand('state',sum(100*clock));
    w = randn(1,N);
    x(i,:) = sqrt(20)*sin(2*Pi*0.2*(n+i-1))+sqrt(2)*sin(2*Pi*0.213*(n+i-1))+w;
end
Rx = x*x'/N;

%[U,d,V]=svd(R);U为左奇异矩阵，d为对角阵
[V,d,U]=svd(Rx);
for i=1:m
    PP(i)=d(i,i);
end
for i=1:m
    if (PP(i)/PP(1))<0.05
        p=i-1;
        break;
    end
end
S = U(:,1:p);
G = U(:,p+1:m);
sigma = sum(PP(p+1:m))/(m-p);
U1=zeros(m,m);
for i=1:p
    U1 = U1+S(:,i)*(S(:,i)')*PP(i)/(sigma-PP(i))^2;
end
U1 = U1*sigma;

dw=2*Pi*0.001;
df=dw/(2*Pi);
L = round(Pi/dw);
%function a = geta(w,m)
for i=1:L
    a = geta((i-1)*dw,m);
    P(i)=(a'*U1*a)/(a'*G*G'*a);
end
P=abs(P);
num=1;

Pmax = max(P);
for i=2:L-1
    if  P(i)>0.1*Pmax & P(i)>P(i-1) & P(i)>P(i+1)
        fsin(k,num) = (i-1)*df;
        num=num+1;
    end
end
end

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

uf
df
realf=[0.2 0.213]
offset=abs(realf-uf)


W = [0:L-1]*dw;
f=W/(2*Pi);
%plot(f,P);
%axis([0.15,0.25,0,10]);