互质阵列DOA估计matlab

最新推荐文章于 2023-09-19 13:50:47 发布

weixin_46150551

最新推荐文章于 2023-09-19 13:50:47 发布

阅读量1.3k

点赞数 1

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/weixin_46150551/article/details/128860657

版权

该文探讨了使用稀疏嵌套阵列进行DOA（方向-of-arrival）估计的方法，通过MATLAB代码实现。在无噪声情况下，DOA估计效果良好，但当加入噪声后，估计性能下降。代码中展示了阵列结构、信号源生成、噪声添加及音乐谱方法用于DOA估计的过程。

摘要由CSDN通过智能技术生成

阵列结构

在这里插入图片描述

来自 1 J. Li, Y. He, P. Ma, X. Zhang and Q. Wu, “Direction of Arrival Estimation Using Sparse Nested Arrays With Coprime Displacement,” in IEEE Sensors Journal, vol. 21, no. 4, pp. 5282-5291, 15 Feb.15, 2021, doi: 10.1109/JSEN.2020.3034761.

matlab code

clc;
clear;
close all;
YYY=zeros(16,125);
l=200;
m=125;
L=4;%天线数
M=3;%子阵间隔
N=7;%初始间隔
P=3;%信源数
snap=200;%快拍数
C=3e8;
f0=10e6;
lamda=C/f0;
d=0.5*lamda;
snr=10;
theta0=10;
theta1=20;
theta2=40;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:L-1]';%阵列矢量
derad = pi/180;
% %%%相干信号源
% s0=exp(j*2*pi*f0*t);
% s1=exp(j*2*pi*2*f0*t);
% s2=exp(j*2*pi*5*f0*t);


%阵列1流行矢量
a1_theta0=exp(j*2*pi*d/lamda*a*M*sin(theta0/180*pi));
a1_theta1=exp(j*2*pi*d/lamda*a*M*sin(theta1/180*pi));
a1_theta2=exp(j*2*pi*d/lamda*a*M*sin(theta2/180*pi));

A1=[a1_theta0 a1_theta1 a1_theta2];%子阵的导向矢量，每个子阵完全相同

%阵列2流行矢量
a2_theta0=exp(-j*L*2*pi*d/lamda*a*M*sin(theta0/180*pi));
a2_theta1=exp(-j*L*2*pi*d/lamda*a*M*sin(theta1/180*pi));
a2_theta2=exp(-j*L*2*pi*d/lamda*a*M*sin(theta2/180*pi));

A2=[a2_theta0 a2_theta1 a2_theta2];%子阵的导向矢量，每个子阵完全相同
%偏移因子
beta=[exp(j*pi*N*sin(theta0/180*pi)),exp(j*pi*N*sin(theta1/180*pi)),exp(j*pi*N*sin(theta2/180*pi))];

fai=diag(beta);

% S=[s0;s1;s2];
S=randn(M,l*m);%信号源 
%Khatri-Rao空间 构建数据矩阵
for k=1:m
    X1=A1*fai*S(:,1+(k-1)*200:200*k);
    X2=A2*S(:,1+(k-1)*200:200*k);
    X1=awgn(X1,snr,'measured');
    X2=awgn(X2,snr,'measured');%加噪声
    R_c=X1*X2'/200;
    r=reshape(R_c,[],1);
    YYY(:,k)=r;
end



ppp=eye(m)-(1/m)*ones(m);
yyy=YYY*ppp;
[UU,SS,VV]=svd(yyy);

UN=UU(:,7:end);

for iang = 1:361
    angle(iang)=(iang-181)/2;
    phim=derad*angle(iang);
    a1=exp(j*2*pi*d/lamda*a*M*sin(phim)); 
    a2=exp(-j*L*2*pi*d/lamda*a*M*sin(phim));
    beta1=exp(j*pi*N*sin(phim));
    ka = kron(conj(a2),a1*beta1)
    Pmusic(iang)=1/(ka'*UN*UN'*ka);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic);
%Pmusic=Pmusic/Pmmax
Pmusic=10*log10(Pmusic/Pmmax);            % 归一化处理
h=plot(angle,Pmusic);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;