http://bk.treeblog.top/index.php/category/xyjl/ 关于学习记录的回忆 http://bk.treeblog.top/index.php/archives/173/ 2025-06-20T23:09:44+08:00 主要学习来源于这篇文章：《Joint DOD and DOA estimation of bistatic MIMO radar in the presence of unknown mutual coupling》大致思路是：先回忆了下 ULA 的 ESPRIT 方法；然后是 MIMO 下的 ESPRIT 方法然后是看文章看不懂，然后 AI 帮我解析，先读出大致的框架，然后丰富细节，最后对自己的理解进行确认，然后画出代码仿真的关键点下面在这个网盘下面提供 2 个 PDF 和一个代码方便之后回顾学习https://a.siyouyun.ren:30597/OutSaveFile2/RadarXG/ESPRIT_Double_mimomatlab 代码预览clc,clear all,close all; %% Parameters of the radar system M = 7; % number of Tx N = 7; % number of Rx d = 0.5; % inter-element space % theta = [15 -10 5]; % DOA % phi = [20 35 60]; theta = [-50 20 60]; % DOD % DOA phi = [20 40 60]; K = length(theta); % number of target % dopplor frequency shift L = 100; % number of sanpshot SNR = 30; % signal-to-noise ratio before matched filtering Geo_Rx = [0:N-1]; % geometry of Rx Geo_Tx = [0:M-1]; % geometry of Tx At = exp(-j*2*pi*d*Geo_Tx.'*sind(phi)); % transmitting direction matrix Ar = exp(-j*2*pi*d*Geo_Rx.'*sind(theta)); % receiving direction matrix item = 10; %% 设置互耦系数; mc_K = 2; mc_nb = mc_K+1; C_t = generateMCMmat(M,mc_nb); C_r = generateMCMmat(N,mc_nb); C = kron(C_t,C_r); %% 设置选择Pt和Pt矩阵 Pt = [zeros(M-2*mc_K,mc_K),eye(M-2*mc_K),zeros(M-2*mc_K,mc_K)]; Pr = [zeros(N-2*mc_K,mc_K),eye(N-2*mc_K),zeros(N-2*mc_K,mc_K)]; M1 = M-2*mc_K; N1 = N-2*mc_K; P = kron(Pr,Pt); % number of trials %% Selective matrices JM1 = [eye(M1-1),zeros(M1-1,1)]; JM2 = [zeros(M1-1,1),eye(M1-1)]; JN1 = [eye(N1-1),zeros(N1-1,1)]; JN2 = [zeros(N1-1,1),eye(N1-1)]; Jt1 = kron(JM1,eye(N1)); Jt2 = kron(JM2,eye(N1)); Jr1 = kron(eye(M1),JN1); Jr2 = kron(eye(M1),JN2); a_kronrao = zeros(M*N,K); for i=1:K a_kronrao(:,i) = kron(At(:,i),Ar(:,i)); end for item_num = 1:item disp(['SNR = ',num2str(SNR),' dB, ',num2str(item_num), ' # try : ']); %% Matched Filtering S = randn(K,L)+1i*randn(K,L); X_0 = C*a_kronrao * S; X = awgn(X_0,SNR,"measured","dB"); Y = P * X; %% Eigen decomposition %Rx = (X*X')/L; % Estimated covariance matrix Rx = (Y*Y')/L; [Es,D] = eigs(Rx,K,'LM'); % Signal subspace %% Rotational invariant property Vt = pinv(Jt2*Es)*Jt1*Es; %Vt = pinv(Jt1*Es)*Jt2*Es; [T,Phit] = eig(Vt); Phir = inv(T)*pinv(Jr2*Es)*Jr1*Es*T; %Phir = inv(T)*pinv(Jr1*Es)*Jr2*Es*T; E_theta = asind(angle(diag(Phir))/pi).'; E_phi = asind(angle(diag(Phit))/pi).'; plot(E_theta,E_phi,'k*');hold on; end H(1)=plot(E_theta,E_phi,'k*');hold on; H(2)=plot(theta,phi,'rx','markersize',28);grid on; xlabel('DOA'),ylabel('DOD'); legend([H(1),H(2)],'Esimated','Ture')闲聊没想到 matlab 还可以画出这样的图片，“赛博生命”，hhhhttps://mp.weixin.qq.com/s/GS-jSHnasVxrrIn5C2E3WQfigure('Position',[300,50,900,900], 'Color','k'); axes(gcf, 'NextPlot','add', 'Position',[0,0,1,1], 'Color','k'); axis([0, 400, 0, 400]) SHdl = scatter([], [], 2, 'filled','o','w', 'MarkerEdgeColor','none', 'MarkerFaceAlpha',.4); t = 0; i = 0:2e4; x = mod(i, 100); y = floor(i./100); k = x./4 - 12.5; e = y./9 + 5; o = vecnorm([k; e])./9; while true     t = t + pi/90;     q = x + 99 + tan(1./k) + o.*k.*(cos(e.*9)./4 + cos(y./2)).*sin(o.*4 - t);     c = o.*e./30 - t./8;     SHdl.XData = (q.*0.7.*sin(c)) + 9.*cos(y./19 + t) + 200;      SHdl.YData = 200 + (q./2.*cos(c));     drawnow end