chapter7_4.asv

阵列信号处理中各中方位来波方向估计

%%%%chapter7_4%%%%%
clc;clear;
P=2;%信号源个数
fai=[30,40]*pi/180;%信号的方向转为弧度
fai=pi*sin(fai);
theta=linspace(-100,100,1000);
L=8;
N=100;%N为快拍数
%%%产生A_theta%%%%
A_theta=zeros(L,2);
for k=1:L
    A_theta(k,1)=exp(j*fai(1)*(k-1));
    A_theta(k,2)=exp(j*fai(2)*(k-1));
end
%%%% 产生信号%%%%
Sig=zeros(2,N);
for k=1:N
    Sig(1,k)=exp(j*pi/4*(k-1));
    Sig(2,k)=exp(j*pi/6*(k-1));
end
SNR=[-5,0,5];%分别是三个信噪比
R_guji=zeros(L,L);
y_gce=zeros(L,1);%实际的观测值
sums=zeros(L,L);
for k=1:N
    y_gce=1/sqrt(2)*sqrt(10.^(SNR(3)/10))*A_theta*Sig(:,k)+randn(L,1);%SNR(i)分别代表三个信噪比的取值
    sums=sums+y_gce*y_gce';
end
R_guji=sums/N;
%%%对估计的相关函数进行svd分解求出特征值，求得的u和v是一样的特征向量
[u,s,v]=svd(R_guji);
%%%算术平均变化率判别法%%%%
% lamda=zeros(1,L);
% for k=1:L
%     lamda(k)=s(k,k)
% end

deta_lamda1=(s(1,1)-s(L,L))/(L-1);
ik=zeros(1,L);
k=1;
for n=1:L-1
    data_lamda_1=s(n,n)-s(n+1,n+1);
    if data_lamda_1<=deta_lamda1
        ik(k)=n;%求出所对应的序列集
        k=k+1;
    end
end
for i=1:length(ik)-1
    if ik(i+1)==ik(i)+1
       p1=ik(i)-1%p为要求的信号源数目
         break;
    end 
end
%%%%几何平均变化率判别法%%%%
deta_lamda2=log(s(1,1)/s(n,n))/(L-1);
ik2=zeros(1,L);
k=1;
for n=1:L-1
    data_lamda_2=log(s(n,n)/s(n+1,n+1));
    if data_lamda_2<=deta_lamda2
        ik2(k)=n;%求出所对应的序列集
        k=k+1;
    end
end
for i=1:length(ik2)-1
    if ik2(i+1)==ik2(i)+1
       p2=ik2(i)-1%p2为要求的信号源数目
         break;
    end 
end
%%%特征值概率分布准则判别法%%%
gama_gate=0.01;
 lamda=zeros(1,L);
for k=1:L
    lamda(k)=s(k,k);
end
 e=0.02;
 p3=0;
for p=0:L
     x=[p3+1:L];
     y=lamda(p3+1:L);
     y1=polyfit(x,y,2);
     y_nihe=polyval(y1,x);
     e=sqrt(sum((y-y_nihe).^2)/N);
     if e>gama_gate
         p3=p3+1;
     else break;
     end
end
p3