MIMO雷达DOA估计-MATLAB仿真源代码

%% ---- * ---- * ---- * ---- * 发射圆阵的双基地 MIMO 雷达角度估计 * ---- * ---- * ---- * ---- % Create Time :2022 03 14 % Finish Time :2022 03 15 % Author :Xu Y. B.（CSDN ID:在路上，正出发） % MATLAB Edition:R2021a % Sketch : % -1- 信号模型建立（波达方向、非相干信源数目固定为 1 ） % -2- 单目标快速参数估计与定位 % Reference ； % [1]陈浩. MIMO雷达参数估计与数据融合方法研究[D].南京大学,2018. % 主要参考文献 3.2 小节 %% Clear clc ; clearvars ; close all ; set(0,'defaultfigurecolor','w') ; %% 信号模型建立 % -1- 阵列信号参数配置 % -1.1- 电磁波参数 fc = 1e9; % 电磁波频率 单位：Hz c = 3e8; % 电磁波传播速度 单位：m/s lambda = c/fc; % 电磁波波长 单位：m % -1.2- 阵元参数 M_t = 32; % 发射阵元数 -圆阵 r = 1; % 发射阵元的分布半径 单位：m 【注】r/lambda >1/4 避免相位模糊 N_r = 16; % 接收阵元数 -线阵 d = lambda/2; % 接收阵元间距 单位：m La = 10e3; % 发射站圆阵中心与接收站参考点的间距 单位：m % -1.3- 目标参数 alpha = 50; % 目标DOA 单位：度 范围：[0,90] theta = 10; % 目标高度角 单位：度 范围：[0,90] phi = -40; % 目标方位角 单位：度 范围：[-180,180] v = 10; % 目标径向速度 单位：米每秒 fd = 2*v/lambda; % 目标的多普勒频移 单位：Hz A = 1; % 目标幅度 单位：V % -1.4- 时域参数（信号波形设为LFM波形） N = 1024; % 快拍（时域样点）总数 fs = 10e3; % 时域采样率 单位：Hz sigma2_n = 1e-3; % 高斯白噪声的功率 σ平方（建议不要设置过大，此参数设置不合理会严重影响角度估计值） % -2- 阵列信号构建过程 % -2.1- 基带信号构建 t = (0:N-1)/fs; % 时间轴 单位：秒 S = A*exp(1j*2*pi*fd*t); % 信号的时域波形矩阵 % -2.2- 矩阵 A 构建 a_alpha = exp(1j*2*pi*(d/lambda)*((0:(N_r-1)).')*cos(alpha/180*pi)); b_theta_phi = exp(1j*2*pi*r*sin(theta/180*pi)*cos(phi/180*pi-2*pi*((0:(M_t-1)).')/M_t)/lambda); A = Kronecker_Mult_(a_alpha,b_theta_phi); % -2.3- 无噪声阵列信号信号构建 X = A*S; % -2.4- 阵列信号加噪声 n = sqrt(sigma2_n/2)*randn(size(X))+1j*sqrt(sigma2_n/2)*randn(size(X)); % 复高斯噪声矩阵 X_array = X+n; %% 单目标快速参数估计与定位 % -1.1- 获取 Es Cx = (X_array*X_array')/N; % 式3.8 [Cx_V,Cx_eig] = eig(Cx); % 特征值分解 Cx_eig = diag(Cx_eig); % 取特征值 转为列向量 [eigen_sort,eigen_index] = sort(abs(Cx_eig),'descend'); % 特征值降序排列 Es = Cx_V(:,eigen_index(1)); % 与最大特征值对应的特征向量 % -1.2- 目标DOA估计 % -1.2.1- 定义选择矩阵 J1 = Kronecker_Mult_([eye(N_r-1),zeros(N_r-1,1)],eye(M_t)); % 式3.9 J2 = Kronecker_Mult_([zeros(N_r-1,1),eye(N_r-1)],eye(M_t)); % 式3.10 % -1.2.2- 矩阵选择 E1 = J1*Es; E2 = J2*Es; % -1.2.3- 求最小二乘解 Phi_r = (pinv(E1))*E2; % -1.2.4- 得到DOA估计值 alpha_Est = (acos(lambda*phase(Phi_r)/(2*pi*d)))/pi*180; % 单位：度 % -1.3- 目标的俯仰、方位角估计 % -1.3.1- 定义选择矩阵 Ar = conj(exp(1j*2*pi*(d/lambda)*((0:(N_r-1)).')*cos(alpha_Est/180*pi))); J3 = Kronecker_Mult_(Ar,eye(M_t)); % -1.3.2- 计算Xb及其协方差矩阵 Xb = (J3.')*X_array-mean(eigen_sort(2:end)); R = (Xb*Xb')/N; % -1.3.3- 求解 F q_p = 4; % 可变 H = [cos(2*pi*((0:M_t-1-q_p).')/M_t)-cos(2*pi*((0:M_t-1-q_p).'+q_p)/M_t),sin(2*pi*((0:M_t-1-q_p).')/M_t)-sin(2*pi*((0:M_t-1-q_p).'+q_p)/M_t)]; U = zeros(1,M_t-q_p); for i = 1:M_t-q_p U(i) = phase(R(i,i+q_p)); end F = (inv(H.'*H))*H.'*U.'; % -1.3.4- 估计θ、φ theta_Est = asin(lambda*abs(F(1)+1j*F(2))/(2*pi*r))/pi*180; phi_Est = phase(F(1)+1j*F(2))/pi*180; %% 【售后】 % 如果对代码有疑问不理解的，可在 CSDN 私信我，有时间都会看私信，会回复的。 % 另外欢迎关注，会不定时上传 #信号处理算法# 的MATLAB源码资源。 % ——在路上，正出发 % 2022-03-15