登堂入室：毫米波雷达开发手册之信号模型
严肃声明
严格禁止未经本人允许转载本文工作成果（包括引言、证明、代码分析等），一经发现被用于学术或商业用途，将予以法律警告。欢迎和笔者深度合作，探讨学术话题。
本人发现部分网站（如https://it.cha138.com/等野鸡网站）未经本人允许私自将本人博客内容转发公开，本人严禁杜绝此类行为，一经发现将投诉举报相关行为并追究违法责任！！！
写在前面
深知新手在接触毫米波雷达板硬件时需要花费的沉没成本，因此在行将告别毫米波雷达之际，总结这两年以来在毫米波雷达上的一些经验和教训。
本文档用于为实现基于AWR1243BOOST等单板毫米波雷达开发提供参考指南与解决方案，主要包括硬件配置、基础参数、信号模型、应用DEMO开发以及可深入研究方向思考等；为更好地匹配后续级联雷达应用的学习路线，在本手册中会尽可能同化单板雷达和级联雷达中的相关表述。
本指南作者信息：Xl，联系方式：12331100@zju.edu.cn。未经本人允许，请勿用于商业和学术用途。
希望后者在使用本指南时可以考虑引用作者在毫米波雷达旅途中的相关工作,如下所列：
[1] X. Yu, Z. Cao, Z. Wu, C. Song, J. Zhu and Z. Xu, "A Novel Potential Drowning Detection System Based on Millimeter-Wave Radar," 2022 17th International Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, Singapore, 2022, pp. 659-664, doi: 10.1109/ICARCV57592.2022.10004245.
[2] X. Yu, Z. Cao, Z. Wu, C. Song and Z. Xu, "Sample Intercorrelation-Based Multidomain Fusion Network for Aquatic Human Activity Recognition Using Millimeter-Wave Radar," in IEEE Geoscience and Remote Sensing Letters, vol. 20, pp. 1-5, 2023, Art no. 3505205, doi: 10.1109/LGRS.2023.3284395.
[3] Z. Wu, Z. Cao, X. Yu, J. Zhu, C. Song and Z. Xu, "A Novel Multiperson Activity Recognition Algorithm Based on Point Clouds Measured by Millimeter-Wave MIMO Radar," in IEEE Sensors Journal, vol. 23, no. 17, pp. 19509-19523, 1 Sept.1, 2023, doi: 10.1109/JSEN.2023.3283778.
[4] Z. Cao, Z. Wu, X. Yu, C. Song and Z. Xu, "People Flow Detection Algorithm Based on a Multiencoder-Classifier Cotraining Architecture for FMCW Radar," in IEEE Transactions on Geoscience and Remote Sensing, vol. 61, pp. 1-18, 2023, Art no. 5108418, doi: 10.1109/TGRS.2023.3305551.
本章节为信号模型章节，主要解读时域信号模型、空域信号模型和信号处理栈。
欢迎各位读者通过邮件形式与笔者交流讨论，本章节完整程序请私信笔者，希望使用本代码时能够提供一份引用和Star，以表示对笔者工作的尊重，谢谢！在后续将定时维护更新。
https://github.com/DingdongD/TDMA-MIMO
往期内容：
初出茅庐：毫米波雷达开发手册之基础参数
扬帆起航：毫米波雷达开发手册之硬件配置
眼观四海：自动驾驶&4D成像毫米波雷达 如今几何？
信号模型分析
时域信号模型
如图所示，雷达收发器（TRX）目前通常被实现为单个单片微波集成电路（MMIC）或甚至包括数字功能的片上系统（SoC），其生成雷达发射（TX）波形并经由连接到其TX信道的天线辐射该波形。电磁波通过空间传播，直到它命中一个或多个目标。假设目标的各向同性辐射，与目标雷达截面（RCS）成比例的入射波通量密度的一部分被捕获并朝向雷达再辐射。根据经典理论，目标通常被假设为孤立的点状散射体即点目标，如汽车、行人或环境之类的物理对象由若干点状目标组成，并且各个目标响应的总和则是总体接收（RX）信号。
通过比较雷达的TX和RX信号，可以估计出感知任务的基本参数：目标数目\(K_{tgt}\)和每个目标\(k_{tgt}\)的目标距离\(r_{k_{tgt}}\)，目标径向速度\(v_{r_{k_{tgt}}}\)，方位角\(\phi_{k}\)和俯仰角\(\epsilon_k\)。在每个测量周期中检测到的目标集合被称为目标列表，并且可以被解释为从雷达测量获得的点云。后续的毫米波雷达信号处理步骤主要有是聚类（即识别属于同一物理目标）、跟踪（即从目标列表的序列中导出物理上存在的对象的信息和统计数据）。本章节的核心是介绍TDMA-MIMO框架下的信号模型仿真及实现方法。
本文所研究的模型与快速线性调频调频连续波（FMCW）雷达相关，这是目前智能驾驶、智能家居领域中所使用的主要雷达波形。如下图所示为典型的由Q个chirp组成的chirp序列构成的发射波形（绿色部分）和来自单个目标的接受波形信号（蓝色部分）。
发射信号\(s（t）\)由以chirp-chirp间隔的为\(T_{cc}\)的\(Q\)个线性调频脉冲波形\(s_c（t）\)的序列给出：
\[s(t)=\sum\limits_{q=1}^Q s_{\mathrm{c}}(t-q T_{\mathrm{cc}})
\]
其中脉冲持续时间为\(T_c\)的每个chirp的信号表示为：
\[s_\mathrm{c}(t)=e^{j(2\pi f_1t+\pi\mu t^2+\varphi_0)}\mathrm{rect}\left(\frac{t-T_\mathrm{c}/2}{T_\mathrm{c}}\right).
\]
其中，chirp斜率为\(\mu\)，起始频率为\(f_1\)，截止频率为\(f_h=f_l+W\)，\(W=\mu T_c\)表示chirp的带宽，chirp的中心频率\(f_c\)定义为\(f_1+W/2\)，一帧chirp序列的长度或称帧周期为\(T_{seq}=QT_{cc}-(T_{cc}-T_c)\)。
每个点目标作为发射信号\(s(t)\)的时延副本以加权形式参与到接收信号\(r(t)\)中，其表示形式如下：
\[r(t)=\sum_{k_{t g t}=1}^{K_{t g t}}a_{k_{t g t}}s(t-T_{\mathrm{d},k}(t))
\]
其中，目标衰减系数\(a_{k_{tgt}}\)是由于自由空间传播和目标RCS引起的系数加权，往返时间\(T_{d,k}(t)\)是时间相关的，因为这考虑了目标的移动，通常帧周期为几毫秒或更小，与目标速度相比较短，无法形成积累效应。因此，\(T_{\mathrm{d},k}(t))\)可以近似表示为：
\[T_{\mathrm{d},k}(t)=\frac{2r_{k_\mathrm{tgt}}+2v_{r,k_\mathrm{tgt}}t}{c_0}.
\]
FMCW雷达的基本原理是，发射信号\(s（t）\)和接收信号\(r（t）\)相乘混频，混频后通过低通滤波器LPF以滤除在和频（约两倍载波频率）处振荡的和分量，保留低通滤波信号即中频信号IF Signal，然后与发射信号定时同步地ADC采样。对于每个chirp，ADC采样将生成P个样本的向量，这意味着\(Q\)个线性调频脉冲的采样将生成\(Q\)个向量，并且每个向量长度为\(P\)，以二维矩阵\({\mathcal{B}}=(b_{p,q}),{\mathcal{B}}\in{{\mathbb{C}}^{P\times Q}}\)的形式存储，每个元素为\(b_{p,q}\)。\(p\)所在的维度称为快时间维，\(q\)所在的维度称为慢时间维度，这在上一节内容中已经提及。
基于所给的调制参数，可以将二维矩阵\({\mathcal{B}}\)的每个元素\(b_{p,q}\)展开如下：
\[b[p,q]\approx\sum_{k_1\mathfrak{g}=1}^{k_1}a_{k_1\mathfrak{g}n}e^{j(\lambda_{k_1\mathfrak{g}n}p+\zeta_{k_1\mathfrak{g}n}q)}e^{j\varphi_{k_1\mathfrak{g}n}}+\omega[p,q].
\]
其中，对应距离维度的归一化角频率可表示如下：
\[\lambda_{k_\mathrm{tgt}}=2\pi\left(\frac{2\mu r_{k_\mathrm{tgt}}}{c_0}+f_{\mathrm{d},k}\right)T_\mathrm{s},
\]
对应多普勒维度的归一化角频率可表示如下：
\[\xi_{k_{\mathrm{tgt}}}=2\pi\left(\frac{2f_{\mathrm{c}}v_{r,k_{\mathrm{tgt}}}}{c_{0}}\right)T_{\mathrm{cc}}
\]
目标的初始相位则可表示如下：
\[\varphi_{k_\mathrm{tgt}}=\frac{2f_c r_{k_\mathrm{tgt}}}{c_0}.
\]
空域信号模型
为了建立与目标方位角和仰角的关系，可以使用\(N_{RX}\)个接收天线的阵列，这将被称为单输入多输出（SIMO）配置；为了进一步提高角分辨率，在共同定位的相干MIMO配置中的若干发射和接收天线可用于现代毫米波雷达中，形成\(N=N_{TX}N_{RX}\)的虚拟接收阵列。
下图提供了一种基于三发四收的物理天线配置可能和MIMO后的虚拟阵列效果，左边天线标记为接收天线，右边天线标记为发射天线，可以很明显地看到TX和TX1/TX2在俯仰维度上存在偏移量，RX2和其他三根接收天线在俯仰维度上存在偏移量。
为更好地理解MIMO机制，下面将手把手推导如何得到虚拟阵列图(假设每个方格代表阵元间距d)。
-
TX1在四根接收天线上形成的水平相位为[0 w 2w 3w]，俯仰相位为[0 w 0 0] -
TX2在四根接收天线上形成的水平相位为[w 2w 3w 4w]，俯仰相位为[0 w 0 0] -
TX3在四根接收天线上形成的水平相位为[2w 3w 4w 5w]，俯仰相位为[2w 3w 2w 2w] -
匹配上面的推导和下图中的右下边虚拟阵列，可以发现结果一致。这种推导在小阵列情况下容易实现，但是随着阵元数目的增加则容易出错，这是因为可能存在冗余阵元，如
TI AWR2243Cascade具有12个发射天线和16个接收天线，实际在方位维度仅有86个虚拟阵元。为了排列准确，笔者写了以下程序来实现物理阵元到虚拟阵列的精确映射。function cfgOut = ConfigureRTX() % 下面这个程序实现了不同级联雷达的物理发射阵元和接收阵元配置，并考虑了天线发射顺序（这对虚拟阵元的排列影响不大，但会影响adc数据中的数据排列）。 % 实际的阵元间距可能与下面的存在误差，以下板子参数多为本人估算所得。 % 具体代码结合注释和上下文语句进行分析，这里不再赘述。 % By Xuliang cfgOut.Mode = 4; % 选择需要仿真的雷达模式 这里提供了6种毫米波雷达配置，分别为TI if cfgOut.Mode == 1 disp(strcat(['=====','已选择TI AWR1243-单板雷达模式','====='])); % 单板雷达模式3T4R cfgOut.numTx = 3; cfgOut.numRx = 4; cfgOut.PosTX_X = [0 2 4]; % 发射天线水平方向阵元间距 cfgOut.PosTX_Y = [0 1 0]; % 发射天线垂直方向阵元间距 cfgOut.PosTX_BOARD_ID = [1 2 3]; % 发射天线板上顺序 后面没用到 cfgOut.PosTX_Trans_ID = [1 3 2]; % 发射天线发射顺序 cfgOut.PosRX_X = [0 1 2 3]; % 接收天线水平方向阵元间距 cfgOut.PosRX_Y = [0 0 0 0]; % 接收天线垂直方向阵元间距 cfgOut.PosRX_BOARD_ID = [1 2 3 4]; % 接收天线板上顺序 elseif cfgOut.Mode == 2 disp(strcat(['=====','已选择TI AWR2243-级联雷达模式','====='])); % 级联雷达模式12T16R cfgOut.numTx = 12; cfgOut.numRx = 16; cfgOut.PosTX_X = [11 10 9 32 28 24 20 16 12 8 4 0]; cfgOut.PosTX_Y = [6 4 1 0 0 0 0 0 0 0 0 0]; cfgOut.PosTX_BOARD_ID = [12 11 10 3 2 1 9 8 7 6 5 4]; cfgOut.PosTX_Trans_ID = [12:-1:1]; cfgOut.PosRX_X = [11 12 13 14 50 51 52 53 46 47 48 49 0 1 2 3]; cfgOut.PosRX_Y = zeros(1,cfgOut.numRx); cfgOut.PosRX_BOARD_ID = [13 14 15 16 1 2 3 4 9 10 11 12 5 6 7 8]; elseif cfgOut.Mode == 3 % 这个模式为特斯拉6T8R级联板 disp(strcat(['=====','已选择特斯拉6T8R-级联雷达模式','====='])); cfgOut.numTx = 6; cfgOut.numRx = 8; cfgOut.PosTX_X = [0 17 34 41 48 55]; cfgOut.PosTX_Y = [0 0 0 0 4 8 12]; cfgOut.PosTX_BOARD_ID = [1 2 3 4 5 6]; cfgOut.PosTX_Trans_ID = [1 2 3 4 5 6]; cfgOut.PosRX_X = [0 2 3 5 8 11 14 17]; cfgOut.PosRX_Y = zeros(1,cfgOut.numRx); cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; elseif cfgOut.Mode == 4 % 这个模式为加特兰 4T4R单板 disp(strcat(['=====','已选择加特兰4T4R-单板雷达模式','====='])); cfgOut.numTx = 4; cfgOut.numRx = 4; cfgOut.PosTX_X = [0 5 12 18]; cfgOut.PosTX_Y = [5 0 0 0]; cfgOut.PosTX_BOARD_ID = [1 2 3 4]; cfgOut.PosTX_Trans_ID = [4 3 2 1]; cfgOut.PosRX_X = [0 2 5 7]; cfgOut.PosRX_Y = [0 0 0 2]; cfgOut.PosRX_BOARD_ID = [1 2 3 4]; elseif cfgOut.Mode == 5 % 这个模式为福瑞泰克 6T8R级联板 disp(strcat(['=====','已选择福瑞泰克6T8R-级联雷达模式','====='])); cfgOut.numTx = 6; cfgOut.numRx = 8; cfgOut.PosTX_X = [0 2 4 6 8 10]; cfgOut.PosTX_Y = [0 0 0 0 0 0]; cfgOut.PosTX_BOARD_ID = [1 2 3 4 5 6]; cfgOut.PosTX_Trans_ID = [6 5 4 3 2 1]; cfgOut.PosRX_X = [0 2 5 13 13+ones(1,4)*1] + [0 0 0 0 0 2 5 13]; % 暂时不确定双级联间距 cfgOut.PosRX_Y = [1 0 0 0 1 0 0 0 0]; cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; elseif cfgOut.Mode == 6 % 这个模式为TI双级联毫米波雷达 disp(strcat(['=====','已选择TI 6T8R-级联雷达模式','====='])); cfgOut.numTx = 6; cfgOut.numRx = 8; cfgOut.PosTX_X = [0 2 4 7 11 15]; cfgOut.PosTX_Y = [0 1.6 0 0 0 0]; cfgOut.PosTX_BOARD_ID = [1 3 2 4 5 6]; cfgOut.PosTX_Trans_ID = [6 5 4 3 2 1]; cfgOut.PosRX_X = [0 1 2 3 19 20 21 22]; cfgOut.PosRX_Y = [0 0 0 0 0 0 0 0 0]; cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; end virtual_azi_arr = repmat(cfgOut.PosTX_X(cfgOut.PosTX_Trans_ID), cfgOut.numRx, 1) + repmat(cfgOut.PosRX_X(cfgOut.PosRX_BOARD_ID), cfgOut.numTx, 1).'; virtual_ele_arr = repmat(cfgOut.PosTX_Y(cfgOut.PosTX_Trans_ID), cfgOut.numRx, 1) + repmat(cfgOut.PosRX_Y(cfgOut.PosRX_BOARD_ID), cfgOut.numTx, 1).'; virtual_array.azi_arr = reshape(virtual_azi_arr,[],1); virtual_array.ele_arr = reshape(virtual_ele_arr,[],1); virtual_array.tx_id = reshape(repmat(cfgOut.PosTX_Trans_ID, cfgOut.numRx, 1),[],1); virtual_array.rx_id = reshape(repmat(cfgOut.PosRX_BOARD_ID, cfgOut.numTx, 1).',[],1); virtual_array.virtual_arr = cat(2,virtual_array.azi_arr,virtual_array.ele_arr,virtual_array.rx_id,virtual_array.tx_id); cfgOut.virtual_array = virtual_array; % 将virtual_array参数存入更新 end% 测试文件：可视化虚拟阵元排列 % By Xuliang clc;clear;close all cfgOut = ConfigureRTX(); plot(cfgOut.virtual_array.virtual_arr(:,1), cfgOut.virtual_array.virtual_arr(:,2),'ro'); xlabel('方位维');ylabel('俯仰维')；grid minor； % 这里不再进行可视化。
在上一章节中，我们可以得到单信道下的时域信号模型为
\(b[p,q]\approx\sum_{k_1\mathfrak{g}=1}^{k_1}a_{k_1\mathfrak{g}n}e^{j(\lambda_{k_1\mathfrak{g}n}p+\zeta_{k_1\mathfrak{g}n}q)}e^{j\varphi_{k_1\mathfrak{g}n}}+\omega[p,q].\)$
假设目标满足远场窄带信号假设，那么可进一步扩展时域信号至空域信号如下：
\[b[p,q,n]\approx\sum_{k\text{g}s=1}^{K_\text{g}t}a_{k\text{g}t}e^{i(\lambda_{k\text{g}t}p+\zeta_{k\text{g}s}q)}e^{i\varphi_{k\text{g}t}}e^{j\frac{2\pi}{\lambda}(\mathfrak{u}_{k\text{g}s}\mathfrak{p}_n)}+\omega[p,q,n].
\]
上述信号对应于雷达3D tensor或radar cube数据的元素，根据ISO8000:2标准定义的球面坐标系，目标\(k_{tgt}\)对应的单位长度波达角矢量可以表示如下：
\[\mathbf{u}_k=\mathbf{r}_{k_\textrm{tgt}}/\|\mathbf{r}_{k_\textrm{tgt}}\|=\begin{bmatrix}\cos\phi_{k_\textrm{tgt}}\sin\theta_{k_\textrm{t}}\\ \sin\phi_{k_\textrm{gt}}\sin\theta_{k_\textrm{tgt}}\\ \cos\theta_{k_\textrm{tgt}}\end{bmatrix}
\]
其中，\(\phi_{k_\textrm{tgt}}\)和\(\theta_{k_\textrm{tgt}}\)分别表示目标的方位角和倾斜角，俯仰角则表示为\(\epsilon_{k_{tgt}}=\pi-\theta_{k_{tgt}}\)；虚拟天线的位置\(p_n\)可以表示如下（\(n=0....N-1\)）
\[\mathbf{p}_n=\mathbf{p}_{\mathrm{RX},n_{\mathrm{RX}}}+\mathbf{p}_{\mathrm{TX},n_{\mathrm{TX}}},n_{\mathrm{RX}}=n\bmod N_{\mathrm{RX}},\quad n_{\mathrm{TX}}=n\bmod N_{\mathrm{TX}}
\]
至此，我们可以很明确地建立目标角度和模型\(e^{(2\pi/\lambda)\mathbf{u}_{k_\mathrm{tgt}}\mathbf{p}_n}\)有关，上面提到的目标初始相位\(\varphi_{k_\mathrm{tgt}}=\frac{2f_c r_{k_\mathrm{tgt}}}{c_0}\)可以衡量目标在参考坐标系原点处出现的相位，不大影响角度的估计。
基于以上信号模型，我们可以生成TDMAI-MIMO下的radar cube数据文件，核心代码如下所示：
function [signal] = GenerateSigIQ(tarOut, cfgOut)
%% 本文件用于毫米波雷达的ADC-IQ信号
%% By Xuliang
tarNum = tarOut.nums; % 目标数目
numTx = cfgOut.numTx; % 发射天线数目
numRx = cfgOut.numRx; % 接收天线数目
ADCNum = cfgOut.ADCNum; % ADC采样数目
Frame = cfgOut.Frame; % 帧数
ChirpNum = cfgOut.ChirpNum; % 每帧发射Chirp数目
TotalChirpNum = Frame * ChirpNum; % 总Chirp数目
% 信号参数
fc = cfgOut.fc; % 载频 Hz
fs = cfgOut.fs; % ADC采样频率 Hz
Ramptime = cfgOut.Ramptime; % 工作时间
Idletime = cfgOut.Idletime; % 空闲时间
Slope = cfgOut.Slope; % Chirp斜率
validB = cfgOut.validB; % 实际有效带宽
totalB = cfgOut.totalB; % 完整带宽
Tc = cfgOut.Tc; % 单Chirp周期
ValidTc = cfgOut.ValidTc; % 单Chirp内ADC有效采样时间
% 硬件参数
Pt = cfgOut.Pt; % dbm 发射功率
Fn = cfgOut.Fn; % 噪声系数
Ls = cfgOut.Ls; % 系统损耗
% 天线阵列
antennaPhase = cfgOut.antennaPhase; % 天线相位
virtual_array = cfgOut.virtual_array; % 虚拟阵列 这里的虚拟阵列是struct virtual_arr中包含了发射方位相位和俯仰相位 接收天线 发射天线顺序
arr = virtual_array.virtual_arr;
% 初阶版本顺序发射可行
% arr = cfgOut.array; % 虚拟阵元排列
arrNum = size(arr, 1); % 虚拟阵列数目
arrdelay = zeros(1, arrNum) + Tc * reshape(repmat([0:numTx-1], numRx,1), 1, arrNum); % [0 0 0 0 1 1 1 1 2 2 2 2]
% arrdelay = zeros(1, arrNum); % 这一项相当于常数项 本来考虑的是天线间时延
delayTx = Tc * numTx; % 发射天线维度积累的时延
arrdx = cfgOut.arrdx; % 归一化阵元间距
arrdy = cfgOut.arrdy;
c = physconst('LightSpeed'); % 光速
lambda = c / fc; % 波长
Ts = 1 / fs; % ADC采样时间
% 信号模型
signal = zeros(ADCNum, TotalChirpNum, arrNum); % 信号
noiseI = normrnd(0, 1, ADCNum, TotalChirpNum, arrNum); % 正态噪声
noiseQ = normrnd(0, 1, ADCNum, TotalChirpNum, arrNum); % 正态噪声
% noiseI = exprnd(1, ADCNum, TotalChirpNum, arrNum); % 指数分布噪声
% noiseQ = exprnd(1, ADCNum, TotalChirpNum, arrNum); % 指数分布噪声
for tarId = 1 : tarNum
disp(strcat(['正在获取目标',num2str(tarId),'的数据']));
% 目标参数
targetR = tarOut.range(tarId); % 距离 m
targetV = tarOut.velocity(tarId); % 速度 m/s
targetAzi = tarOut.Azi(tarId); % 目标方位角 度
targetEle = tarOut.Ele(tarId); % 目标俯仰角 度
targetRCS = tarOut.RCS(tarId); % 目标RCS值
targetGt = tarOut.Gt(tarId); % 发射增益
targetGr = tarOut.Gr(tarId); % 接收增益
[tarSNR] = CalculateSNR(targetR, targetRCS, targetGt, targetGr, lambda, Pt, Fn, Ls, validB); % 信噪比
A = sqrt(2 * db2pow(tarSNR)); % 信号幅度
targPhi0 = rand * 2 * pi; % 产生随机相位[0 2*pi] 这里的目标初始相位也可以根据上面定义的形式来写 但影响不大
tempSigI = zeros(ADCNum, TotalChirpNum, arrNum); % I路缓存信号
tempSigQ = zeros(ADCNum, TotalChirpNum, arrNum); % Q路缓存信号
tempSig = zeros(ADCNum, TotalChirpNum, arrNum); % I路缓存信号
for channelId = 1 : arrNum
for chirpId = 1 : TotalChirpNum
for adcId = 1 : ADCNum
% 动目标时延
tarDelay = 2 * targetR / c + 2 * targetV * ((chirpId - 1) * delayTx + arrdelay(channelId) + adcId * Ts) / c;
% 控制来波方向
% 初阶版本的目标相位表达
% tarPhi = targPhi0 + 2 * pi * (arr(1, channelId) * sind(targetAzi) * arrdx + ...
% arr(2, channelId) * sind(targetEle) * arrdy) + deg2rad(antennaPhase(channelId)); % 考虑阵元初始、随机相位和目标波达角
% tarPhi0 = 2 * targetR * fc / c; % 目标初始相位 也可以用随机相位
tarPhi = targPhi0 + 2 * pi * (arr(channelId) * sind(targetAzi) * arrdx + ...
arr(arrNum+channelId) * sind(targetEle) * arrdy) + deg2rad(antennaPhase(channelId)); % 这里考虑了不同发射顺序天线相位
% 中频信号：exp[1j * 2 *pi * fc * tau + 2 * pi * S * t * tau - pi * tau * tau] t的范围为0-Tc tau的范围为t+nTc
tempSigI(adcId, chirpId, channelId) = A * cos(2 * pi * (fc * tarDelay + Slope * tarDelay * adcId * Ts - Slope * tarDelay * tarDelay / 2) + tarPhi);
tempSigQ(adcId, chirpId, channelId) = A * sin(2 * pi * (fc * tarDelay + Slope * tarDelay * adcId * Ts - Slope * tarDelay * tarDelay / 2) + tarPhi);
tempSig(adcId, chirpId, channelId) = tempSigI(adcId, chirpId, channelId) + 1j * tempSigQ(adcId, chirpId, channelId);
end
end
end
signal = signal + tempSig; % 多目标信号相加
end
signal = signal - (noiseI + noiseQ); % 考虑噪声
% signal = reshape(signal, size(signal,1), size(signal,2), numRx, numTx);
end
信号处理栈
在上一章节，我们已经得到了radar cube数据，这个数据和毫米波雷达ADC数据是异曲同工的。那么，针对现有数据可以展开信号处理栈的分析，简言之，传统的信号处理栈为RDA，如本文开篇的图所示。
首先，沿着radar cube数据的快时间和慢时间维度分别使用FFT（快速傅立叶变换）的到距离多普勒域数据RDD，这可以表示为如下形式：
\[x_n(\lambda,\zeta)=\sum_{k_{tgt}=1}^{K_\text{tgt}}a_{k_{tgt}}{k_{tgt}}W_{\text{ft}}(\lambda-\lambda_{k_{tgt}})W_{\text{st}}(\zeta-\zeta_{k_{tgt}})e^{j\varphi_{k_{tgt}}}\cdot e^{j\frac{2\pi}{\lambda}(\mathbf{u}_{k\mathrm{tgt}}\cdot\mathbf{p}_n)}+\Omega_n(\lambda,\zeta)
\]
其中\(W_{ft}(\lambda)\)和\(W_{st}(\zeta)\)分别表示在快时间和慢时间维度上应用的窗函数傅里叶变换，这实现了信号移位至目标对应的距离维和多普勒峰值频率处：\(\lambda_{k_{tgt}}\)和\(\zeta_{k_{tgt}}\)。目标初始相位项\(e^{j\varphi_{k_{tgt}}}\)和相位分布项\(e^{j\frac{2\pi}{\lambda}(\mathbf{u}_{k\mathrm{tgt}}\cdot\mathbf{p}_n)}\)保持不变，这是因为沿距离-多普勒维度的傅里叶变换并不影响天线维度的相位信息。\(\Omega_n(\lambda,\zeta)\)则表示频域噪声。
为了进一步对距离-多普勒变换后的信号处理，将上述混合信号表示为如下形式：
设距离多普勒域的radar cube的元素为\(x_{l,m,n}~~=~~x[l,m,n]\)，快时间和慢时间维度的FFT点数对应\({P_{\mathrm{FFT}}}\)和\({Q_{\mathrm{FFT}}}\)
\[\lambda=\frac{2\pi}{P_{\mathrm{FFT}}}l,\quad\zeta=\frac{2\pi}{Q_{\mathrm{FFT}}}m,\quad x[l,m,n]=x_n\left(\frac{2\pi}{P_{\mathrm{FFT}}}l,\frac{2\pi}{Q_{\mathrm{FFT}}}m\right).
\]
其次，我们需要对距离-多普勒域的radar cube进行目标检测，这通常用恒虚警率检测算法(cfar, constant false alarm rate)来实现，在正式进行检测前，我们需要对radar cube沿天线维度展开，通常可以选择一个天线通道对应的距离多普勒数据或使用所有天线通道的距离多普勒数据进行非相干积累处理（波束成形）得到cfar检测输入数据，这个过程可以表示如下：(下面第一式为方法1，第二式为方法2，\(N_{det}\)为单天线通道，\(N^\prime\)为全部天线通道)
\[X^\prime=|x'[l,m]|^2=|x[l,m,n=N_{\det}]|^2
\]
\[X^\prime=|x'[l,m]|^2=\sum\limits_{n\in N'}|x[l,m,n]|^2
\]
那么，现在我们可以对处理后的距离多普勒矩阵\(X^\prime\)进行cfar检测。通常，cfar检测针对数据维度可以分为1D-cfar和2D-cfar，针对处理思想可以分为均值类、有序统计类、自动筛选类和自适应类（这部分详细细节可以参考何友院士和关键老师的雷达目标检测与恒虚警处理书，后续会整理我的本科毕设中关于这方面的内容）。实际上，我们可以考虑两次1D-CFAR即对距离多普勒矩阵\(X^\prime\)的距离和多普勒维度分别作cfar并关联匹配目标，或者考虑1次2d-cfar即直接在二维距离多普勒矩阵\(X^\prime\)上作二维检测，这里我们主要用2d-cfar来检测目标。
在下面我们通过2D-CFAR示意图来更好地理解原理。浅蓝色单元构成的大矩形区域为输入的距离-多普勒矩阵\(X^\prime\)深红色单元为检测单元，浅红色单元为保护单元（防止旁瓣能量信号泄露主瓣形成掩蔽效应），浅绿色为参考单元，由深红色单元、浅红色单元和浅绿色单元构成的小矩形区域作为滑动窗窗；依次移动滑动窗在大矩形区域内遍历检测目标，移动过程中需要保证检测单元的中心从黄色的起始中心单元开始依次完成自左向右、自上而下的二维扫描，直至完成整个大矩形区域的遍历。【注意到，我们在这里强调了"从黄色的起始中心单元"，这部分需要在程序中重点处理，否则可能会导致部分区域是没有参与检测的。】
以单元平均恒虚警检测算法(ca-cfar，cell-averaging cfar)为例，首先需要根据预设虚警概率和滑动窗尺寸等先验参数计算门限因子\(\alpha\)；在每次扫描中，需求求解参考窗内所有检测单元对应的信号电平的和并取平均值，这一估计值可以表示为\(\bar{X}\)，将估计阈值\(\bar{X}\cdot \alpha\)和检测单元下的信号电平进行比较，如果检测单元下信号电平大于估计阈值\(\bar{X}\cdot \alpha\)，则认为检测单元出为目标信号，否则认为噪声信号。
在这里我们通过以下程序来更好地帮助理解2d-cfar过程：
function [cfarOut] = CFAR_2D(RDM, Pfa, TestCells, GuardCells)
% rd_map：输入数据
% Pfa：虚警概率
% TestCells：测试单元 [4 4]
% GuardCells：保护单元 [4 4]
% By Xuliang
% CFAR检测器
detector = phased.CFARDetector2D('TrainingBandSize',TestCells, ...
'ThresholdFactor','Auto','GuardBandSize',GuardCells, ...
'ProbabilityFalseAlarm',Pfa,'Method','SOCA','ThresholdOutputPort',true, 'NoisePowerOutputPort',true);
N_x = size(RDM,1); % 距离单元
N_y = size(RDM,2); % 多普勒单元
% 获取二维滑动窗口的size
Ngc = detector.GuardBandSize(2); % 保护窗列向
Ngr = detector.GuardBandSize(1); % 保护窗行向
Ntc = detector.TrainingBandSize(2); % 参考窗行向
Ntr = detector.TrainingBandSize(1); % 参考窗行向
cutidx = [];
colstart = Ntc + Ngc + 1; % 列窗首
colend = N_y + ( Ntc + Ngc); % 列窗尾
rowstart = Ntr + Ngr + 1; % 行窗首
rowend = N_x + ( Ntr + Ngr); % 行窗尾
for m = colstart:colend
for n = rowstart:rowend
cutidx = [cutidx,[n;m]]; % 完整二维窗
end
end
ncutcells = size(cutidx,2); % 获取窗口遍历区间
rd_map_padding = repmat(RDM, 3, 3); % 对RDM补零
chosen_rd_map = rd_map_padding(N_x+1-Ntr-Ngr:2*N_x+Ntr+Ngr,N_y+1-Ntc-Ngc:2*N_y+Ntc+Ngc);
[dets, ~, noise] = detector(chosen_rd_map,cutidx); % 输出CFAR结果
cfar_out = zeros(size(chosen_rd_map)); % 检测点输出
noise_out = zeros(size(chosen_rd_map)); % 噪声
snr_out = zeros(size(chosen_rd_map)); % 信噪比
for k = 1:size(dets,1)
if dets(k) == 1
cfar_out(cutidx(1,k),cutidx(2,k)) = dets(k);
noise_out(cutidx(1,k),cutidx(2,k)) = noise(k);
snr_out(cutidx(1,k),cutidx(2,k)) = chosen_rd_map(cutidx(1,k),cutidx(2,k));
end
end
cfarOut = {};
cfarOut.cfarMap = cfar_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y);
cfarOut.snrOut = snr_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y);
cfarOut.noiseOut = noise_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y);
cfarOut.snrOut = cfarOut.snrOut ./ (eps + cfarOut.noiseOut);
end
经过cfar检测后，我们可以得到一组由\(K_{det}\)个检测目标构成的集合：
\[\mathrm{cfar}\left(\left|x'[l,m]\right|^2\right)=\{(r_1,v_{\mathrm{r},1}),(r_2,v_{\mathrm{r},2}),\ldots,(r_{k_{\det}},v_{\mathrm r,k_{\det}}),\ldots,(r_{K_{\det}},v_{\mathrm r,K_{\det}})\}
\]
其中，每个检测目标由与其关联的参数元组进行表示，每个参数元组\((r_{K_{\det}},v_{\mathrm r,K_{\det}})\)中包含了目标的距离和速度，这可以通过距离或多普勒维度的峰值索引通过物理映射实现距离和速度解析，这里需要用到距离分辨率和速度分辨率的概念（可以阅读先前的博文进行回顾）。
下面，我们可以进入doa估计；不过我们需要准备好用于空间谱估计的输入信号，这需要用到由上面缓存的radar cube数据和经过cfar检测得到的目标距离-速度集合\(\mathrm{cfar}\left(\left|x'[l,m]\right|^2\right)\)来提取空域信号\(\mathbf{x}_{k\mathrm{det}}\)，表示如下：（\(N\)表示阵元数目）
\[\mathbf{x}_{k\mathrm{det}}=x[l_{k\mathrm{det}},m_{k\mathrm{det}},n]=[x_0,\ldots,x_{N-1}]^T
\]
基于经过距离多普勒FFT的混合信号模型\(x_n(\lambda,\zeta)\)，我们可以建立空域信号模型为：
\[\mathbf{x}_{k_{\det}}=\sum_{k_{\text{doa}}=1}^{K_{\text{doa}}}s_{k_{\text{doa}}}\mathbf{a}(\phi_{k_{\text{doa}}},\epsilon_{k_{\text{doa}}})+\boldsymbol{\omega}
\]
不难发现上述信号表示是标准的空间谱估计信号模型X=AS+N，其中A为流形，S为信号空间，N为噪声空间（上面用\(\omega\)表示噪声项，通常假设满足复高斯白噪声分布），X为观测空间；在这里，我们也可以将A视为完备基或字典，字典原子或称导向矢量\(\mathbf{a}(\phi_{k_{\text{doa}}},\epsilon_{k_{\text{doa}}})\)满足线性无关，那么目标的来波信号可以由字典中的原子通过复权重\(s_{k_{doa}}\)加权得到，导向矢量和复权重的表示分布如下：
\[\mathbf{a}(\phi_{k_\text{doa}},\epsilon_{k_\text{doa}})=e^{j(2\pi/\lambda)\mathbf{u}_{k_\text{doa}}\mathbf{p}_n}
\]
\[s_{k\text{doa}}=a_{k\text{doa}}W_{\text{ft}}(\lambda-\lambda_{k\text{doa}})W_{\text{st}}(\zeta-\zeta_{k\text{doa}})e^{j\varphi_{k}\text{doa}}
\]
doa估计的本意是在给定峰值目标空域信号\(\mathbf{x}_{k\mathrm{det}}\)和阵列集合结构\(p_n\)的基础上估计出目标数目\(K_{doa}\)、目标来波方向\(\phi_{k_{doa}}\)和\(\epsilon_{k_{doa}}\)，即如下目标：\((r_{k_\mathrm{det}},v_{\mathrm{r},k_\mathrm{det}},\mathbf{x}_{k_\mathrm{det}})\to(r_{k_\mathrm{tgt}},v_{r,k_\mathrm{tgt}},\phi_{k_\mathrm{tgt}},\epsilon_{k_\mathrm{tgt}})\)。在本处，我们可以对上述阵列信号完成对应维度的天线映射后，首先对方位天线维度做FFT处理，并认为方位维空间最大谱峰所在处为对应的目标方位角；其次，基于目标方位天线的峰值索引，我们沿其俯仰维度做FFT处理，并认为俯仰维空间最大谱峰所在处为对应的目标俯仰角，以此实现目标信号方位和俯仰维度的匹配关联。这种方法可以理解为“两步走”战略，和上面提及的两次1D-CFAR是类似的。
当然，我们也可以使用波束成形CBF、线性滤波Capon、子空间方法（MUSIC，ESPRIT，DML等）和压缩感知方法(OMP，IAA和LISTA等)，这些内容我们将在后面的章节中讨论，在这里需要注意一点，我们要保证方位天线维度上的空间谱估计不能够破坏俯仰天线维度观测信号的相位信息，否则会导致俯仰角的估计错误。
至此，我们已经基本完成了RDA信号处理栈的基本物理量提取并最终得到点云表示\((r_{k_\mathrm{tgt}},v_{r,k_\mathrm{tgt}},\phi_{k_\mathrm{tgt}},\epsilon_{k_\mathrm{tgt}})\)。所述过程仍存在一些细节性内容没有解释，这部分可以结合代码更好地去理解，下面给出核心信号处理实现部分。
for frame_id = 1 : Frame % Frame
adcData = squeeze(RawData(:, :, frame_id, :));
%% 距离维FFT和多普勒FFT
disp(strcat(['=====','RD-MAP生成','=====']));
tic
fftOut = rdFFT(adcData, IQFlag);
toc
rangeFFTOut = fftOut.rangeFFT;
dopplerFFTOut = fftOut.dopplerFFT;
if ShowRange
figure(2);
set(gcf,'unit','centimeters','position',[10,0,10,10])
plot(rangeIndex, db(rangeFFTOut(:,1,1)));
xlabel('Range(m)');ylabel('Amplitude(dB)');
title(strcat(['第',num2str(frame_id),'帧-目标距离分布']));grid minor;
pause(0.1);
end
if ShowRD
figure(3);
set(gcf,'unit','centimeters','position',[20,12,10,10])
imagesc(rangeIndex, velocityIndex, db(dopplerFFTOut(:,:,1).'));
xlabel('Range(m)');ylabel('Velocity(m/s)'); colormap('jet');
title(strcat(['第',num2str(frame_id),'帧目标-距离多普勒分布']));
grid minor; axis xy;
pause(0.1);
end
%% 非相干积累
disp(strcat(['=====','非相干积累','=====']));
RDM = dopplerFFTOut;
tic
[accumulateRD] = incoherent_accumulation(RDM);
toc
%% CFAR检测器
disp(strcat(['=====','恒虚警检测','=====']));
Pfa = 1e-3; % 虚警概率
TestCells = [8, 8]; % 参考窗
GuardCells = [2, 2]; % 保护窗
tic
[cfarOut] = CFAR_2D(accumulateRD, Pfa, TestCells, GuardCells);
toc
cfarMap = cfarOut.cfarMap; % 检测点输出
noiseOut = cfarOut.noiseOut; % 噪声输出
snrOut = cfarOut.snrOut; % 信噪比输出
if ShowCFAR
figure(4);
set(gcf,'unit','centimeters','position',[20,0,10,10])
imagesc(rangeIndex, velocityIndex, cfarMap.');
xlabel('Range(m)');ylabel('Velocity(m/s)'); colormap('jet');
title(strcat(['第',num2str(frame_id),'帧目标-CFAR检测结果']));
grid minor;axis xy;
pause(0.1);
end
%% 峰值聚合-获取点目标
disp(strcat(['=====','峰值聚合','=====']));
[range_idx, doppler_idx] = find(cfarMap);
cfar_out_idx = [range_idx doppler_idx]; % 获取CFAR输出结果的行列索引
tic
[rd_peak_list, rd_peak] = peakFocus(db(accumulateRD), cfar_out_idx);
toc
peakMap = zeros(size(cfarMap)); % 峰值聚合结果矩阵
for peak_idx = 1 :size(rd_peak_list, 2)
peakMap(rd_peak_list(1,peak_idx), rd_peak_list(2,peak_idx)) = 1;
end
if ShowPeak
figure(5);
set(gcf,'unit','centimeters','position',[30,12,10,10])
imagesc(rangeIndex, velocityIndex, peakMap.');
xlabel('Range(m)');ylabel('Velocity(m/s)'); colormap('jet');
title(strcat(['第',num2str(frame_id),'帧目标-峰值聚合结果']));
grid minor;axis xy;
pause(0.1);
end
%% DOA/AOA估计
disp(strcat(['=====','DOA/AOA估计','=====']));
cfgDOA.FFTNum = 180; % FFT点数
% 这里虽然封装了不同的DOA算法 但是需要注意的是 可用的算法有限 在本套代码里 建议使用的是FFT-FFT和FFT-MUSIC等
% 因为MUSIC-MUSIC的使用会导致在方位维度空间谱估计时破坏相位信号 有损俯仰维度的相位估计
cfgDOA.AziMethod = 'FFT'; % 方位维度DOA估计方法
cfgDOA.EleMethod = 'MUSIC'; % 俯仰维度DOA估计方法
cfgDOA.thetaGrids = linspace(-90, 90, cfgDOA.FFTNum); % 空间网格
cfgDOA.AzisigNum = 1; % 约束每个RD-CELL上的信源数目
cfgDOA.ElesigNum = 1; % 约束每个方位谱峰上的信源数目
aziNum = length(virtual_array.noredundant_aziarr); % 方位天线数目
% Aset = exp(-1j * 2 * pi * arrDx * [0:aziNum]' * sind(cfgDOA.thetaGrids)); % 稀疏字典设计
targetPerFrame = {};
targetPerFrame.rangeSet = [];
targetPerFrame.velocitySet = [];
targetPerFrame.snrSet = [];
targetPerFrame.azimuthSet = [];
targetPerFrame.elevationSet = [];
if ~isempty(rd_peak_list) % 非空表示检测到目标
rangeVal = (rd_peak_list(1, :) - 1) * range_res; % 目标距离
speedVal = (rd_peak_list(2, :) - ChirpNum / 2 - 1) * doppler_res; % 目标速度
doaInput = zeros(size(rd_peak_list, 2), arrNum);
for tar_idx = 1 :size(rd_peak_list, 2)
doaInput(tar_idx, :) = squeeze(dopplerFFTOut(rd_peak_list(1, tar_idx), rd_peak_list(2, tar_idx), :)); % tarNum * arrNum
end
doaInput = reshape(doaInput, [], numRx, numTx);
% 方位角估计前需要考虑多普勒补偿
[com_dopplerFFTOut] = compensate_doppler(doaInput, cfgOut, rd_peak_list(2, :), speedVal, rangeVal);
tic
for peak_idx = 1 : size(rd_peak_list, 2) % 遍历检测到的每个目标
snrVal = mag2db(snrOut(rd_peak_list(1, peak_idx), rd_peak_list(2, peak_idx))); % 信噪比的提升是由于chirpNum*ADCNum的积累
tarData = squeeze(com_dopplerFFTOut(peak_idx, :,:));
% aziarr = unique(arr(1,arr(2,:)==0)); % 初阶版本获取方位维度天线排列
% aziArrData = arrData(aziarr+1); % 获取方位维度信号
% 方位角解析
sig = tarData;
sig_space = zeros(max(virtual_array.azi_arr)+1,max(virtual_array.ele_arr)+1); % 初始化信号子空间
for trx_id = 1 : size(cfgOut.sigIdx,2)
sig_space(cfgOut.sigSpaceIdx(1, trx_id), cfgOut.sigSpaceIdx(2,trx_id)) = sig(cfgOut.sigIdx(1,trx_id), cfgOut.sigIdx(2,trx_id)); % 重排后的信号空间
end
% 至此我们生成的信号子空间维度为 方位虚拟天线数目 * 俯仰虚拟天线数目
eleArrData = zeros(cfgDOA.FFTNum, size(sig_space,2)); % 俯仰维度数据
for ele_idx = 1 : size(sig_space, 2) % 这里采取遍历是为了适配不同的空间谱估计方法
tmpAziData = sig_space(:, ele_idx);
[azidoaOut] = azimuthDOA(tmpAziData, cfgDOA); % 提取第一列方位维度天线信息进行方位角估计
eleArrData(:, ele_idx) = azidoaOut.spectrum(:); % 复空间谱
end
for azi_peak_idx = 1 : length(azidoaOut.angleVal) % 对方位维度检测的谱峰进行检索
tmpEleData = eleArrData(azidoaOut.angleIdx(azi_peak_idx), :).'; % 获取与方位维目标关联的信号
[eledoaOut] = elevationDOA(tmpEleData, cfgDOA); % 进行俯仰角估计
% 关联目标的距离、多普勒、信噪比、方位和俯仰信息
aziVal = azidoaOut.angleVal;
eleVal = eledoaOut.angleVal;
targetPerFrame.rangeSet = [targetPerFrame.rangeSet, rangeVal(peak_idx)];
targetPerFrame.velocitySet = [targetPerFrame.velocitySet, speedVal(peak_idx)];
targetPerFrame.snrSet = [targetPerFrame.snrSet, snrVal];
targetPerFrame.azimuthSet = [targetPerFrame.azimuthSet,aziVal];
targetPerFrame.elevationSet = [targetPerFrame.elevationSet,eleVal];
end
end
toc
%% 点云生成
disp(strcat(['=====','点云生成','=====']));
tic
pcd_x = targetPerFrame.rangeSet .* cosd(targetPerFrame.elevationSet) .* sind(targetPerFrame.azimuthSet);
pcd_y = targetPerFrame.rangeSet .* cosd(targetPerFrame.elevationSet) .* cosd(targetPerFrame.azimuthSet);
pcd_z = targetPerFrame.rangeSet .* sind(targetPerFrame.elevationSet);
PointSet = [pcd_x.', pcd_y.', pcd_z.', targetPerFrame.velocitySet.', targetPerFrame.snrSet.'];
toc
%% 点云聚类
eps = 1.1; % 邻域半径
minPointsInCluster = 3; % 簇内最小点数阈值
xFactor = 1; % 变大控制距离变大，变小分类距离变小 椭圆
yFactor = 1; % 变大控制角度变大，变小分类距离变小 椭圆
figure(6);
set(gcf,'unit','centimeters','position',[30,0,10,10])
disp(strcat(['=====','点云聚类','=====']));
tic
[sumClu] = dbscanClustering(eps,PointSet,xFactor,yFactor,minPointsInCluster,frame_id); % DBSCAN聚类
toc
end
end
在本文的尾巴，本文将展示本程序所实现的可视化功能。在以下仿真中，目标距离为[5 10 15 18]，目标速度为[0.4 -0.3 0 0.6]，目标方位角为[15 -2 30 -25]，目标俯仰角为[0 15 -5 8]，目标RCS为[20 20 10 20]，仿真实验平台为TI AWR1243BOOST，其他平台均可实现类似效果。
TI级联雷达例程走读
import time
import numpy as np
import cupy as cp
import scipy.io as scio
class NormalModeProcess(object):
def __init__(self, infos, data_real, data_imag, calib, apply_vmax_extend=False):
self.infos = infos
# data_real.shape[2]=data_imag.shape[2]=16: ordered by dev[0,1,2,3]_rx[0,1,2,3],
# data_real.shape[3]=data_imag.shape[3]=12: ordered by self.tx_id_transfer_order
self.data_raw = data_real + data_imag * 1j # 恢复实部和复部数据 256*128*16*12
self.rangeFFT_size = self.data_raw.shape[0] # 距离FFT点数
self.dopplerFFT_size = self.data_raw.shape[1] # 多普勒FFT点数
self.azimuthFFT_size = 256 # 方位FFT点数
self.elevationFFT_size = 256 # 俯仰FFT点数
self.TI_Cascade_Antenna_DesignFreq_GHz = 76.8
# doa_unitDis = 0.5 * (self.infos['startFreqConst_GHz'] + 256 / self.infos['digOutSampleRate_ksps']*self.infos['freqSlopeConst_MHz_usec']/2) / TI_Cascade_Antenna_DesignFreq_GHz
self.doa_unitDis = 0.5 * self.infos['freq_center_GHz'] / self.TI_Cascade_Antenna_DesignFreq_GHz # 这里应该是0.5 表示的是归一化阵元间距 半波长
self.range_bins = np.arange(self.rangeFFT_size) * self.infos['range_res_m'] # 距离
self.doppler_bins = np.arange(-self.dopplerFFT_size / 2, self.dopplerFFT_size / 2) * self.infos[
'velocity_res_m_sec'] # 多普勒速度
self.azimuth_bins = np.arange(-self.azimuthFFT_size, self.azimuthFFT_size, 2) * np.pi / self.azimuthFFT_size # 方位角度
self.elevation_bins = np.arange(-self.elevationFFT_size, self.elevationFFT_size,
2) * np.pi / self.elevationFFT_size # 俯仰角度
self.calib = calib # 校准矩阵的读取
self.num_rx = self.infos['numRXPerDevice'] * self.infos['numDevices'] # 发射天线数目
self.num_tx = self.infos['numTXPerDevice'] * self.infos['numDevices'] # 接收天线数目
self.rx_id = np.arange(self.num_rx) # 发射天线顺序
# 接收和发射天线的阵元排列分布
self.rx_id_onboard = np.array([13, 14, 15, 16, 1, 2, 3, 4, 9, 10, 11, 12, 5, 6, 7, 8]) - 1
self.rx_position_azimuth = np.array([11, 12, 13, 14, 50, 51, 52, 53, 46, 47, 48, 49, 0, 1, 2, 3])
self.rx_position_elevation = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
self.tx_id_transfer_order = np.array([12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]) - 1
self.tx_id = np.arange(self.num_tx)
self.tx_id_onboard = np.array([12, 11, 10, 3, 2, 1, 9, 8, 7, 6, 5, 4]) - 1
self.tx_position_azimuth = np.array([11, 10, 9, 32, 28, 24, 20, 16, 12, 8, 4, 0])
self.tx_position_elevation = np.array([6, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
self.apply_vmax_extend = apply_vmax_extend # 速度扩展？
if self.apply_vmax_extend:
self.min_dis_apply_vmax_extend = 10 # 最小速度扩展值？
self.pcd = None
self.heatmapBEV = None
self.heatmap4D = None
self.scale_factor = np.array([0.2500, 0.1250, 0.0625, 0.0312, 0.0156, 0.0078, 0.0039, 0.0020]) # 1/2 * [1/2,1/4,...]
def __timer__(self, fn):
start_time = time.time()
fn()
end_time = time.time()
print('{}() function consume: {:.3f} s'.format(fn.__name__, end_time - start_time))
def run(self, generate_pcd, generate_heatmapBEV, generate_heatmap4D):
self.generate_pcd = generate_pcd
self.generate_heatmapBEV = generate_heatmapBEV
self.generate_heatmap4D = generate_heatmap4D
if self.generate_pcd or self.generate_heatmapBEV or self.generate_heatmap4D:
self.__timer__(self.__calibrate_rawdata__)
self.__timer__(self.__re_order__)
self.__timer__(self.__calculate_antenna_array__)
self.__timer__(self.__rangeFFT__)
self.__timer__(self.__dopplerFFT__)
if self.generate_pcd:
self.__timer__(self.__rdm_cfar__)
if self.res_rdm_cfar is not None:
self.__timer__(self.__doa__)
if self.res_doa is not None:
self.__timer__(self.__generate_pcd__)
if self.generate_heatmapBEV:
self.__timer__(self.__generate_heatmapBEV__)
if self.generate_heatmap4D:
self.__timer__(self.__generate_heatmap4D__)
def get_range_bins(self):
return self.range_bins
def get_doppler_bins(self):
return self.doppler_bins
def get_azimuth_bins(self, unit_degree=True):
if unit_degree:
return np.arcsin(self.azimuth_bins / 2 / np.pi / self.doa_unitDis) / np.pi * 180
else:
return np.arcsin(self.azimuth_bins / 2 / np.pi / self.doa_unitDis)
def get_elevation_bins(self, unit_degree=True):
if unit_degree:
return np.arcsin(self.elevation_bins / 2 / np.pi / self.doa_unitDis) / np.pi * 180
else:
return np.arcsin(self.elevation_bins / 2 / np.pi / self.doa_unitDis)
def get_pcd(self):
return self.pcd
def get_heatmapBEV(self):
return self.heatmapBEV
def get_heatmap4D(self):
return self.heatmap4D
def __calibrate_rawdata__(self):
adc_sample_rate = self.infos['digOutSampleRate_ksps'] * 1e3 # adc采样率
chirp_slope = self.infos['freqSlopeConst_MHz_usec'] * 1e6 / 1e-6 # chirp斜率
num_sample = self.infos['numAdcSamples'] # ADC数目 256
num_loop = self.infos['numLoops']
range_mat = self.calib['calibResult']['RangeMat'] # 12 * 16
fs_calib = self.calib['params']['Sampling_Rate_sps']
slope_calib = self.calib['params']['Slope_MHzperus'] * 1e6 / 1e-6
calibration_interp = 5
peak_val_mat = self.calib['calibResult']['PeakValMat']
phase_calib_only = 1
tx_id_ref = self.tx_id_transfer_order[0] # 以某个通道为基准进行校准
# 构建频率补偿矩阵
freq_calib = 2 * np.pi * (
(range_mat[self.tx_id_transfer_order, :] - range_mat[tx_id_ref, 0])
* fs_calib / adc_sample_rate * chirp_slope / slope_calib # 12 * 16
) / (num_sample * calibration_interp) # calibration_interp对应前面求校准矩阵时距离FFT点数的放大倍数 实际adc是256 # 12 * 16
correction_vec = np.exp(
1j * np.arange(num_sample).reshape((-1, 1, 1)) * np.expand_dims(freq_calib.T, axis=0)
).conj() # 256 * 1 * 1, 1 * 16 * 12 ---> 256 * 1 * 16 * 12
freq_correction_mat = np.expand_dims(correction_vec, axis=1) # 256 * 1 * 16 * 12
# 构建相位校准矩阵
phase_calib = peak_val_mat[tx_id_ref, 0] / peak_val_mat[self.tx_id_transfer_order, :] # 参考通道的值/其余通道的值 12 * 16
# remove amplitude calibration
if phase_calib_only == 1:
phase_calib = phase_calib / np.abs(phase_calib) # 只校准相位
phase_correction_mat = np.expand_dims(np.expand_dims(phase_calib.T, axis=0), axis=0) # 1 * 1 * 16 * 12
self.data_calib = self.data_raw * freq_correction_mat * phase_correction_mat # 256 * 128 * 16 * 12 完成校准
def __re_order__(self): # 重排后的数据 256 * 128 * 16 * 12
# data_reordered.shape[2]=16: ordered by self.rx_id_onboard
# data_reordered.shape[3]=12: ordered by self.tx_id_transfer_order
self.data_reordered = self.data_calib[:, :, self.rx_id_onboard, :] # 发射天线维度进行重排
def __calculate_antenna_array__(self): # 计算天线孔径
# 收发天线在方位向形成的相位 tile进行拷贝扩展
self.virtual_array_azimuth = np.tile(self.tx_position_azimuth[self.tx_id_transfer_order], (self.num_rx, 1)) + \
np.tile(self.rx_position_azimuth[self.rx_id_onboard], (self.num_tx, 1)).T # 16*12
# 收发天线在俯仰向形成的相位
self.virtual_array_elevation = np.tile(self.tx_position_elevation[self.tx_id_transfer_order],
(self.num_rx, 1)) + \
np.tile(self.rx_position_elevation[self.rx_id_onboard], (self.num_tx, 1)).T # 16*12
self.virtual_array_tx_id = np.tile(self.tx_id_transfer_order.reshape(1, -1), (self.num_rx, 1)) # 12 * 16
self.virtual_array_rx_id = np.tile(self.rx_id_onboard.reshape(-1, 1), (1, self.num_tx)) # 16 * 12
# azimuth, elevation, rx_id, tx_id
self.virtual_array = np.hstack(
(
self.virtual_array_azimuth.reshape((-1, 1), order='F'), # 192*1 优先按列走
self.virtual_array_elevation.reshape((-1, 1), order='F'),
self.virtual_array_rx_id.reshape((-1, 1), order='F'),
self.virtual_array_tx_id.reshape((-1, 1), order='F')
)
)
# get antenna_noredundant
_, self.virtual_array_index_noredundant = np.unique(
self.virtual_array[:, :2], axis=0, return_index=True # 返回下标值
) # 134个冗余，返回的是冗余阵元的索引
self.virtual_array_noredundant = self.virtual_array[self.virtual_array_index_noredundant, :] # 获取去冗余后的方位和俯仰的虚拟阵元
# 方位 俯仰 接收天线序号 发射天线序号
# get antenna_redundant
self.virtual_array_index_redundant = np.setxor1d(
np.arange(self.virtual_array.shape[0]), # 192*1
self.virtual_array_index_noredundant # 134*1
) # 找到两个数组的异或集合 即未用的58个虚拟阵元【冗余虚拟阵元】
self.virtual_array_redundant = self.virtual_array[self.virtual_array_index_redundant, :] # 58*4 返回未用阵元索引
# 查找并关联重叠的Rx_Tx对
virtual_array_info_overlaped_associate = []
for i in range(self.virtual_array_index_redundant.shape[0]): # 遍历每个未用阵元 58
mask = (self.virtual_array_noredundant == self.virtual_array_redundant[i])
# 这里是对冗余阵元和非冗余阵元进行匹配 输出结果为[[T,F,F,T],...]等形式 维度为134*4
mask = np.logical_and(mask[:, 0], mask[:, 1]) # 比较发射天线和接收天线是否有一致的 如果存在一致则掩码为true 134*1
info_associate = self.virtual_array_noredundant[mask][0] # 1*4 找到冗余和非冗余中接收和发射天线重叠的元素
info_overlaped = self.virtual_array_redundant[i] # 当前冗余阵元值 1*4
virtual_array_info_overlaped_associate.append(
np.concatenate((info_associate, info_overlaped)).tolist() # 将有关联的非冗余阵元和冗余阵元加入列表 重叠的意思指的是由不同TX-RX天线对形成的相位差是相同的
)
# azimuth, elevation, rx_associated, tx_associated, azimuth, elevation, rx_overlaped, tx_overlaped
virtual_array_info_overlaped_associate = np.array(virtual_array_info_overlaped_associate) # 58 * 8, 前4个为匹配的非冗余阵元（134） 后四个为冗余阵元（58）
diff_tx = abs(virtual_array_info_overlaped_associate[:, 7] - virtual_array_info_overlaped_associate[:, 3]) # 计算接收天线的位置差
virtual_array_info_overlaped_diff1tx = virtual_array_info_overlaped_associate[diff_tx == 1] # 筛选出发射天线位置差为1的重叠阵元
sorted_index = np.argsort(virtual_array_info_overlaped_diff1tx[:, 0]) # 按升序排序 返回排序后的元素的索引值
virtual_array_info_overlaped_diff1tx = virtual_array_info_overlaped_diff1tx[sorted_index, :]
self.virtual_array_info_overlaped_diff1tx = virtual_array_info_overlaped_diff1tx # 这有啥意义？
# 找到方位维度的虚拟阵元排列
self.virtual_array_noredundant_row1 = self.virtual_array_noredundant[self.virtual_array_noredundant[:, 1] == 0]
# 86 * 4, 方位维度非冗余阵元在俯仰向相位差为0 方位 俯仰 发射位置 接收位置
def __rangeFFT__(self, window_enable=True, scale_on=False):
self.data_rangeFFT = self.data_reordered # 获取重排后的数据 256 * 128 * 16 * 12
# 直流滤波 和距离维相减
self.data_rangeFFT = self.data_rangeFFT - self.data_rangeFFT.mean(axis=0)
if window_enable: # 加窗
window_coeff_vec = np.hanning(self.rangeFFT_size + 2)[1:-1]
# 使用汉宁窗 【为什么截掉首尾两点 是因为Python汉宁窗首尾为0 需要多加两个点并去掉】
self.data_rangeFFT = self.data_rangeFFT * window_coeff_vec.reshape((-1, 1, 1, 1)) # 窗size 256*1*1*1 广播相乘
self.data_rangeFFT = np.fft.fft(self.data_rangeFFT, n=self.rangeFFT_size, axis=0) # 距离维FFT
if scale_on:
scale_factor = self.scale_factor[int(np.log2(self.rangeFFT_size) - 4)]
self.data_rangeFFT *= scale_factor # 2 * 2 / N的补偿，2是加汉宁窗的恢复系数 2/N是距离维FFT的恢复系数
# scio.savemat('rfft2.mat', {'data': self.data_rangeFFT}) # 写距离fft数据
def __dopplerFFT__(self, window_enable=False, scale_on=False, clutter_remove=False):
self.data_dopplerFFT = self.data_rangeFFT # 获取距离FFT后的数据 256 * 128 * 16 * 12
if window_enable:
window_coeff_vec = np.hanning(self.dopplerFFT_size + 2)[1:-1] # 使用汉宁窗
self.data_dopplerFFT = self.data_dopplerFFT * window_coeff_vec.reshape((1, -1, 1, 1))
if clutter_remove: # MTI 滤除静态杂波
self.data_dopplerFFT = self.data_dopplerFFT - np.expand_dims(self.data_dopplerFFT.mean(axis=1), axis=1)
self.data_dopplerFFT = np.fft.fft(self.data_dopplerFFT, n=self.dopplerFFT_size, axis=1)
self.data_dopplerFFT = np.fft.fftshift(self.data_dopplerFFT, axes=1)
if scale_on:
scale_factor = self.scale_factor[int(np.log2(self.dopplerFFT_size) - 4)]
self.data_rangeFFT *= scale_factor # 这里同样考虑和距离维FFT的增益补偿
def __range_cfar_os__(self):
refWinSize = 8 # 参考单元
guardWinSize = 8 # 保护单元
K0 = 5 # 门限系数 这里设的有点硬 应该自适应调整
discardCellLeft = 10 # 去掉首尾部分数据 可能存在近场泄露
discardCellRight = 20
maxEnable = 0
sortSelectFactor = 0.75 # 有序统计阈值
gaptot = refWinSize + guardWinSize # 滑动窗
n_obj = 0 # 目标数目
index_obj = [] # 目标距离存放矩阵
intensity_obj = [] # 目标强度存放矩阵
noise_obj = [] # 噪声估计存放矩阵
snr_obj = [] # 信噪比存放矩阵
n_range = self.sig_integrate.shape[0] # 256
n_doppler = self.sig_integrate.shape[1] # 128
for i_doppler in range(n_doppler): # 对每个多普勒单元进行遍历
sigv = self.sig_integrate[:, i_doppler] # 获得一维距离信号 256*1
vecMid = sigv[discardCellLeft: n_range - discardCellRight] # 获取截断后的数据 近场信号泄露较强
vecLeft = vecMid[0: gaptot]
vecRight = vecMid[-gaptot:]
vec = np.hstack((vecLeft, vecMid, vecRight)) # 填充截断的数据 为什么要加这一段
# 如果不加这一段的话检测直接从vecMid的第gaptot个单元开始到倒数第gaptot单元 漏掉了两端gaptot的检测 通过填充以后则可以避免这种情况
for j in range(len(vecMid)):
index_cur = j + gaptot
index_left = list(range(index_cur - gaptot, index_cur - guardWinSize))
index_right = list(range(index_cur + guardWinSize + 1, index_cur + gaptot + 1))
sorted_res = np.sort(np.hstack((vec[index_left], vec[index_right])), axis=0)
value_selected = sorted_res[int(np.ceil(sortSelectFactor * len(sorted_res)) - 1)] # 有序选择门限
if maxEnable == 1:
# whether value_selected is the local max value
value_local = vec[index_cur - gaptot: index_cur + gaptot + 1]
value_max = value_local.max()
if vec[index_cur] >= K0 * value_selected and vec[index_cur] >= value_max:
# 既要满足OS又要满足检测单元大于滑动窗局部最大值 意义是什么？这样点不是更少了吗
n_obj += 1
index_obj.append([discardCellLeft + j, i_doppler]) # 将检测到的目标为止进行存储更新 目标位置需要加上截断的部分
intensity_obj.append(vec[index_cur]) # 功率值
noise_obj.append(value_selected) # 更新门限值 在后面被用作了噪声估计值 非相干积累值
snr_obj.append(vec[index_cur] / value_selected) # 更新信噪比 这里用OS的有序选择门限做信噪比是不是不太合理？
else:
if vec[index_cur] >= K0 * value_selected:
n_obj += 1
index_obj.append([discardCellLeft + j, i_doppler])
intensity_obj.append(vec[index_cur])
noise_obj.append(value_selected)
snr_obj.append(vec[index_cur] / value_selected)
self.res_range_cfar = {
'n_obj': n_obj,
'index_obj': np.array(index_obj, dtype='int'),
'intensity_obj': np.array(intensity_obj, dtype='float'),
'noise_obj': np.array(noise_obj, dtype='float'),
'snr_obj': np.array(snr_obj, dtype='float')
}
def __doppler_cfar_os_cyclicity__(self):
refWinSize = 4 # 参考窗大小为4
guardWinSize = 0 # 多普勒保护窗设为0?
K0 = 0.5
maxEnable = 0
sortSelectFactor = 0.75
gaptot = refWinSize + guardWinSize
n_obj = 0
index_obj = []
intensity_obj = []
noise_obj = []
snr_obj = []
index_obj_range = self.res_range_cfar['index_obj'] # 这里要用距离维cfar检测到的点来搜索多普勒维度的目标
index_obj_range_unique = np.unique(index_obj_range[:, 0])
for i_range in index_obj_range_unique:
sigv = self.sig_integrate[i_range, :] # 1*128 获取每个距离单元下的一维多普勒信号
# cyclicity
vecMid = sigv
vecLeft = sigv[-gaptot:]
vecRight = sigv[0: gaptot]
vec = np.hstack((vecLeft, vecMid, vecRight))
for j in range(len(vecMid)):
index_cur = j + gaptot
index_left = list(range(index_cur - gaptot, index_cur - guardWinSize))
index_right = list(range(index_cur + guardWinSize + 1, index_cur + gaptot + 1))
sorted_res = np.sort(np.hstack((vec[index_left], vec[index_right])), axis=0)
value_selected = sorted_res[int(np.ceil(sortSelectFactor * len(sorted_res)) - 1)]
if maxEnable == 1:
# whether value_selected is the local max value
value_local = vec[index_cur - gaptot: index_cur + gaptot + 1]
value_max = value_local.max()
if vec[index_cur] >= K0 * value_selected and vec[index_cur] >= value_max:
if j in index_obj_range[
index_obj_range[:, 0] == i_range, 1]: # 检查当前多普勒单元是否在距离cfar时也被加入了
n_obj += 1
index_obj.append([i_range, j])
intensity_obj.append(vec[index_cur])
noise_obj.append(value_selected)
snr_obj.append(vec[index_cur] / value_selected)
else:
if vec[index_cur] >= K0 * value_selected:
if j in index_obj_range[index_obj_range[:, 0] == i_range, 1]:
n_obj += 1
index_obj.append([i_range, j])
intensity_obj.append(vec[index_cur])
noise_obj.append(value_selected)
snr_obj.append(vec[index_cur] / value_selected)
self.res_doppler_cfar = {
'n_obj': n_obj,
'index_obj': np.array(index_obj, dtype='int'),
'intensity_obj': np.array(intensity_obj, dtype='float'),
'noise_obj': np.array(noise_obj, dtype='float'),
'snr_obj': np.array(snr_obj, dtype='float')
}
def __rdm_cfar__(self): # 真正的距离多普勒 CFAR流程
self.sig_integrate = np.sum(
np.power(
np.abs(
self.data_dopplerFFT.reshape(
(self.data_dopplerFFT.shape[0], self.data_dopplerFFT.shape[1], -1)
)
), 2
),
axis=2
) + 1 # 非相干积累 不太理解的是这里加1？
self.sig_integrate = np.asarray(self.sig_integrate) # 与sig_intergrate共享同块内存
# do CFAR on range doppler map
self.__range_cfar_os__() # 距离维cfar
if self.res_range_cfar['n_obj'] > 0: # 距离维cfar目标检测数目大于0 则进行多普勒为cfar检测
self.__doppler_cfar_os_cyclicity__()
n_obj = self.res_doppler_cfar['n_obj'] # 获取多普勒cfar后得到的目标数目
range_index = self.res_doppler_cfar['index_obj'][:, 0] # 获取目标距离单元索引
range_obj = range_index * self.infos['range_res_m'] # 转换真实距离
doppler_index = self.res_doppler_cfar['index_obj'][:, 1]
doppler_obj = (doppler_index - self.data_dopplerFFT.shape[1] / 2) * self.infos['velocity_res_m_sec'] # 转换速度索引
# use the first noise estimate, apply to the second detection
index1 = self.res_range_cfar['index_obj'] # 获取检测目标的距离
index2 = self.res_doppler_cfar['index_obj'] # 获取检测目标的速度
index_noise = [
np.argwhere(np.logical_and(index1[:, 0] == index2[i, 0], index1[:, 1] == index2[i, 1]))[0][0]
for i in range(index2.shape[0]) # 采用逻辑与判断 遍历多普勒cfar中的每个检测点和距离cfar检测点进行匹配 逻辑1则输出位置
] # 最后输出的长度为多普勒cfar后的点数
noise = self.res_range_cfar['noise_obj'][index_noise] # 获取满足逻辑条件的噪声门限值 长度和多普勒cfar后点数一致
bin_val = self.data_dopplerFFT[range_index, doppler_index, :, :] # 获取cfar后的3217*16*12
intensity = np.power(np.abs(bin_val.reshape((bin_val.shape[0], -1))), 2).sum(axis=1) # 3217*1获取距离维通道的功率并求和
snr = 10 * np.log10(intensity / noise) # 求解信噪比
self.res_rdm_cfar = {
'n_obj': n_obj,
'range_index': range_index,
'range': range_obj,
'doppler_index': doppler_index,
'doppler': doppler_obj,
'doppler_index_origin': doppler_index,
'doppler_origin': doppler_obj,
'noise': noise,
'bin_val': bin_val,
'intensity': intensity,
'snr': snr
}
# 考虑最大速度扩展的多普勒相移补偿
if self.apply_vmax_extend: # 这里考虑的是基于多普勒相偏补偿的速度解模糊方法
doppler_index_unwrap = doppler_index.reshape((-1, 1)) + self.dopplerFFT_size * \
(np.arange(self.num_tx).reshape((1, -1)) - self.num_tx / 2 + 1 * (
doppler_obj <= 0).reshape((-1, 1))) # 3217 * 12
#
doppler_index_unwrap = doppler_index_unwrap.astype(np.int)
sig_bin_org = bin_val # 3217 * 16 * 12 cfar后检测的点
# TDM MIMO导致的多普勒相位校正
delta_phi = 2 * np.pi * (doppler_index_unwrap - self.dopplerFFT_size / 2) / (
self.num_tx * self.dopplerFFT_size) # 3217 * 12
# 根据可能的假设数量，构建所有可能的信号向量
correct_matrix = np.exp(
-1j * np.arange(self.num_tx).reshape((1, 1, -1, 1)) *
np.expand_dims(np.expand_dims(delta_phi, axis=1), axis=1)
) # 3217 * 1 * 12 * 12
sig_bin = np.expand_dims(sig_bin_org, axis=-1) * correct_matrix # 3217 * 16 * 12 * 12
# 使用重叠天线来做最大的多普勒解缠绕
# 解析未冗余阵元和冗余阵元的发射接收阵列
# 为什么只用未冗余阵元和冗余阵元的匹配对来做解缠绕？而不是对发射阵元对应的全部阵元解缠绕？
index_antenna_overlaped_diff1tx = np.stack(
( # 发射天线相位差为一的冗余阵元? 2、3表示未冗余的接收id和发射id;6、7表示冗余的接收和发射id
np.array([np.argwhere(self.rx_id_onboard == i)[0][0] for i in
self.virtual_array_info_overlaped_diff1tx[:, 2]], dtype='int'),
np.array([np.argwhere(self.tx_id_transfer_order == i)[0][0] for i in
self.virtual_array_info_overlaped_diff1tx[:, 3]], dtype='int'),
np.array([np.argwhere(self.rx_id_onboard == i)[0][0] for i in
self.virtual_array_info_overlaped_diff1tx[:, 6]], dtype='int'),
np.array([np.argwhere(self.tx_id_transfer_order == i)[0][0] for i in
self.virtual_array_info_overlaped_diff1tx[:, 7]], dtype='int'),
),
axis=1
) # 32 * 4
# 获取未冗余阵元和冗余阵元对应的信号
sig_overlap = np.stack(
(
sig_bin_org[:, index_antenna_overlaped_diff1tx[:, 0], index_antenna_overlaped_diff1tx[:, 1]],
sig_bin_org[:, index_antenna_overlaped_diff1tx[:, 2], index_antenna_overlaped_diff1tx[:, 3]]
),
axis=2
) # 3217 * 32 * 2
# 检查每个假设的每个重叠天线对的相位差
angle_sum_test = np.zeros((sig_overlap.shape[0], sig_overlap.shape[1], delta_phi.shape[1])) # 3217 * 32 * 12
for i_sig in range(angle_sum_test.shape[1]): # 32
signal2 = np.matmul(
np.expand_dims(sig_overlap[:, :i_sig + 1, 1], axis=2), # 3217 * ？ * 1
np.expand_dims(np.exp(-1j * delta_phi), axis=1) # 3217 * 1 * 12
) # 3217 * 32 * 12
angle_sum_test[:, i_sig, :] = np.angle(
np.sum(
np.expand_dims(sig_overlap[:, :i_sig + 1, 0], axis=2) * signal2.conj(),
axis=1
)
) # 3217 * 32 * 12
# 选择相位差最小的假设来估计解系数
doppler_unwrap_integ_overlap_index = np.argmin(np.abs(angle_sum_test), axis=2) # 3217 * 12
# 角度维FFT信噪比？
index_antenna_noredundant_row1 = np.stack(
(
np.array([np.argwhere(self.rx_id_onboard == i)[0][0] for i in
self.virtual_array_noredundant_row1[:, 2]], dtype='int'), # 获取接收位置
np.array([np.argwhere(self.tx_id_transfer_order == i)[0][0] for i in
self.virtual_array_noredundant_row1[:, 3]], dtype='int') # 获取发射位置
),
axis=1
) # 86 * 2
# scio.savemat('array.mat', {'data': index_antenna_noredundant_row1}) # 写通道排列顺序
# 3217 * 86 * 12 获取方位维度的天线排列
sig_bin_row1 = sig_bin[:, index_antenna_noredundant_row1[:, 0], index_antenna_noredundant_row1[:, 1], :]
angleFFT_size = 128
sig_bin_row1_fft = np.fft.fftshift(
np.fft.fft(sig_bin_row1, n=angleFFT_size, axis=1),
axes=1
) # 3217 * 128 * 12 这里做了角度维FFT
angle_bin_skip_left = 4
angle_bin_skip_right = 4
sig_bin_row1_fft_cut = np.abs(
sig_bin_row1_fft[:, angle_bin_skip_left:angleFFT_size - angle_bin_skip_right, :]
)
doppler_unwrap_integ_FFT_index = np.argmax(
np.max(sig_bin_row1_fft_cut, axis=1),
axis=1
)
doppler_unwrap_integ_index = np.array([
np.argmax(np.bincount(doppler_unwrap_integ_overlap_index[i, :]))
for i in range(doppler_unwrap_integ_overlap_index.shape[0])
]) # 3217
self.res_rdm_cfar['bin_val_correct'] = sig_bin[np.arange(sig_bin.shape[0]), :, :,
doppler_unwrap_integ_index] # 3217 * 16 * 12 * 12
# corret doppler after applying the integer value
doppler_index_correct = doppler_index_unwrap[
np.arange(doppler_index_unwrap.shape[0]), doppler_unwrap_integ_index] # 3217 * 1
doppler_index_correct_FFT = doppler_index_unwrap[
np.arange(doppler_index_unwrap.shape[0]), doppler_unwrap_integ_FFT_index] # 3217 * 1
doppler_index_correct_overlap = doppler_index_unwrap[
np.arange(doppler_index_unwrap.shape[0]), doppler_unwrap_integ_index] # 3217 * 1
self.res_rdm_cfar['doppler_index_correct'] = doppler_index_correct
self.res_rdm_cfar['doppler_correct'] = (doppler_index_correct - self.data_dopplerFFT.shape[1] / 2) * \
self.infos['velocity_res_m_sec']
self.res_rdm_cfar['doppler_index_correct_FFT'] = doppler_index_correct_FFT
self.res_rdm_cfar['doppler_correct_FFT'] = (doppler_index_correct_FFT - self.data_dopplerFFT.shape[
1] / 2) * \
self.infos['velocity_res_m_sec']
self.res_rdm_cfar['doppler_index_correct_overlap'] = doppler_index_correct_overlap
self.res_rdm_cfar['doppler_correct_overlap'] = (doppler_index_correct_overlap -
self.data_dopplerFFT.shape[1] / 2) * \
self.infos['velocity_res_m_sec']
# replace doppler_index: doppler_index_correct, doppler: doppler_correct
self.res_rdm_cfar['doppler_index'] = self.res_rdm_cfar['doppler_index_correct']
self.res_rdm_cfar['doppler'] = self.res_rdm_cfar['doppler_correct']
# 最大速度扩展的最小距离约束
index_noapply = np.argwhere(range_obj <= self.min_dis_apply_vmax_extend)[:, 0] # 438
# 未使用速度扩展的多普勒相位纠正？
delta_phi_noapply = (
2 * np.pi * (doppler_index - self.dopplerFFT_size / 2) / (self.num_tx * self.dopplerFFT_size) # 3217 * 1
).reshape((-1, 1)) # 3217 * 1
i_tx_noapply = np.arange(self.num_tx).reshape((1, -1)) # 1 * 12
sig_bin_noapply = bin_val * np.exp( # 394*16*12
-1j * np.matmul(delta_phi_noapply, i_tx_noapply)
).reshape((delta_phi_noapply.shape[0], 1, i_tx_noapply.shape[1])) # 3217 * 16 * 12
self.res_rdm_cfar['bin_val_correct'][index_noapply, :, :] = sig_bin_noapply[index_noapply, :, :] # 3217 * 16 * 12
self.res_rdm_cfar['doppler_index_correct'][index_noapply] = self.res_rdm_cfar['doppler_index_origin'][
index_noapply]
self.res_rdm_cfar['doppler_correct'][index_noapply] = self.res_rdm_cfar['doppler_origin'][index_noapply]
self.res_rdm_cfar['doppler_index_correct_FFT'][index_noapply] = \
self.res_rdm_cfar['doppler_index_origin'][index_noapply]
self.res_rdm_cfar['doppler_correct_FFT'][index_noapply] = self.res_rdm_cfar['doppler_origin'][
index_noapply]
self.res_rdm_cfar['doppler_index_correct_overlap'][index_noapply] = \
self.res_rdm_cfar['doppler_index_origin'][index_noapply]
self.res_rdm_cfar['doppler_correct_overlap'][index_noapply] = self.res_rdm_cfar['doppler_origin'][
index_noapply]
self.res_rdm_cfar['doppler_index'][index_noapply] = self.res_rdm_cfar['doppler_index_origin'][
index_noapply]
self.res_rdm_cfar['doppler'][index_noapply] = self.res_rdm_cfar['doppler_origin'][index_noapply]
else:
# 未考虑最大速度扩展的多普勒相移补偿
sig_bin_org = bin_val
delta_phi = (
2 * np.pi * (doppler_index - self.data_dopplerFFT.shape[1] / 2) / \
(self.data_dopplerFFT.shape[3] * self.data_dopplerFFT.shape[1])
).reshape((-1, 1))
i_tx = np.arange(self.data_dopplerFFT.shape[3]).reshape((1, -1))
correct_matrix = np.exp(-1j * np.matmul(delta_phi, i_tx)).reshape(
(delta_phi.shape[0], 1, i_tx.shape[1])
)
sig_bin_correct = sig_bin_org * correct_matrix
self.res_rdm_cfar['bin_val_correct'] = sig_bin_correct
else:
self.res_rdm_cfar = None
def __single_obj_signal_space_mapping__(self, sig):
sig_space = np.zeros(
(self.virtual_array_azimuth.max() + 1, self.virtual_array_elevation.max() + 1),
dtype=sig.dtype
) # 获取方位维度和俯仰维度的最大孔径 初始化
# 获取非冗余阵列的方位、俯仰维度的虚拟阵元
sig_space_index0 = self.virtual_array_noredundant[:, 0] # 134
sig_space_index1 = self.virtual_array_noredundant[:, 1] # 134
sig_index0 = np.array(
[np.argwhere(self.rx_id_onboard == i_rx)[0][0] for i_rx in self.virtual_array_noredundant[:, 2]]) # 接收阵元位置
# np.argwhere 返回行列坐标 1*16 返回对应接收天线位置
sig_index1 = np.array(
[np.argwhere(self.tx_id_transfer_order == i_tx)[0][0] for i_tx in self.virtual_array_noredundant[:, 3]]) # 发射阵元位置
sig_space[sig_space_index0, sig_space_index1] = sig[sig_index0, sig_index1] # 获取重排后的阵元空间
return sig_space
def __DOA_BF_PeakDet_loc__(self, input_data, sidelobeLevel_dB):
# 类似于findpeaks函数？
gamma = np.power(10, 0.02)
min_val = np.inf
max_val = -np.inf
max_loc = 0
max_data = []
locate_max = 0
num_max = 0
extend_loc = 0
init_stage = 1
abs_max_value = 0
i = -1
N = len(input_data) # 输入数据长度 256
while (i < (N + extend_loc - 1)):
i += 1
i_loc = np.mod(i, N)
current_val = input_data[i_loc] # 当前值
if current_val > abs_max_value:
abs_max_value = current_val
if current_val > max_val:
max_val = current_val
max_loc = i_loc
max_loc_r = i
if current_val < min_val:
min_val = current_val
if locate_max == 1:
if current_val < max_val / gamma:
num_max += 1
bwidth = i - max_loc_r
max_data.append([max_loc, max_val, bwidth, max_loc_r])
min_val = current_val
locate_max = 0
else:
if current_val > min_val * gamma:
locate_max = 1
max_val = current_val
max_loc = i_loc
max_loc_r = i
if init_stage == 1:
extend_loc = i
init_stage = 0
max_data = np.array(max_data)
if len(max_data.shape) < 2:
max_data = np.zeros((0, 4))
max_data = max_data[max_data[:, 1] >= abs_max_value * pow(10, -sidelobeLevel_dB / 10), :] # 约束了峰值条件
peak_val = max_data[:, 1]
peak_loc = max_data[:, 0].astype('int')
return peak_val, peak_loc
def __DOA_beamformingFFT_2D__(self, sig):
sidelobeLevel_dB_azim = 1
sidelobeLevel_dB_elev = 0
doa_fov_azim = [-70, 70] # 方位维度的FOV
doa_fov_elev = [-20, 20] # 俯仰维度的FOV
data_space = self.__single_obj_signal_space_mapping__(sig) # 返回86*7 7表示的是俯仰维度的空间
data_azimuthFFT = np.fft.fftshift(np.fft.fft(data_space, n=self.azimuthFFT_size, axis=0), axes=0) # 方位维度的FFT 256*7
data_elevationFFT = np.fft.fftshift(np.fft.fft(data_azimuthFFT, n=self.elevationFFT_size, axis=1), axes=1) # 俯仰维度的FFT 256*256
spec_azim = np.abs(data_azimuthFFT[:, 0]) # 256
_, peak_loc_azim = self.__DOA_BF_PeakDet_loc__(spec_azim, sidelobeLevel_dB_azim) # 返回检测到的方位维度的峰值索引
n_obj = 0
azimuth_obj = []
elevation_obj = []
azimuth_index = []
elevation_index = []
for i in range(len(peak_loc_azim)): # 对每个峰值进行遍历
spec_elev = abs(data_elevationFFT[peak_loc_azim[i], :]) # 256, 俯仰维度的空间谱
peak_val_elev, peak_loc_elev = self.__DOA_BF_PeakDet_loc__(spec_elev, sidelobeLevel_dB_elev)
# 获取俯仰维度的峰值和索引
for j in range(len(peak_loc_elev)): # 对俯仰维度的峰值索引进行遍历
est_azimuth = np.arcsin(
self.azimuth_bins[peak_loc_azim[i]] / 2 / np.pi / self.doa_unitDis) / np.pi * 180 # 方位角解析
est_elevation = np.arcsin(
self.elevation_bins[peak_loc_elev[j]] / 2 / np.pi / self.doa_unitDis) / np.pi * 180 # 俯仰角解析
if est_azimuth >= doa_fov_azim[0] and est_azimuth <= doa_fov_azim[1] and est_elevation >= doa_fov_elev[
0] and est_elevation <= doa_fov_elev[1]: # 如果测得的角度在俯仰和方位FOV内则认为正确 加入目标角度
n_obj += 1
azimuth_obj.append(est_azimuth) #
elevation_obj.append(est_elevation)
azimuth_index.append(peak_loc_azim[i])
elevation_index.append(peak_loc_elev[j])
# 更新目标角度和峰值索引
azimuth_obj = np.array(azimuth_obj, dtype='float')
elevation_obj = np.array(elevation_obj, dtype='float')
azimuth_index = np.array(azimuth_index, dtype='int')
elevation_index = np.array(elevation_index, dtype='int')
res = {
'data_space': data_space,
'data_azimuthFFT': data_azimuthFFT,
'data_elevationFFT': data_elevationFFT,
'n_obj': n_obj,
'azimuth_obj': azimuth_obj,
'elevation_obj': elevation_obj,
'azimuth_index': azimuth_index,
'elevation_index': elevation_index
}
return res
def __doa__(self):
n_obj = 0
range_index = []
range_val = []
doppler_index = []
doppler = []
if self.apply_vmax_extend:
doppler_index_origin = []
doppler_origin = []
azimuth_index = []
azimuth = []
elevation_index = []
elevation = []
snr = []
intensity = []
noise = []
scio.savemat('cfar_data.mat', {'data': self.res_rdm_cfar})
for i in range(self.res_rdm_cfar['n_obj']):
x = self.res_rdm_cfar['bin_val_correct'][i, :, :] # 遍历每个目标的空域信号 16*12
res_DOA_beamformingFFT_2D = self.__DOA_beamformingFFT_2D__(x) # 返回doa估计的内容
for j in range(res_DOA_beamformingFFT_2D['n_obj']): # 遍历估计目标 可能会估计出多个目标 这个时候需要更新下标
n_obj += 1 # 信源数统计
range_index.append(self.res_rdm_cfar['range_index'][i]) # 获取距离单元索引
range_val.append(self.res_rdm_cfar['range'][i]) # 获取距离值
doppler_index.append(self.res_rdm_cfar['doppler_index'][i]) # 获取多普勒单元索引
doppler.append(self.res_rdm_cfar['doppler'][i]) # 获取多普勒速度
if self.apply_vmax_extend:
doppler_index_origin.append(self.res_rdm_cfar['doppler_index_origin'][i])
doppler_origin.append(self.res_rdm_cfar['doppler_origin'][i])
azimuth_index.append(res_DOA_beamformingFFT_2D['azimuth_index'][j]) # 方位角峰值索引
azimuth.append(res_DOA_beamformingFFT_2D['azimuth_obj'][j]) # 方位角目标角度
elevation_index.append(res_DOA_beamformingFFT_2D['elevation_index'][j]) # 俯仰角峰值索引
elevation.append(res_DOA_beamformingFFT_2D['elevation_obj'][j]) # 俯仰角目标角度
snr.append(self.res_rdm_cfar['snr'][i]) # 信噪比
intensity.append(self.res_rdm_cfar['intensity'][i]) # 强度积累
noise.append(self.res_rdm_cfar['noise'][i]) # 噪声门限
if n_obj > 0:
self.res_doa = {
'n_obj': n_obj,
'range_index': np.array(range_index, dtype='int'),
'range': np.array(range_val, dtype='float'),
'doppler_index': np.array(doppler_index, dtype='int'),
'doppler': np.array(doppler, dtype='float'),
'azimuth_index': np.array(azimuth_index, dtype='int'),
'azimuth': np.array(azimuth, dtype='float'),
'elevation_index': np.array(elevation_index, dtype='int'),
'elevation': np.array(elevation, dtype='float'),
'snr': np.array(snr, dtype='float'),
'intensity': np.array(intensity, dtype='float'),
'noise': np.array(noise, dtype='float')
}
if self.apply_vmax_extend:
self.res_doa['doppler_index_origin'] = np.array(doppler_index_origin, dtype='int')
self.res_doa['doppler_origin'] = np.array(doppler_origin, dtype='float')
else:
self.res_doa = None
def __generate_pcd__(self):
r = self.res_doa['range']
azi = self.res_doa['azimuth']
ele = self.res_doa['elevation']
x = r * np.cos(ele / 180 * np.pi) * np.sin(azi / 180 * np.pi)
y = r * np.cos(ele / 180 * np.pi) * np.cos(azi / 180 * np.pi)
z = r * np.sin(ele / 180 * np.pi)
doppler = self.res_doa['doppler']
snr = self.res_doa['snr']
intensity = self.res_doa['intensity']
self.pcd = {
'x': x,
'y': y,
'z': z,
'doppler': doppler,
'snr': snr,
'intensity': intensity
}
# scio.savemat('pcd.mat', {'data': self.pcd})
def __generate_heatmapBEV__(self, doppler_correction=False):
sig = self.data_dopplerFFT
if doppler_correction:
# add Doppler correction before generating the heatmap
delta_phi = 2 * np.pi * (np.arange(self.dopplerFFT_size) - self.dopplerFFT_size / 2) / \
(self.num_tx * self.dopplerFFT_size)
cor_vec = np.exp(-1j * delta_phi.reshape((-1, 1)) * np.arange(self.num_tx).reshape((1, -1)))
sig = sig * np.expand_dims(np.expand_dims(cor_vec, axis=1), axis=0)
index_antenna_noredundant_row1 = np.stack(
(
np.array(
[np.argwhere(self.rx_id_onboard == i)[0][0] for i in self.virtual_array_noredundant_row1[:, 2]],
dtype='int'),
np.array([np.argwhere(self.tx_id_transfer_order == i)[0][0] for i in
self.virtual_array_noredundant_row1[:, 3]], dtype='int')
),
axis=1
)
sig_row1 = sig[:, :, index_antenna_noredundant_row1[:, 0], index_antenna_noredundant_row1[:, 1]]
sig_row1_azimuthFFT = np.fft.fftshift(
np.fft.fft(sig_row1, n=self.azimuthFFT_size, axis=2),
axes=2
) # 方位维FFT
#
heatmapBEV_static = np.abs(sig_row1_azimuthFFT[:, self.dopplerFFT_size // 2, :]) # 256*128*256
heatmapBEV_static = np.hstack((heatmapBEV_static, heatmapBEV_static[:, 0:1]))
heatmapBEV_dynamic = np.sum(
np.abs(sig_row1_azimuthFFT[:, np.arange(self.dopplerFFT_size) != self.dopplerFFT_size // 2, :]),
axis=1
)
heatmapBEV_dynamic = np.hstack((heatmapBEV_dynamic, heatmapBEV_dynamic[:, 0:1]))
sine_theta = np.arange(1, -1 - 2 / self.azimuthFFT_size, -2 / self.azimuthFFT_size)
cos_theta = np.sqrt(1 - np.power(sine_theta, 2))
x = np.matmul(
self.range_bins.reshape((-1, 1)),
sine_theta.reshape((1, -1))
)
y = np.matmul(
self.range_bins.reshape((-1, 1)),
cos_theta.reshape((1, -1))
)
self.heatmapBEV = {
'heatmapBEV_static': heatmapBEV_static,
'heatmapBEV_dynamic': heatmapBEV_dynamic,
'x': x,
'y': y
}
# scio.savemat('heatmap.mat', {'data': self.heatmapBEV})
import numpy as np
import scipy.io as scio
from signal_process import NormalModeProcess
def load_TIRadar_npz(path):
data = np.load(path, allow_pickle=True)
res = {
'data_imag': data['data_imag'],
'data_real': data['data_real'],
'mode_infos': data['mode_infos'][()]
}
return res
def load_TIRadar_calibmat(path):
data = scio.loadmat(path)
calibResult = dict()
names = data['calibResult'].dtype.names
values = data['calibResult'].tolist()[0][0]
for i in range(len(values)):
calibResult[names[i]] = values[i]
calibResult['RangeMat'] = calibResult['RangeMat'].astype(np.int)
params = dict()
names = data['params'].dtype.names
values = data['params'].tolist()[0][0]
for i in range(len(values)):
params[names[i]] = values[i][0, 0]
res = {
'calibResult': calibResult,
'params': params
}
return res
def main():
'''
path_data = 'data/data1/1671267290.996.npz'
path_calib = 'data/data1/mode1_256x128.mat'
path_json = 'data/data1/1671267290.996.json'
'''
path_data = 'data/data2/1671267291.096.npz'
path_calib = 'data/data2/mode1_256x128.mat'
path_json = 'data/data2/1671267291.096.json'
radar_data = load_TIRadar_npz(path_data)
radar_calib = load_TIRadar_calibmat(path_calib)
nmp = NormalModeProcess(radar_data['mode_infos'], radar_data['data_real'], radar_data['data_imag'], radar_calib)
nmp.run(generate_pcd=True, generate_heatmapBEV=True, generate_heatmap4D=False)
pcd = nmp.get_pcd()
print(pcd)
if __name__ == '__main__':
main()
参考文献
[1] X. Yu, Z. Cao, Z. Wu, C. Song, J. Zhu and Z. Xu, "A Novel Potential Drowning Detection System Based on Millimeter-Wave Radar," 2022 17th International Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, Singapore, 2022, pp. 659-664, doi: 10.1109/ICARCV57592.2022.10004245.
[2] J. Fuchs, M. Gardill, M. Lübke, A. Dubey and F. Lurz, "A Machine Learning Perspective on Automotive Radar Direction of Arrival Estimation," in IEEE Access, vol. 10, pp. 6775-6797, 2022, doi: 10.1109/ACCESS.2022.3141587.
[3] 初出茅庐：毫米波雷达开发手册之基础参数 https://www.cnblogs.com/yuxuliang/p/MyRadar_4.html
[4] https://ieee-aess.org/presentation/webinar/introduction-mimo-radar-waveforms
