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\');