DOA估计的三种方法

clear all;clc;
%**************输入参数**************************
M=8;%阵元数目
% D=0.5;%阵元间距
D=0.5;%阵元间距
d=2;%入射源数目
theta1=-5/180*pi;%入射角度
theta2=5/180*pi;
f=300*10^6;%中心频率
N=100;%快拍次数
SNR=10;%信噪比
lamda=1;%波长
%*************************入射信号产生*********
%产生s1
T1=randint(1,N*2);
T_IQ1=reshape(T1,2,N)';
symbol1=bi2de(T_IQ1,'left-msb');
Table=(1/sqrt(2))*[-1-j -1+j 1-j 1+j];
s1=Table(symbol1+1);
%产生s2
T2=randint(1,N*2);
T_IQ2=reshape(T2,2,N)';
symbol2=bi2de(T_IQ2,'left-msb');
s2=Table(symbol2+1);
%产生复高斯白噪声
NN=zeros(M,N);
for k=1:M
NN(k,:)=1/(sqrt(2*SNR))*[randn(1,N)+j*randn(1,N)];
% NN(k,:)=1/(sqrt(2*SNR))*[randn(1,N)];
end
%*****************产生接收信号******************
w=2*pi*f;
ta1=D*sin(theta1)/(lamda*f);
ta2=D*sin(theta2)/(lamda*f);
A=zeros(M,2);
% A=zeros(M,1);
X=zeros(M,N); %接收到的信号
for i=1:M
A(i,:)=[exp(j*(i-1)*w*ta1) exp(j*(i-1)*w*ta2)];
% A(i,:)=[exp(j*(i-1)*w*ta1) ];
end
S=[s1;s2];
% S=[s1];
for i=1:M
for k=1:N
X(i,k)=A(i,:)*S(:,k)+NN(i,k);
end
end
%计算自相关矩阵
sum=zeros(M,M);
for i=1:N
sum=sum+X(:,i)*X(:,i)';
end
R_xx=sum/N;
%*****************周期图法***********************
%*****************周期图法***********************
theta=[-90:90]';
tmp=j*w*D*sin(theta'*pi/180)/(lamda*f);
tmp2=[0:M-1]';
e=tmp2*tmp;
E=exp(e);

for i=1:181
P_BT(i)=E(:,i)'*R_xx*E(:,i)/M;
end

S_theta=P_BT;
S=10*log10(S_theta/max(S_theta));%换算成db的输出形式 % 1*181
%画出谱峰图形
figure(2)
plot(-90:90,S,'r*');
grid on;
title('DOA估计图周期图法');
xlabel('方位角（度）');
ylabel('谱峰（dB）');
hold on
%谱峰与信号强度无关,只是反映a_theta与噪声子空间的正交性
%*****************capon法***********************
theta=[-90:90]';
tmp=j*w*D*sin(theta'*pi/180)/(lamda*f);
tmp2=[0:M-1]';
e=tmp2*tmp;
E=exp(e);
invR_xx = inv(R_xx)
for i=1:181
P_BT(i)=E(:,i)'*invR_xx*E(:,i);
end
doa = 1./P_BT;
S_theta=doa;
S=10*log10(S_theta/max(S_theta));%换算成db的输出形式 % 1*181
%画出谱峰图形
% figure(2)
plot(-90:90,S,'g-');
grid on;
title('DOA估计图CAPON方法');
xlabel('方位角（度）');
ylabel('谱峰（dB）');
hold on
%谱峰与信号强度无关,只是反映a_theta与噪声子空间的正交性
%*****************MUSIC法***********************
theta=[-90:90]';
tmp=j*w*D*sin(theta'*pi/180)/(lamda*f);
tmp2=[0:M-1]';
e=tmp2*tmp;
E=exp(e);
p=2;
[V,D] = eig(R_xx);
%分解信号子空间和噪声子空间
Sp=V(:,M-p+1:M)
Nn=V(:,1:M-p) %噪声子空间8*6
for i=1:181
P_BT(i)=(E(:,i)'*Nn)*(E(:,i)'*Nn)';
end
doa = 1./P_BT;
S_theta=doa;
S=10*log10(S_theta/max(S_theta));%换算成db的输出形式 % 1*181
%画出谱峰图形
% figure(2)
plot(-90:90,S,'b-.');
grid on;
title('DOA估计图MUSIC方法');
xlabel('方位角（度）');
ylabel('谱峰（dB）');
%谱峰与信号强度无关,只是反映a_theta与噪声子空间的正交性
% semilogy(theta,doa,'r');