【雷达仿真 | FMCW TDMA-MIMO毫米波雷达信号处理仿真（可修改为DDMA-MIMO）】

在这里插入图片描述

本文编辑：调皮哥的小助理

本文引用了CSDN雷达博主@XXXiaojie的文章源码（https://blog.csdn.net/Xiao_Jie1），加以修改和注释，全面地、详细地阐述了FMCW TDM-MIMO毫米波雷达的工作原理，同时配套MATLA仿真实现方法，非常适合于雷达刚入门的同学参考学习，并引导大家基于TDMA-MIMO扩展到DDMA-MIMO，进而在宏观上认识雷达，从微观上掌握雷达，形成雷达学习过程中战略与战术的统一。

本文尤其感谢CSDN雷达博主@XXXiaojie在雷达领域贡献的技术文章，帮助了很多人。XXXiaojie是桂林电子科技大学的硕士研究生，本专业不是雷达，而是从计算机专业转到雷达这边的，其代码的阅读和分析能力非常强，对于TI毫米波雷达架构和原理的理解也是十分到位的，感兴趣的朋友可以搜索他的CSDN文章，XXXiaojie博主也在雷达技术交流【1】群中，感兴趣可以相互交流。

立足于本文，雷达初学者下可掌握雷达基本原理，上可继续在本文的基础之上扩展DDMA或其他调制波形（只需修改发射信号模型即可），以及聚类、跟踪等雷达数据处理，甚至是雷达通信一体化（通感一体化）等算法的仿真。在本文分享的MATLAB仿真代码程序中会引出如何生成多帧的雷达数据，用于后续的雷达数据处理。

然而，雷达信号处理和数据处理的理论远远不止限于此，但雷达的科普一直持续在入门阶段，我们不能只停留在这里，因此本文算是雷达入门科普的一个闭环文章，以后希望随着我个人见闻的增加，能够分享更加有深度的雷达科普内容。

〇、基础知识介绍

1、FMCW叫做调频连续波，调频连续波有步进调频连续波（SFMCW）、线性调频连续波（LFMCW），以及其他非线性调频连续波。LFMCW有锯齿波和三角波等等；锯齿波有恒定调制周期的锯齿波和非恒定周期的锯齿波。关于FMCW雷达目标检测的基本原理可以参考Tide官方文档：Introduction to mmwave sensing:FMCW radars.pdf和The fundamentals of millimeter wave sensors。

2、MIMO是多输入多输出（Multiple input Multiple output），也就是多发多收。关于MIMO的基本原理可以参考TI的官方文档：MIMO radar.pdf

TDM-MIMO也就是时分复用的多发多收，也称为TDMA-MIMO。DDMA-MIMO是多普勒多通道分离的多发多收。除了TDMA、DDMA外，还有CDMA等MIMO发射波形，总结如下表：

在这里插入图片描述
（图来自参考论文[1]）

TDM-MIMO也就是多个发射天线分时间轮流发射信号，但每个发射信号之间信号模型都是一样的；但DDMA是多个天线同时发射信号，每个发射信号之间的初始相位是按照chirp数和发射天线数量变化的。具体关于DDMA的原理可以参考TI官方文档：基于AWR2944的汽车雷达DDMA波形的原理和实现.pdf。

上述文章是本文的基本理论基础，如果还要更加深层次了解雷达的基本原理，那么推荐看下面这两篇文章：

1.调皮连续波：雷达入门系列文章（3）| 毫米波雷达信号处理入门教程

2.调皮连续波：雷达原理 | 用MATLAB信号处理是如何解算目标的距离和速度信息的？

好了，理论研究得差不多了，那么我们就开始讲仿真吧。

一、总结概述

本文的仿真流程主要如下所示：

1、雷达参数生成；
2、信号模型构建；
3、IQ通道校正；
4、雷达信号处理之距离估计（1D-FFT）；
5、雷达信号处理之速度估计（2D-FFT）；
6、雷达信号处理之角度估计（3D-FFT）；
7、多通道间的非相干积累；
8、CFAR检测；
9、峰值检索；
10、多普勒补偿
11、DOA估计。

上述过程可以采用多帧重复运行，连续多帧多目标运动的演示效果如下所示：

在这里插入图片描述

二、雷达参数设置

雷达参数设置包括雷达波形参数设计、目标参数设计、射频前端参数选择等等，这些参数都是可以任意修改的，也就是可以随着雷达系统设计人员的需求而尽情发挥，并不只是限于本文。比如77G可以修改为24G或者35G，原理都是相通的，只是频段不同。

其中，需要注意采样点数与采样时间、采样率之间的匹配关系，修改的时候要匹配好，不然会错先错误。同样，雷达探测的最大距离、最大不模糊速度也要设置好，否则，如果没有特别的算法来解决，那就会和预想的结果不一样。

parameter.c = 3e8;             %光速
parameter.stratFreq = 76.5e9;  %起始频率


parameter.Tr = 10e-6;          %扫频时间 也就是周期
parameter.Samples = 256;       %采样点
parameter.Fs = 25.6e6;         %采样率

parameter.rangeBin = parameter.Samples ;      %rangebin
parameter.Chirps = 512;        %chirp数
parameter.dopplerBin = parameter.Chirps;    %dopplerbin

parameter.Slope = 30e12;       %chirp斜率
parameter.Bandwidth = parameter.Slope * parameter.Tr ; %发射信号有效带宽
parameter.BandwidthValid = parameter.Samples/parameter.Fs*parameter.Slope; %发射信号带宽
parameter.centerFreq = parameter.stratFreq + parameter.Bandwidth / 2; %中心频率
parameter.lambda = parameter.c / parameter.centerFreq; %波长

parameter.txAntenna = ones(1,3); %发射天线个数
parameter.rxAntenna = ones(1,4); %接收天线个数
parameter.txNum = length(parameter.txAntenna);
parameter.rxNum = length(parameter.rxAntenna);
parameter.virtualAntenna = length(parameter.txAntenna) * length(parameter.rxAntenna);

parameter.dz = parameter.lambda / 2; %接收天线俯仰间距
parameter.dx = parameter.lambda / 2; %接收天线水平间距
parameter.doaMethod = 2; %测角方法选择 1-dbf  2-fft 3-capon
parameter.target = [
                    100 -20 0; %target1 range speed angle
                    0  10  -30;  %target2 range speed angle
                    0  20   30;  %target2 range speed angle
                    ];

设置雷达参数这部分内容就是属于雷达系统设计了，客户或者我们自己需要设计一个什么样的雷达，满足什么样的指标需求，所涉及到的参数都在这里设置好，因此这些参数完全可以按照我们的想法设计。比如雷达的收发天线可以设计为单发单收、单发多收、多发多收（如12T16R）。

上述代码中，我们设置了三个目标，如同在上述gif动图中见到的那样，目标的速度、距离需要符合实际。如果有同学需要生成许多点云，那就可以把目标的个数搞多一点，还可以加上测量俯仰角的天线，模拟人体三维点云。

TDM-MIMO信号模型

根据之前设置好的雷达系统参数，在这里就可以生成发射信号模型、接收信号模型、混频，以及中频信号模型。

发射信号模型（复信号形式）：

$s_t(t)=\exp \left\{2 \pi\left(f_0 t+u t^2 / 2\right)+\phi_0\right\}, \quad t \in\left[0, T_r\right]$

在这个程序中的信号是按照天线和脉冲依次生成的，所以下面的代码和上述的公式的写法有些不同，这个要明白。

St = exp((1i2pi)(centerFreq(t+( chirpId-1)Tr)+slope/2t.^2)); %发射信号

接收信号模型：
$S_r(t)=K \exp \left\{2 \pi\left[f_0(t-\tau)+u(t-\tau)^2 / 2\right]+\phi_0\right\}, t \in\left[0, T_r\right]$

Sr = 10exp((1i2pi)((centerFreq-fd)(t-tau+( chirpId-1) * Tr)+slope/2(t-tau).^2 + wx)); %回波信号

在程序中，增加了不同天线的相位差项。以后不管是运动目标还是静止目标，信号模型都统一写为一个，其中静止目标只是速度为零罢了。如下列代码中targetSpeed =0。而当遇到多目标时，直接对接收信号求和即可。

tau=2*(targetRange+targetSpeed*(txId-1)*Tr)/c;

然后可以增加点噪声：

Sif=awgn(Sif,20);%叠加20dB高斯白噪声

生成信号模型的代码如下，这里的信号模型也就是指后续处理需要的数据，其中在目标参数方面我修改了一些，主要是让多帧连续运动的目标看起来更贴近实际运动。比如说，不能让两个速度不同的目标，看起来跑的一样快，再者目标原理雷达和靠近雷达的速度方向是不同的，一般在连续波雷达中目标靠近雷达是速度是负的，而远离雷达目标速度是正的，这有点违反了人的直觉，但实事就是如此。关于这个现象的解释，我在这里做了个人的分析：

调皮连续波：雷达原理 | 讨论调频连续波雷达目标运动方向与速度正负的关系？

t = 0:1/fs:Tr-(1/fs); %chirp采样的时间序列
for chirpId = 1:chirps
   for txId = 1:txNum 
        St = exp((1i*2*pi)*(centerFreq*(t+( chirpId-1)*Tr)+slope/2*t.^2)); %发射信号
        for rxId = 1:rxNum
            Sif = zeros(1,rangeBin);
            for targetId = 1:targetNum

                %%连续帧 目标设置，如果不需要连续帧，令Parameter.frame=0，即可。
                if targetId==1
                    targetRange = target(targetId,1)-Parameter.frame; 
                    targetSpeed = target(targetId,2); 
                    targetAngle = target(targetId,3);
                elseif targetId==2
                    targetRange = target(targetId,1)+0.5*Parameter.frame; 
                    targetSpeed = target(targetId,2); 
                    targetAngle = target(targetId,3);
                elseif targetId==3
                    targetRange = target(targetId,1)+Parameter.frame; 
                    targetSpeed = target(targetId,2); 
                    targetAngle = target(targetId,3);
                end

                tau = 2 * (targetRange + targetSpeed * (txId - 1) * Tr) / c;
                fd = 2 * targetSpeed / lambda;
                wx = ((txId-1) * rxNum + rxId) / lambda * dx * sind(targetAngle);
                Sr = 10*exp((1i*2*pi)*((centerFreq-fd)*(t-tau+( chirpId-1) * Tr)+slope/2*(t-tau).^2 + wx));  %回波信号
                Sif = Sif + St .* conj(Sr);
                %叠加20dB高斯白噪声
                Sif = awgn(Sif,20);
            end
            rawData((txId-1) * rxNum + rxId,:,chirpId) = Sif;
        end
    end
end

三、回波信号分析

回波信号分析这里，我们其实可以假设，原始数据就是来自雷达射频前端，经过混频器、带通或者低通滤波器，然后被ADC正交采样。雷达其实可以采用正交采样（复采样），也可以采用实采样，二者各自有优缺点，选择什么却决于雷达系统设计人员的设计，以及实际的需求，但本文还是选择正交采样。关于这点，可以看这篇文章：

调皮连续波：FMCW雷达系统中的复基带架构（中英文翻译对照）

这里我们采用一个chirp来验证我们的想法，假设接收到的信号是IQ不平衡的，我们拆分IQ的实数部分和虚数部分，然后把实数部分增加幅度因子和相位因子，人为使得IQ不平衡。

%% IQ通道校正
firstChirp_I = imag(firstChirp);
firstChirp_Q = real(firstChirp);
%设定幅度误差因子
alpha = 2;
%设定相位误差因子
phi = 5;
firstChirp_Q_1 = (alpha+1)firstChirp_Qexp(phi2pi/360);

然后采用IQ不平衡补偿算法来修正，具体的算法原理可以看这篇文章：

调皮连续波：雷达信号处理中I/Q正交信号不平衡补偿算法（含matlab代码）（文章中有个公式和代码写错了，注意区分，在本文中的代码修改过来了）。本文一谈到原理，我就叫大家去看别的文章，其实很多基础知识之前都分享个人的见解，在这里方便引用，就不再重复说了。

%% IQ通道不平衡效果
firstChirp_IQ_1 = complex(firstChirp_Q_1,firstChirp_I);
%IQ通道校正算法
%求均值
firstChirp_I = imag(firstChirp_IQ_1);
firstChirp_Q = real(firstChirp_IQ_1);
I_before_correction =firstChirp_I-mean(firstChirp_I);
Q_before_correction=firstChirp_Q -mean(firstChirp_Q);
%估计参数
alphi = sqrt(mean(Q_before_correction.*Q_before_correction)/mean(I_before_correction.*I_before_correction))-1;
phi = -asin(mean(I_before_correction.*Q_before_correction)/sqrt(mean(I_before_correction.*I_before_correction)*mean(Q_before_correction.Q_before_correction)));
%P矩阵求解
P = [1,0;tan(phi),1/((1+alphi)cos(phi))];
%计算IQ
IQ = P[I_before_correction;Q_before_correction];
%重组信号
I =IQ(1,:);
Q =IQ(2,:);
signal_IQ =Q+I1j;
%图形绘制
% figure(2)
% subplot(2,1,1);
% plot((abs(fft(firstChirp_IQ_1))));
% xlabel(‘距离(m)’); ylabel(‘幅值’);title(‘距离维FFT(校正前)’);
% subplot(2,1,2);
% plot((abs(fft(signal_IQ))));
% xlabel(‘距离(m)’); ylabel(‘幅值’);title(‘距离维FFT(校正后)’);

IQ补偿前后的效果如下图所示，可见镜像频率还是很明显的，抑制后频谱干净了许多。本文采取的是窄带IQ不平衡补偿，因为中频信号的带宽比较窄，如果是换做通信系统中，则是需要宽带IQ不平衡补偿算法，那时会采取一种滤波器组的方法实现。

在这里插入图片描述
时域信号如下：

在这里插入图片描述

四、距离估计

现在要解算目标的距离了，原理很简单，之前的文章看过就清楚了，无非就是对一个chirp内的中频信号做一次FFT，这里以单个chirp为例子说明：

%%1D FFT
fft1dData = fft(firstChirp);
% figure(3);
% plot(rangeIndex,db(abs(fft1dData)./max(abs(fft1dData))));
% xlabel(‘距离(m)’); ylabel(‘幅值(dB)’);title(‘距离维FFT’);

在这里插入图片描述

五、速度估计

计算速度和计算距离的顺序可以交替，就像我在之前的文章（知乎答疑 | 雷达信号处理中的距离维FFT和速度维FFT的执行顺序能够交换嘛？）中谈到的那样。

这里需要对所有接收通道都做FFT，其实如果时单纯计算目标的速度和距离的话，利用一个通道就可以了，但是为了后续进行CFAR检测之前做非相干积累，这里必须多所有通道处理。

%% 2D FFT
%% 距离-多普勒谱
channelNum = size(rawData,1);
rangebinNum = size(rawData,2);
dopplerbinNum = size(rawData,3);
fft2dDataPower= zeros(size(rawData));
fft2dDataDB = zeros(size(rawData));
fftRADataPower= zeros(size(rawData));
for chanId = 1:1:channelNum
fft2dDataPower(chanId,:, : ) = RDfftMatrix(rawData(chanId,:,: ));
end

% figure(4);
% mesh(dopplerIndex’,rangeIndex,db(abs(squeeze(fft2dDataPower(chanId,:,: )))));
% view(2);
% xlabel(‘速度(m/s)’); ylabel(‘距离(m)’); zlabel(‘幅值’);
% title(‘距离-多普勒谱’);
% mesh(abs(squeeze(fft2dDataPower(chanId,:,: ))));

速度估计得到的时距离-多普勒谱，也就是RD图。这里距离和速度估计统一做FFT，代码如下，并且还采用了hanning窗加窗处理（加权）。关于雷达信号处理加窗的详细原理，在文章（雷达信号处理中的离散傅里叶变换（DFT）以及加窗）中有说明。

rawData = squeeze(rawData);
[rangeBin,dopplerBin] = size(rawData);
rangeWin = hanning(rangeBin);
rangeWin2D = repmat(rangeWin,1,dopplerBin);
dopplerWin = hanning(dopplerBin)';
dopplerWin2D = repmat(dopplerWin,rangeBin,1);
rawDataWin = rawData .* rangeWin2D;
fft1dData = fft(rawDataWin,rangeBin,1);
fft1dDataWin = fft1dData .* dopplerWin2D;
fft2dData = fftshift(fft(fft1dDataWin,dopplerBin,2),2);

距离处理不需要fftshift，速度需要fftshift处理把零速度通道（零频分量）移到中点，重新排列傅里叶变换。

在这里插入图片描述

六、角度估计

角度估计和距离估计也是放在一起的，这里为了方便体现二者之间的关系。步骤与距离和速度估计一样，代码如下，也是采用了hanning窗加窗，窗函数可以选择其他的，比如泰勒窗、布莱克曼窗等等，具体看信噪比以及其他方面的需求。

rawData = squeeze(rawData);
[angleBin,rangeBin] = size(rawData);

angleWin = hanning(angleBin);
angleWin2D = repmat(angleWin,1,rangeBin);

rangeWin = hanning(rangeBin)';
rangeWin2D = repmat(rangeWin,angleBin,1);

rawDataWin = rawData .* angleWin2D;
fft1dData = fftshift(fft(rawDataWin,angleBin,1));

fft1dDataWin = fft1dData .* rangeWin2D;
fft2dData = fft(fft1dDataWin,rangeBin,2);

计算结果如下图所示，符合目标参数中的角度设置，但是角度还是有些不太精确，还要看后续进行多普勒补偿后再看看结果如何。

parameter.target = [100 -20 0; %target1 range speed angle0 10 -30; %target2 range speed angle0 20 30; %target2 range speed angle ];

在这里插入图片描述

七、非相干累积

相干累积速度估计时做FFT时体现的，利用多个chirp做FFT就是相干积累，非相干积累就是纯粹的幅值叠加，二者之间的信噪比改善是不同的。

[channelNum,rangeBin,dopplerBin] = size(fft2dDataDB);
accumulateRD = zeros(rangeBin,dopplerBin);

for channelId = 1:channelNum
    accumulateRD = accumulateRD + (abs(squeeze(fft2dDataDB(channelId,:,:))));
end

在这里插入图片描述

八、CFAR检测

CFAR检测需要预先根据目标特性设置CFAR参数，如参考单元、保护单元、阈值，以及CFAR模式选择（CA-CFAR、SO-CFAR、GO-CFAR、OS-CFAR）等。

parameter.dopplerMethod = 1; %1：ca-cfar 2：so-cfar 3：go-cfar
parameter.dopplerSNR = 10;
parameter.dopplerWinGuardLen = 2;
parameter.dopplerWinTrainLen = 8;

parameter.rangeMethod = 1; %1：ca-cfar 2：so-cfar 3：go-cfar
parameter.rangeSNR = 10;
parameter.rangeWinGuardLen = 2;
parameter.rangeWinTrainLen = 4;

CFAR检测采用两次一维CFAR组成二维CFAR，首先再速度为进行一次CA-CFAR，然后再距离维再进行一次CA-CFAR，这样可以减少计算量。

上述中CFAR阈值也被称为门限因子，通常门限因子是10进制的常数： $\geqslant K_0 Z$ ，Z是左右参考单元平均后的平均值，K0是常数。

但是某些情况下也需要转为对数（dB），转换公式为：

$\log (Z)+10 * \log (K 0)$

代码中设置为T=10dB（（程序中是parameter.dopplerSNR））:

$\left.T=20 * \log 10(T 0), \quad T 0=10^{(} T / 20\right), T 2=\log 2(T 0)=\log 2\left(10^{\prime} T / 20\right)=T / 6$

在MATLAB中采用log10()好计算，但是在硬件平台采用log2()更好，因为计算机都是二进制，计算效率高一些。TI的程序代码中CFAR阈值设置的是一个十进制数（Tcli =5120），转换为对数后T=15dB，转化公式为：

$T_{c l i}=256 \times$ numVirtualAntennas $\times \mathrm{T}(\mathrm{dB}) / 6$

CFAR检测的结果如下图所示：

在这里插入图片描述

九、峰值聚合

峰值聚合的目的是为了对CFAR检测之后的目标点进行一次精细搜索，本次仿真所涉及到的原理如下图所示。

在这里插入图片描述
（峰值聚合模型）

对照代码分析，简单一句话就是判断某一个点是否比周围四个点大，如果大于，就是需要的点，如果不是那就舍弃，如代码中第二个if。

peakSearchList = [];
[rangeLen, dopplerLen] = size(RD_cfar);
RD_pearkSearch = zeros(rangeLen, dopplerLen);
length = size(cfarTargetList,2);

for targetIdx = 1:length

    rangeIdx   = cfarTargetList(1,targetIdx); 
    dopplerIdx = cfarTargetList(2,targetIdx); %坐标

    if rangeIdx > 1 && rangeIdx < rangeLen && dopplerIdx > 1 && dopplerIdx < dopplerLen %边界点不考虑  
       
        if RD_cfar(rangeIdx,dopplerIdx) > RD_cfar(rangeIdx - 1,dopplerIdx) && ...
                RD_cfar(rangeIdx,dopplerIdx) > RD_cfar(rangeIdx + 1,dopplerIdx) && ...
                RD_cfar(rangeIdx,dopplerIdx) > RD_cfar(rangeIdx,dopplerIdx - 1) && ...
                RD_cfar(rangeIdx,dopplerIdx) > RD_cfar(rangeIdx,dopplerIdx + 1)

                RD_pearkSearch(rangeIdx,dopplerIdx) = RD_cfar(rangeIdx,dopplerIdx);

                cfarTarget = [rangeIdx ; dopplerIdx];

                peakSearchList = [peakSearchList cfarTarget];
        end   
    end
end

峰值聚合后的结果如下图所示：

在这里插入图片描述

峰值聚合的算法不止上面这种，另外还有很多，感兴趣的同学可以自己检索。

十、多普勒补偿

之前的角度估计是利用3D-FFT来实现的，由于FMCW是利用相邻chirp间的相位差的变化来估计目标速度的，并没有考虑运动目标的多普勒频移，根据接收信号模型：

Sr = 10exp((1i2pi)((centerFreq-fd)(t-tau+( chirpId-1) * Tr)+slope/2(t-tau).^2 + wx)); %回波信号

假设，雷达在2T4R TDM-MIMO的工作方式下，TX2形成的4个虚拟天线会由于运动目标的多普勒效应相对于TX1形成的4个虚拟天线将产生一个多普勒相位偏移，这种多普勒偏移引起的相位差会映射到角度上，引起测角不准确。因此，对目标进行角度估计前须对相位偏移进行校正，也称作多普勒相位补偿，保证角度估计的准确性。

在这里插入图片描述
（多普勒偏移）

关于多普勒相位补偿和速度扩展算法的原理，可以阅读TI官方的文档：基于AWR1642汽车雷达的速度扩展算法研究.pdf。多普勒补偿的程序也比较简单，主要就是利用之前计算出来的速度值，把相位偏移计算出来：

$\Delta \varphi=2 \pi \Delta f T_c$

其中 $\Delta f$ 为多普勒频偏， $T_c$ 为chirp调制时间，上述公式在代码中被转化为单纯的多普勒单元做运算了，具体转换公式如下所示：

$\Delta \varphi=2 \pi \Delta f T_c=2 \pi \frac{2 V_r}{\lambda} T_c=2 \pi \frac{2 T_c}{\lambda} \Delta V_r * n=2 \pi \frac{2 T_c}{\lambda} * \frac{\lambda}{2 N * T_c} * n=2 \pi * \frac{n}{N}$

其中，n就是当前速度所在的多普勒单元索引，N就是全部多普勒单元索引，然后再补偿回去就ok了，具体看下面的程序就一目了然。

compCoffVec = zeros(parameter.virtualAntenna,1);
txNum = length(parameter.txAntenna);
rxNum = length(parameter.rxAntenna);
phi = 2 * pi * (speedBin - parameter.dopplerBin/2 - 1) / parameter.dopplerBin;
delta = phi 
for txId = 1:txNum
    for rxId = 1:rxNum
        compCoffVec((txId-1) * rxNum + rxId) = exp(-1i * (txId-1) * delta);
    end
end

经过多普勒补偿前后的角度结果对比如下图所示，上图为整体角度差异，下图为局部角度差异，可见多普勒补偿后角度更加接近实际值：

在这里插入图片描述

（局部图，左边为未补偿，右边为补偿后）

十一、DOA估计

终于到最后一节了，不容易啊，经过多普勒补偿之后角度会更加准确，在仿真程序中采用了三种DOA估计算法，分别是FFT、DBF以及Capon，三种算法都可以选择。其中DBF关键在于计算加权系数，然后加权求和：

在这里插入图片描述

for degscan = deg
    for txId = 1:txNum
        for rxId = 1:rxNum
            dphi = ((txId-1) * rxNum + rxId - 1) * 2 * pi / lambda * dx * sind(degscan);
            weightVec((txId-1) * rxNum + rxId) = exp(-1i * dphi);
        end
    end
    doa_dbf(kk) = antVec'*weightVec;
    kk = kk + 1;
end
doa_abs = abs(doa_dbf);
[pk,loc]=max(doa_abs);
angle = deg(loc);

FFT算法跟之前一样，比较简单：

angleIndex = asin((-512:1:512-1)/512) * 180 / pi;
doa_fft=fftshift(fft(antVec,1024));
doa_abs=abs(doa_fft);
[pk,loc]=max(doa_abs);
angle = angleIndex(loc);

Capon公式稍微多一些，但也不复杂，三个算法本质上都是空间滤波。

kk = 1;
for degscan = deg
    for txId = 1:txNum
        for rxId = 1:rxNum
            virtualAntennaId = (txId-1) * rxNum + rxId - 1;
            dphi = 2 * pi / lambda * dx * virtualAntennaId * sind(degscan);
            a((txId-1) * rxNum + rxId) = exp(-1i * dphi);
        end
    end
    RxInv = inv(Rx);
    P_out(kk) = 1/(a'*RxInv*a);
    kk = kk + 1;
end

DOA估计部分，读者可以自由发挥，各种超分辨算法（如MUSIC）、快速单快拍FFT DOA都可以加进来仿真看效果，这里给予大家的自由度是非常高的。

参考资料

【1】Analysis and Comparison of MIMO Radar Waveforms
【2】https://blog.csdn.net/Xiao_Jie123/article/details/126538880

结语

本文花费了大量时间来编程序、读程序、写文章，但是为了感谢各位读者长期以来的关注和支持，所以就一瓶冰红茶钱，基于本文的程序，各位读者还可以将TDMA-MIMO模式修改为DDMA-MIMO模式，参考TI的文档实现DDMA-MIMO加入空带（EmptyBand）解速度模糊，挺好玩的。另外很多算法、方案都可以基于本程序进行探索，祝愿大家一切顺利！

在这里插入图片描述
（DDMA-MIMO步进相位）