社区共建 | 北京理工大学黄涛基于 A100 实现 MATLAB 全链路信号处理工程

作者介绍

黄涛 北京理工大学毫米波与太赫兹波技术在读博士生，研究方向太赫兹光电子学

社区贡献说明

本文为 SuperRadar 开发者贡献共建应用案例，项目源于 CTSAI A100 Issue #7。

贡献者基于 CTSAI A100 原始 ADC 数据，建立了一套可复现、可继续二次开发的 MATLAB 雷达信号处理工程。

项目工程位于：https://github.com/Super-Radar/CTSAI-A100/tree/main/community/projects/a100-matlab-adc-signal-processin

引言

毫米波雷达最终输出一个目标的距离、速度和角度之前，中间究竟经历了什么？

本项目从 CTSAI A100 ADC 采集工具导出的 TXT 原始数据开始，建立了一条完整的 MATLAB 信号处理链：

• 从头字段检查、32-bit packed word 解包开始，建立 Range FFT 与 Doppler FFT 的物理坐标；
• 随后加入慢时间杂波抑制、二维 CFAR、自适应 VI-CFAR、二维亚栅格距离—速度细化；
• 最后利用真实阵列位置和相位校准执行 Angle FFT、DML、MUSIC 与 OMP，输出距离、径向速度、方位角、SNR、CSV/MAT 结果以及目标点云。

后来者不需要再从 ADC 文件格式开始重复搭建基础链路，而是可以直接基于现有工程替换 CFAR、增加新的杂波抑制方法、研究不同 DOA 算法，甚至进一步探索 AI 对 ADC、Range-Doppler 或点云数据的处理。

完整处理链

工程整体处理流程如下：

load_ctsaia100_config
        ↓
derive_radar_parameters
        ↓
load_adc_dataset / read_adc_txt / unpack_uint32_adc
        ↓
range_fft
        ↓
organize_virtual_array
        ↓
suppress_clutter
        ↓
doppler_fft
        ↓
cfar_2d / vi_cfar_map
        ↓
extract_detections
        ↓
refine_detections_subbin
        ↓
estimate_angles_for_detections
        ↓
Angle FFT / DML / MUSIC / OMP
        ↓
CSV / MAT / figures / target point cloud

整个项目由 main.m 作为统一入口，各个关键处理环节被拆分为独立函数。

ADC 数据解析

整条处理链首先要解决的，是如何正确读取 CTSAI A100 输出的 ADC 原始数据。

每个 TXT 文件前三个字段依次表示：

• 接收通道编号；
• 每个 chirp 的快时间采样点数；
• chirp 数量。

以 Pf0 数据为例：

N_s=1024,\qquad N_c=256

每一个 32-bit packed word 包含两个 16-bit ADC 样本，因此标准payload 中 packed word 的数量为：

N_{\mathrm{word}}=\frac{N_sN_c}{2}=\frac{1024\times256}{2}=131072

如果读取出的一个无符号 32-bit 数据字记作 w，其高、低 16 位分别为：

u_H=\left\lfloor\frac{w}{2^{16}}\right\rfloor

u_L=w-2^{16}u_H

ADC 数据使用 16-bit 二进制补码。对于：

u\in[0,2^{16}-1]

恢复为 signed int16 时：

s(u)= \begin{cases} u, & 0\le u<2^{15},\\ u-2^{16}, & 2^{15}\le u<2^{16}. \end{cases}

因此，一个 packed word 最终被还原为两个连续采样：

w\longrightarrow\left[s(u_H),s(u_L)\right]

按照文件中确认的样本顺序重新排列后，每个 RX 通道可以写成：

\mathbf{X}_q\in\mathbb{R}^{N_s\times N_c}

多接收通道进一步组成：

x[n,m,q]

其中对应关系如下：

• n：fast-time sample；
• m：chirp；
• q：RX channel；

先确保输入数据正确，再开始信号处理。

LFMCW 信号模型

在进入 FFT 之前，可以先从 LFMCW 信号模型理解距离和速度信息是如何进入 ADC 数据的。

发射 chirp

设毫米波雷达载频为 f_c，有效调频带宽为 B，chirp 上升沿时间为 T_r。

线性调频斜率为：

S=\frac{B}{T_r}

在一个 chirp 内，瞬时频率可以表示为：

f_{\mathrm{tx}}(t)=f_c+St,\qquad 0\le t\lt T_r

对应复解析发射信号：

s_{\mathrm{tx}}(t) = A_t \exp \left\{ j2\pi \left( f_ct+\frac12St^2 \right) \right\}.

工程中的教学波形仿真如果搬移到复基带，则可以使用：

s_{\mathrm{tx,bb}}(t)=\exp\left(j\pi St^2\right)

直接生成 chirp。

运动目标回波

设目标初始距离为 R_0，径向速度为 v。

在较短的观测区间内，可以采用匀速模型：

R(t)=R_0+vt

对应往返传播时延：

\tau(t)=\frac{2R(t)}{c}=\tau_0+\frac{2v}{c}t

其中：

\tau_0=\frac{2R_0}{c}

忽略高阶传播效应后，接收信号可以表示为：

s_{\mathrm{rx}}(t) = A_r \exp \left\{ j2\pi \left[ f_c(t-\tau) + \frac12S(t-\tau)^2 \right] \right\}.

将发射信号和延迟回波去斜，例如采用：

s_b(t)=s_{\mathrm{tx}}(t)s_{\mathrm{rx}}^{*}(t)

在目标单 chirp 内位移较小的情况下，忽略高阶项，beat frequency 的主要组成可以近似表示为：

f_b\approx S\tau_0+f_D

其中 Doppler frequency 为：

f_D=\frac{2v}{\lambda}

波长为：

\lambda=\frac{c}{f_c}

代入 \tau_0=2R_0/c，得到：

f_b\approx\frac{2SR_0}{c}+\frac{2v}{\lambda}

需要注意的是，Doppler 项的正负号会受到去斜共轭方向和速度符号定义影响，但距离项与速度项的量级关系不会因此改变。

本工程采用的速度方向约定为：

v<0:\ \text{approaching}, \qquad

表示目标靠近。

v>0:\ \text{receding}.

表示目标远离。

如果距离处理时只保留调频对应的距离频率：

f_R=\frac{2SR}{c}

那么：

R=\frac{cf_R}{2S}

实际 beat 中仍包含 f_D，直接把 f_b 全部解释为距离会产生

\hat R = \frac{cf_b}{2S} \approx R + \frac{cf_D}{2S}.

利用 f_D=2v/\lambda：

\hat R-R \approx \frac{cv}{S\lambda}

因此，如果实际 beat frequency 中仍然包含Doppler项，却将全部频率直接解释为距离，就可能引入一定的距离偏差。

对于常见低速目标，这一项通常小于主要距离项，但在高速目标或高精度处理中仍需要考虑。

从连续信号到 Range-Doppler 数据

完成 ADC 采样以后，一个 frame 可以近似表示为二维 fast-time / slow-time 数据。

对于第 p 个理想点目标：

x_p[n,m]=A_p\exp\left(j2\pi f_{R,p}\frac{n}{f_s}\right)\exp\left(j2\pi f_{D,p}mT_{\mathrm{sameTX}}\right)

其中：

f_{R,p}=\frac{2SR_p}{c}

f_{D,p}=\frac{2v_p}{\lambda}

多个目标和噪声叠加以后：

x[n,m]=\sum_{p=1}^{P}A_p\exp\left(j2\pi f_{R,p}\frac{n}{f_s}\right)\exp\left(j2\pi f_{D,p}mT_{\mathrm{sameTX}}\right)+w[n,m]

从这个表达式可以看到：快时间相位主要承载距离信息，慢时间 chirp-to-chirp 相位主要承载径向速度信息。

因此，可以沿两个维度分别执行 FFT。

距离维 FFT

对固定 chirp m 和接收通道 q，距离维处理可以表示为：

X_r[k,m,q]=\sum_{n=0}^{N_s-1}x[n,m,q]w_r[n]\exp\left(-j2\pi\frac{kn}{N_r}\right)

其中：

• w_r[n]：距离维窗函数；
• N_r：Range FFT 点数。

窗函数主要用于降低有限采样截断产生的旁瓣。

若需要保持不同窗之间的幅度尺度一致，可以使用 coherent gain：

G_r=\sum_{n=0}^{N_s-1}w_r[n]

进行归一化：

\widetilde{X}_r[k,m,q]=\frac{X_r[k,m,q]}{G_r}

FFT 第 k 个正频率 Bin 对应：

f_k=\frac{k}{N_r}f_s

结合 FMCW 距离与beat frequency的关系关系：

R=\frac{cf}{2S}

可以得到物理距离轴：

R_k=\frac{cf_s}{2SN_r}k

因此，相邻 FFT 栅格的距离间隔为：

\Delta R_{\mathrm{bin}}=\frac{cf_s}{2SN_r}

如果 beat signal 的正频率使用 Nyquist 区间：

0\le f_b<\frac{f_s}{2}

则最大无混叠距离近似为：

R_{\max}=\frac{cf_s}{4S}

Pf0 参数下的距离轴

Pf0 主要参数如下：

调频斜率：

S=\frac{300\times10^6}{43\times10^{-6}}\approx6.9767\times10^{12}\ \mathrm{Hz/s}

代入距离栅格公式：

\Delta R_{\mathrm{bin}}=\frac{299792458\times25\times10^6}{2\times6.9767\times10^{12}\times1024}\approx0.52454\ \mathrm{m/bin}

对应 Nyquist 距离上限约为：

R_{\max}\approx268.56\ \mathrm{m}

距离分辨率与 FFT 栅格

这里还需要区分两个很容易混淆的概念：

FFT 距离栅格，并不完全等于雷达理论距离分辨率。

完整调频带宽决定经典 FMCW 理论距离分辨率：

\Delta R_{\mathrm{res}}=\frac{c}{2B}

当：

B=300\ \mathrm{MHz}

理论距离分辨率约为：

\Delta R_{\mathrm{res}}\approx0.49965\ \mathrm{m}

而当前 ADC 每个 chirp 的实际观测时间：

T_{\mathrm{obs}}=\frac{N_s}{f_s}=\frac{1024}{25\times10^6}=40.96\ \mu\mathrm{s}

在这一采样窗口内实际覆盖的线性调频带宽约为：

B_{\mathrm{obs}}=ST_{\mathrm{obs}}\approx285.77\ \mathrm{MHz}

因此：

\frac{c}{2B_{\mathrm{obs}}}\approx0.52454\ \mathrm{m}

与当前 Range FFT 栅格基本一致。

所以：约 0.5245 m/bin 描述的是当前采样窗口和 FFT 参数对应的距离坐标间隔，而约 0.4997 m 对应完整 300 MHz 调频带宽下的理论距离分辨能力。

两者数值接近，但物理含义不同。

慢时间 Doppler FFT

完成距离维处理以后，可以对固定距离 Bin 和通道继续观察连续 chirp。

如果目标径向速度近似恒定，其慢时间相位满足：

x[m]\propto\exp\left(j2\pi f_DmT_{\mathrm{sameTX}}\right)

这里使用 T_{\mathrm{sameTX}}，是为了兼容后续 TDM MIMO 场景：

当不同 TX 轮流发射时，同一个 TX 的有效慢时间采样周期可能是多个 chirp period。

Doppler FFT 可以写成：

X_D[k,\ell,q]=\sum_{m=0}^{M-1}X_r[k,m,q]w_d[m]\exp\left(-j2\pi\frac{\ell m}{N_d}\right)

其中：

• M：实际 slow-time 样本数；
• N_d：Doppler FFT 点数。

执行 fftshift 后，零 Doppler 位于频谱中心。

多普勒频率与径向速度满足：

v=\frac{\lambda f_D}{2}

Doppler FFT 的频率栅格：

\Delta f_{D,\mathrm{bin}}=\frac{1}{N_dT_{\mathrm{sameTX}}}

对应速度栅格：

\Delta v_{\mathrm{bin}}=\frac{\lambda}{2N_dT_{\mathrm{sameTX}}}

若实际 slow-time 采集 M 个 chirp，则 coherent processing interval 为：

T_{\mathrm{CPI}}=MT_{\mathrm{sameTX}}

由有限观测时长决定的速度分辨尺度近似为：

\Delta v_{\mathrm{res}}\approx\frac{\lambda}{2T_{\mathrm{CPI}}}=\frac{\lambda}{2MT_{\mathrm{sameTX}}}

因此，当 N_d>M 时，Doppler 零填充可以使频谱栅格变得更加密集，但不会增加 CPI，也就不会等比例提高真实速度分辨能力。

这和距离维中的“FFT 栅格”和“物理分辨率”区别类似。

慢时间采样频率：

f_{s,\mathrm{slow}}=\frac{1}{T_{\mathrm{sameTX}}}

因此无模糊 Doppler 范围近似为：

|f_D|<\frac{1}{2T_{\mathrm{sameTX}}}

转换成速度：

|v| < v_{\max} = \frac{\lambda}{4T_{\mathrm{sameTX}}} .

Pf0 SISO 配置下：

T_{\mathrm{sameTX}}=48\ \mu\mathrm{s}

由：

f_c=76.3\ \mathrm{GHz}

得到：

\lambda=\frac{c}{f_c}\approx3.9291\ \mathrm{mm}

因此：

\Delta v_{\mathrm{bin}}\approx0.15988\ \mathrm{m/s/bin}

v_{\max}\approx20.46\ \mathrm{m/s}

慢时间杂波抑制

真实道路和室内环境中，固定反射面在连续 chirp 之间的位置变化很小，因此其慢时间相位也相对稳定，大量能量会集中在：f_D\approx0 区域附近。

车辆、人体和其他运动目标，则会通过 chirp-to-chirp 相位变化产生非零 Doppler。

本项目在 Doppler FFT 之前提供多种慢时间杂波处理方式：

NONE
MEAN
MTI2
MTI3
SVD
SVD_MEAN
SVD_MTI2

不同方法对“稳定背景”的数学假设并不相同：

• MEAN：直接去除常量慢时间分量；
• MTI：利用差分滤波器在零频构造零点；
• SVD：从整帧数据中寻找能量占主导的低秩相关子空间。

均值消除

对于某个 Range-Channel 单元的慢时间序列：

\mathbf{x}=\begin{bmatrix}x_0&x_1&\cdots&x_{M-1}\end{bmatrix}^{T}

慢时间均值：

\bar{x}=\frac{1}{M}\mathbf{1}^{T}\mathbf{x}

处理后：

\mathbf{y}=\mathbf{x}-\bar{x}\mathbf{1}

也可以写成：

\mathbf{y}=\left(\mathbf{I}-\frac{1}{M}\mathbf{1}\mathbf{1}^{T}\right)\mathbf{x}

记：

\mathbf{P}_{\perp}=\mathbf{I}-\frac{1}{M}\mathbf{1}\mathbf{1}^{T}

由于：

\mathbf{P}_{\perp}\mathbf{1}=\mathbf{0}

任何在一个 CPI 内保持严格常量的慢时间复分量都会被投影掉。

但如果目标速度非常低，一个 CPI 内相位变化同样较小，那么目标分量也可能与常量子空间高度相关。

因此，均值消除并不只会影响“背景”。

MTI2

二脉冲抵消器沿慢时间执行：

y[n]=x[n]-x[n-1]

对应 Z 域传递函数：

H_{\mathrm{MTI2}}(z)=1-z^{-1}

令：

z=e^{j\omega}

得到频率响应：

H_{\mathrm{MTI2}}\left(e^{j\omega}\right)=1-e^{-j\omega}

幅频响应：

\left|H_{\mathrm{MTI2}}\left(e^{j\omega}\right)\right|=2\left|\sin\left(\frac{\omega}{2}\right)\right|

运动目标对应的慢时间数字角频率为：

\omega_D=2\pi f_DT_{\mathrm{sameTX}}=\frac{4\pi vT_{\mathrm{sameTX}}}{\lambda}

因此 MTI2 对不同径向速度的幅度响应可以表示为：

A_2(v)=2\left|\sin\left(\frac{2\pi vT_{\mathrm{sameTX}}}{\lambda}\right)\right|

当目标速度接近零时：

A_2(v)\approx\frac{4\pi T_{\mathrm{sameTX}}}{\lambda}|v|

这也给出了低速目标衰减的直接来源：目标速度越接近零，慢时间频率就越靠近 MTI 的零点。

MTI3

代码中的 MTI3 使用三个相邻脉冲：

y[n]=x[n]-2x[n-1]+x[n-2]

对应：

H_{\mathrm{MTI3}}(z)=1-2z^{-1}+z^{-2}

这里的 “3” 表示三脉冲结构，对应的是二阶差分。

其幅频响应为：

\left|H_{\mathrm{MTI3}}\left(e^{j\omega}\right)\right|=4\sin^2\left(\frac{\omega}{2}\right)

代入目标速度：

A_3(v)=4\sin^2\left(\frac{2\pi vT_{\mathrm{sameTX}}}{\lambda}\right)

当 v\rightarrow0 时：

A_2(v)=O(|v|)

而：

A_3(v)=O(v^2)

因此，三脉冲抵消器在零 Doppler 附近下降更快，可以进一步增强静态背景抑制。

但代价也很明确：极低速目标比使用 MTI2 时更容易受到明显衰减。

SVD 低秩杂波抑制

工程还提供基于整帧数据的 SVD 处理。

将 Range Cube 中的 Range-Channel 维展开：

\mathbf{X}\in\mathbb{C}^{(N_rN_{\mathrm{ch}})\times M}

进行奇异值分解：

\mathbf{X}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{H}

其中：

\sigma_1\ge\sigma_2\ge\cdots

如果环境中存在大面积稳定反射，不同 Range-Channel 单元的慢时间数据往往会形成共同的强相关分量。

这些成分可能集中在较大的奇异值及其对应子空间中。

取前 r 个奇异分量：

\mathbf{X}_r=\mathbf{U}_r\boldsymbol{\Sigma}_r\mathbf{V}_r^{H}

Eckart–Young 定理给出：

\mathbf{X}_r=\operatorname*{arg\,min}_{\operatorname{rank}(\mathbf{Z})\le r}\left\|\mathbf{X}-\mathbf{Z}\right\|_F

若将这一低秩部分作为稳定背景估计：

\mathbf{X}_{\mathrm{clean}}=\mathbf{X}-\mathbf{U}_r\boldsymbol{\Sigma}_r\mathbf{V}_r^{H}

也可以写成：

\mathbf{X}_{\mathrm{clean}}=\mathbf{X}\left(\mathbf{I}-\mathbf{V}_r\mathbf{V}_r^{H}\right)

前 r 个分量对应的能量比例：

\eta_r=\frac{\sum_{i=1}^{r}\sigma_i^2}{\sum_i\sigma_i^2}

如果 r 选择得过大，持续存在或者移动速度非常低的真实目标也可能进入主导子空间。

因此：SVD 并不是“删掉的主成分越多越好”，而是需要结合目标速度范围和实际环境选择合适的秩。

二维 CFAR 检测

完成杂波抑制和 Doppler FFT 以后，检测在 Range-Doppler功率图上进行：

P[r,d]=\left|X_D[r,d]\right|^2

固定门限很难同时适应：

• 近距离强反射；
• 远距离弱背景；
• 局部杂波边缘；
• 多个目标同时进入训练区域。

因此工程采用局部 training cells 估计 CUT 周围的背景功率，再根据目标虚警概率生成门限。

CA-CFAR

在均匀背景假设下，取 N 个 Training Cells：

X_1,X_2,\ldots,X_N

如果复高斯噪声经过平方检波，则单个功率样本服从指数分布。

p_X(x)=\frac{1}{\mu}e^{-x/\mu},\qquad x\ge 0

背景功率估计取均值：

Z=\frac{1}{N}\sum_{i=1}^{N}X_i

门限：

T=\alpha Z

在只有噪声的假设 H_0 下，可以得到：

P_{\mathrm{fa}} = P(X_{\mathrm{CUT}}>\alpha Z).

对于给定 Z=z：

P(X_{\mathrm{CUT}}>\alpha z\mid Z=z) = e^{-\alpha z/\mu}.

对 Z 求期望：

P_{\mathrm{fa}} = E \left[ e^{-\alpha Z/\mu} \right].

令

Y=\sum_{i=1}^{N}X_i,

则 Y 服从 Gamma 分布，且 Z=Y/N。利用 Gamma 分布的拉普拉斯变换：

P_{\mathrm{fa}}=\left(1+\frac{\alpha}{N}\right)^{-N}

因此 CA-CFAR 阈值倍率：

\alpha=N\left(P_{\mathrm{fa}}^{-1/N}-1\right)

最终检测条件：

P_{\mathrm{CUT}}>\alpha Z

也就是：

P_{\mathrm{CUT}}>N\left(P_{\mathrm{fa}}^{-1/N}-1\right)\frac{1}{N}\sum_{i=1}^{N}X_i

这样，每一个待检测位置的门限都会根据局部背景功率动态变化。

OS-CFAR

当 Training Cells 中混入其他强目标时，简单均值容易被高功率样本明显抬高。

OS-CFAR 会先将 training cells 从小到大排序：

X_{(1)}\le X_{(2)}\le\cdots\le X_{(N)}

选择第 k 个次序统计量：

Z_{\mathrm{OS}}=X_{(k)}

检测门限：

T_{\mathrm{OS}}=\alpha_{\mathrm{OS}}X_{(k)}

在指数背景条件下：

P_{\mathrm{fa}}=\prod_{i=0}^{k-1}\frac{N-i}{N-i+\alpha_{\mathrm{OS}}}

因此工程可以根据 N、k 和设定的 P_{\mathrm{fa}} 求解对应的 \alpha_{\mathrm{OS}}。

GOCA、SOCA 与 VI-CFAR

如果 CUT 位于不同背景区域的交界处，训练窗不同方向的平均功率可能存在明显差异。

将训练区域划分成若干个 sector，第 i 个 sector 的平均功率记为：

\mu_i=\frac{1}{N_i}\sum_{j=1}^{N_i}X_{i,j}

GOCA 使用较大的背景估计：

Z_{\mathrm{GOCA}}=\max_i\mu_i

SOCA 使用：

Z_{\mathrm{SOCA}}=\min_i\mu_i

二维 VI-CFAR

ATLAB ADC 全链路工程还实现了二维 variability-index 自适应策略。

对第 i 个 sector 计算均值和方差：：

\mu_i = \frac1{N_i} \sum_{j=1}^{N_i}X_{i,j},

s_i^2 = \operatorname{Var} \{X_{i,1},\ldots,X_{i,N_i}\}.

工程使用

VI_i=\frac{s_i^2}{\mu_i^2+\varepsilon}

其中 s_i^2 表示该区域的功率方差。

这一指标可以描述 sector 内功率起伏程度。该量相当于功率样本变异系数平方的形式。

若

VI_i\le VI_{\mathrm{th}},

当前 sector 被视为相对均匀。

再比较不同 sector 的平均功率：

\rho=\frac{\max_i\mu_i}{\min_i\mu_i+\varepsilon}

工程会结合 VI_i 与 \rho 判断当前背景状态：

• 背景相对均匀时，采用 CA 类背景估计；
• 不同方向均值差异明显时，更接近 clutter edge，可转向 GOCA；
• 只有部分 sector 满足均匀性要求时，仅使用相对可信的局部背景；
• 缺少稳定 sector 时，使用更加保守的 fallback。

最终仍然通过：

T_{\mathrm{VI}}=\alpha Z

与 CUT 进行比较。

简单来说：不是给整张 Range-Doppler 图设置一个固定门槛，而是在每一个待检测位置附近重新判断——这里的背景到底是什么水平。

二维亚栅格距离—速度估计

CFAR 首先获得的是整数 Range-Doppler Bin。

但真实目标并不一定刚好位于某一个 FFT 栅格中心。

如果目标峰值恰好落在两个相邻 Bin 之间，只取最大 FFT Bin 就会带来一定量化误差。

因此，本项目在检测结果附近进一步执行二维亚栅格距离—速度估计。

MATLAB ADC 全链路工程在检测点附近提取 3X3 对数功率 patch，并使用二维二次曲面拟合：

z(\Delta r,\Delta d)=a\Delta r^2+b\Delta d^2+c\Delta r\Delta d+d\Delta r+e\Delta d+f

其中：

• \Delta r：相对粗 Range Bin 的偏移；
• \Delta d：相对粗 Doppler Bin 的偏移。

\mathbf H = \begin{bmatrix} 2a & c\\ c & 2b \end{bmatrix},

线性项写为

\mathbf g = \begin{bmatrix} d\\ e \end{bmatrix}.

局部极值满足：

\mathbf{H}\boldsymbol{\delta}+\mathbf{g}=\mathbf{0}

因此：

\boldsymbol{\delta}^{\star}=-\mathbf{H}^{-1}\mathbf{g}

其中：

\boldsymbol\delta^\star = \begin{bmatrix} \Delta r^\star\\ \Delta d^\star \end{bmatrix}.

令

D=4ab-c^2,

则闭式结果为

\Delta r^\star = \frac{ce-2bd}{D}

\Delta d^\star = \frac{cd-2ae}{D} .

为了确保拟合结果对应局部峰值，而不是鞍点或谷底，工程还会检查 Hessian 的负定性，并限制亚栅格偏移不能距离粗检测 Bin 过远。

二维情况下可检查

a<0,

4ab-c^2>0.

另外工程代码还限制亚栅格偏移不能离粗检测格点过远，并在二维拟合不稳定时回退到一维抛物线插值。

对三个相邻对数功率 L,C,R：

\delta = \frac12 \frac{L-R}{L-2C+R} .

最终映射回物理坐标：

R_{\mathrm{refined}}=R_{k_r}+\Delta r^{\star}\Delta R_{\mathrm{bin}}

v_{\mathrm{refined}}=v_{k_d}+\Delta d^{\star}\Delta v_{\mathrm{bin}}

需要注意：亚栅格算法降低的是单个峰值的位置量化误差，并不会突破雷达带宽和 CPI 决定的真实物理分辨能力。

如果两个目标本身已经小于系统的可分辨尺度，仅靠插值无法把两个不可分目标重新分开。

阵列模型与角度估计

距离和速度主要来自时间维采样，而方位角来自不同天线通道之间的空间相位差。

设第 m 个有效阵元的位置为 x_m，目标方位角为 \theta。

在远场窄带近似下，第 m 个通道相对参考阵元的传播相位可以表示为：

\phi_m(\theta)=\frac{2\pi}{\lambda}x_m\sin\theta

对应 Steering Vector：

\mathbf{a}(\theta)=\begin{bmatrix}\exp\left(j\frac{2\pi}{\lambda}x_1\sin\theta\right)\\exp\left(j\frac{2\pi}{\lambda}x_2\sin\theta\right)\\vdots\\exp\left(j\frac{2\pi}{\lambda}x_M\sin\theta\right)\end{bmatrix}

如果阵元坐标已经以波长为单位：

p_m=\frac{x_m}{\lambda}

则等价写成：

a_m(\theta)=\exp\left(j2\pi p_m\sin\theta\right)

工程使用配置中的阵列位置建立 steering matrix，并对通道相位补偿。

如果第 m 个通道存在固定相位误差 \phi_m^{(\mathrm{cal})}，校准后数据为：

x_m^{(\mathrm{cal})}=x_m\exp\left(-j\phi_m^{(\mathrm{cal})}\right)

因此，角度估计结果不能只从算法本身理解。

阵列位置、通道幅相一致性以及标定质量，都会直接影响后续 DOA 方法。

Angle FFT

对于理想等间距线阵：

x_m=md\lambda

Steering Phase：

a_m(\theta)=\exp\left(j2\pi md\sin\theta\right)

定义空间频率：

u=d\sin\theta

阵列数据沿通道维执行 DFT 后，FFT 峰值对应空间频率 u。

进一步得到：

\theta=\arcsin\left(\frac{u}{d}\right)

本项目面对的 virtual array 并不简单假定为理想连续 ULA，因此 Angle FFT 还会结合实际阵元位置进行排序和空间栅格映射，再执行加窗和零填充 FFT。

DML

单目标窄带模型可以写成：

\mathbf{x}=\mathbf{a}(\theta)s+\mathbf{n}

其中 s 为目标复幅度。

DML 对每一个候选角度求解：

J(\theta,s)=\left|\mathbf{x}-\mathbf{a}(\theta)s\right|_2^2

固定 \theta 后，对 s 求最小二乘解：

\hat{s}(\theta)=\frac{\mathbf{a}^{H}(\theta)\mathbf{x}}{\mathbf{a}^{H}(\theta)\mathbf{a}(\theta)}

代回目标函数：

J(\theta)=\left|\mathbf{x}\right|_2^2-\frac{\left|\mathbf{a}^{H}(\theta)\mathbf{x}\right|^2}{\mathbf{a}^{H}(\theta)\mathbf{a}(\theta)}

如果 Steering Vector 已归一化满足：

\mathbf{a}^{H}(\theta)\mathbf{a}(\theta)=1

则：

J(\theta)=\left|\mathbf{x}\right|_2^2-\left|\mathbf{a}^{H}(\theta)\mathbf{x}\right|^2

工程使用残差倒数形成扫描谱：

P_{\mathrm{DML}}(\theta)=\frac{1}{J(\theta)+\varepsilon}

最大谱峰对应当前单源模型下的角度估计。

MUSIC

MUSIC 不直接使用单个通道向量，而是从多个局部 Snapshot 中估计协方差。

设：

\mathbf{X}=\begin{bmatrix}\mathbf{x}_1&\mathbf{x}_2&\cdots&\mathbf{x}_L\end{bmatrix}\in\mathbb{C}^{M\times L}

样本协方差矩阵：

\widehat{\mathbf{R}}=\frac{1}{L}\mathbf{X}\mathbf{X}^{H}

在有限 Snapshot 条件下，加入 diagonal loading：

\mathbf{R}_L=\widehat{\mathbf{R}}+\gamma\frac{\operatorname{tr}\left(\widehat{\mathbf{R}}\right)}{M}\mathbf{I}

对其进行特征分解：

\mathbf{R}_L=\mathbf{E}\boldsymbol{\Lambda}\mathbf{E}^{H}

如果采用单源模型，将最大特征值对应方向归入信号子空间，其余特征向量组成噪声子空间 \mathbf{E}_n。

理想情况下，真实目标对应的 Steering Vector 与噪声子空间近似正交：

\mathbf{E}_n^{H}\mathbf{a}(\theta_0)\approx\mathbf{0}

因此 MUSIC 伪谱：

P_{\mathrm{MUSIC}}(\theta)=\frac{1}{\mathbf{a}^{H}(\theta)\mathbf{E}_n\mathbf{E}_n^{H}\mathbf{a}(\theta)+\varepsilon}

在阵列模型准确、Snapshot 足够，并且信号子空间与噪声子空间能够有效区分的条件下，真实方向附近会形成明显谱峰。

有限通道、相关多径以及阵列误差都会改变最终结果。

OMP 稀疏角度恢复

将待搜索角度离散为：

\theta_1,\theta_2,\ldots,\theta_G

构造 Steering Dictionary：

\mathbf{A}=\begin{bmatrix}\mathbf{a}(\theta_1)&\mathbf{a}(\theta_2)&\cdots&\mathbf{a}(\theta_G)\end{bmatrix}

角域稀疏模型：

\mathbf{x}=\mathbf{A}\mathbf{s}+\mathbf{n}

其中只有少数角度单元具有非零系数。

理想稀疏恢复可以表示为：

\min_{\mathbf{s}}\left|\mathbf{s}\right|_0

同时满足：

\left|\mathbf{x}-\mathbf{A}\mathbf{s}\right|_2\le\epsilon

OMP 通过逐次选择与当前残差最相关的原子近似求解。

初始化：

\mathbf{r}^{(0)}=\mathbf{x}

\mathcal{S}^{(0)}=\varnothing

第 t 次迭代选择：

k^{\star}=\operatorname*{arg,max}_{k}\left|\mathbf{a}_k^{H}\mathbf{r}^{(t-1)}\right|

更新支持集：

\mathcal S^{(t)} = \mathcal S^{(t-1)} \cup \{k^\star\}.

对已经选择的原子重新进行最小二乘：

\widehat{\mathbf{s}}*{\mathcal{S}^{(t)}}=\operatorname*{arg,min}*{\mathbf{z}}\left|\mathbf{x}-\mathbf{A}_{\mathcal{S}^{(t)}}\mathbf{z}\right|_2^2

残差更新：

\mathbf{r}^{(t)}=\mathbf{x}-\mathbf{A}*{\mathcal{S}^{(t)}}\widehat{\mathbf{s}}*{\mathcal{S}^{(t)}}

当前工程主要将 OMP 用于单主源角度估计，因此它更适合作为 Angle FFT、DML、MUSIC 之外的一种稀疏角度恢复参考。

MIMO、TDM 与 DDMA

在多发多收系统中，最终角度处理面对的不只是物理 RX 通道，还包括不同 TX 与 RX 组合形成的 Virtual Channel。

设第 t 个发射阵元位置为：

x_t^{(\mathrm{TX})}

第 r 个接收阵元位置为：

x_r^{(\mathrm{RX})}

在远场条件下，虚拟阵元位置可以表示为：

x_{t,r}^{(\mathrm{virt})}=x_t^{(\mathrm{TX})}+x_r^{(\mathrm{RX})}

因此，一个 N_{\mathrm{TX}}\times N_{\mathrm{RX}} 系统理论上最多可以形成：

N_{\mathrm{TX}}N_{\mathrm{RX}}

个 TX-RX 组合。

但前提是：不同 TX 的回波能够被正确分离，并且保持所需要的相位关系。

TDM 对速度轴的影响

如果 G 个 TX 轮流发射，那么同一个 TX 每隔 G 个 chirp 才再次出现。

它自己的慢时间采样周期变为：

T_{\mathrm{sameTX}}=GT_c

因此，该 TX 对应的无模糊速度范围：

v_{\max}=\frac{\lambda}{4GT_c}

速度 FFT 栅格也应该按照：

\Delta v_{\mathrm{bin}}=\frac{\lambda}{2N_dGT_c}

计算。

如果仍然错误地使用单个 Chirp Period T_c，速度坐标就会产生比例误差。

DDMA 相位编码

DDMA 通过跨 chirp 的相位编码，将不同 TX 映射到不同的慢时间频率位置。

设第 t 个 TX 的相位编码：

\phi_t[m]=\phi_{0,t}+m\Delta\phi_t

其编码因子：

c_t[m]=\exp\left(j\phi_t[m]\right)

每个 chirp 增加 \Delta\phi_t，对应慢时间编码频率：

f_{\mathrm{code},t}=\frac{\Delta\phi_t}{2\pi T_c}

在 N_d 点 Doppler FFT 中，对应 Bin 偏移：

k_{\mathrm{offset},t}=\frac{\Delta\phi_t}{2\pi}N_d

如果物理目标 Doppler 为 f_D，编码后的慢时间频率位于：

f_{\mathrm{obs},t}=f_D+f_{\mathrm{code},t}

已知相位编码后，可以通过乘以共轭编码进行解调：

x_t^{(\mathrm{demod})}[m]=x_t[m]\exp\left[-j\left(\phi_{0,t}+m\Delta\phi_t\right)\right]

本项目保持了一个非常重要的工程原则：只有在配置中能够明确得到 TX Group 与相位编码关系时，才进行对应的通道组织和 DDMA 解调。

对于缺少明确 DDMA Phase Increment、Doppler Offset 或 Virtual Channel 信息的数据，不根据频谱外观主观反推编码参数。

避免把编码偏移错误解释为真实目标速度。

从距离、速度到目标点云

完成亚栅格距离估计和 DOA 后，每一个检测目标可以表示为：

\left(\widehat{R},\widehat{v},\widehat{\theta}\right)

在二维水平面上：

x=\widehat{R}\sin\widehat{\theta}

y=\widehat{R}\cos\widehat{\theta}

工程进一步使用径向速度作为第三维：

z=\widehat{v}

因此，项目图中的三维坐标实际表示：

\left(x,y,v_r\right)

而不是包含高度信息的 (x,y,z) 三维空间位置。

所以在理解最终点云时需要特别注意：第三维是径向速度，不是目标高度。

代码与数学模型的对齐

derive_radar_parameters.m

代码首先计算

\lambda=\frac{c}{f_c}, \qquad f_s=\frac{f_{\mathrm{ADC}}}{D}, \qquad S=\frac{B}{T_r},

其中 D 为 ADC decimation factor。

Range FFT 的坐标间隔按

\Delta R_{\mathrm{grid}} = \frac{cf_s}{2SN_r}

建立，因此代码变量 range_resolution_m 更准确的物理含义是 range FFT grid spacing。

经典波形分辨率则是

\Delta R_{\mathrm{waveform}} \approx \frac{c}{2B}.

在 Pf0 中，两者分别约为 0.5245\ \mathrm m 与 0.4997\ \mathrm m，数值接近，但来源不同：前者由采样率、斜率与 FFT 点数共同决定，后者由有效调频带宽决定。

对于 Doppler，代码并不固定使用单个 chirp period，而是先解析 TX group。

若相同 TX 每隔 G 个 chirp 才再次出现，则

T_{\mathrm{sameTX}}=GT_c,

于是

\Delta v_{\mathrm{grid}} = \frac{\lambda}{2N_dT_{\mathrm{sameTX}}}, \qquad v_{\max} = \frac{\lambda}{4T_{\mathrm{sameTX}}}.

这一步决定了 TDM 情况下速度轴能否被正确解释。

cfar_2d.m

工程支持 CA、OS、GOCA、SOCA 和 VI 五类检测路径。

• CA：用全部有效训练单元的平均功率估计背景；
• OS：先排序，再取第 k 个 order statistic；代码根据指数背景下的 P_{\mathrm{fa}} 方程数值求解倍率；
• GOCA / SOCA：把二维训练环划分为 sector，分别取 sector 均值的最大值 / 最小值；
• VI：先计算每个 sector 的 VI=s^2/(\mu^2+\varepsilon)，再根据 sector 内是否均匀以及 sector 间均值比选择 CA、GOCA、可信 sector 平均或保守 fallback。

因此 VI-CFAR 的“自适应”发生在 背景估计器选择 上，最终仍然形成

T=\alpha Z

并与 CUT 比较。

angle_fft_doa.m：不把非均匀阵列直接当 ULA

工程先按真实阵元位置排序，并映射到指定空间网格；重复位置做合并，再加窗和零填充 FFT，最终通过空间频率

u=d\sin\theta

映射到角度。

DML、MUSIC 与 OMP

DML 使用归一化 steering dictionary。

对单目标模型，扫描角度的投影能量为

|\mathbf a^H(\theta)\mathbf x|^2,

残差为

J(\theta)=\|\mathbf x\|_2^2-|\mathbf a^H(\theta)\mathbf x|^2,

代码用 1/(J+\varepsilon) 形成谱。

MUSIC 使用多个 snapshot 构造样本协方差，并进行 trace-scaled diagonal loading，再分离信号子空间和噪声子空间。它依赖阵列模型、快拍数与通道校准，谱峰更尖不等于物理孔径变大。

OMP 逐次选择与当前残差最相关的 steering atom，再对当前支持集做最小二乘更新。当前工程主要把它作为单主源角度恢复参考；用于统一绘图的参考相关谱不应被误解为“最终稀疏系数的连续功率谱”。

decode_tx_groups.m：只在编码参数明确时处理 DDMA

代码从 tx_groups 解码真实 TX group。对随附配置，Pf0/Pf2 为 TX3 单发，Pf1/Pf3 为 TX1 单发，不把这些配置误判为 DDMA。

只有配置明确给出 TX DDMA phase increment 或等价 Doppler offset 时，才建立

f_{\mathrm{code},t} = \frac{\Delta\phi_t}{2\pi T_c}, \qquad k_{\mathrm{offset},t} = \frac{\Delta\phi_t}{2\pi}N_d

并进行解调。

若同时发射但缺少可验证的编码信息，代码选择报错，而不是根据频谱外观反推物理编码。

从算法到工程结果

原始 ADC 解包和距离维处理：

二维距离—速度结果：

慢时间杂波抑制对比：

二维 CFAR 与 VI-CFAR 诊断：

整数栅格检测之后的二维亚栅格细化：
![51663a2307068f6b4d2835ea11009732.jpg](https://super-radar-pro.oss-cn-shenzhen.aliyuncs.com/file/jpg/2026-08-25/51663a2307068f6b4d2835ea11009732_2f5582cf.jpg）

四类 DOA 方法：

最终目标点云：

这里的三维点云使用

(x,y,z) = (R\sin\theta,\ R\cos\theta,\ v_r),

第三维是径向速度，不是高度。

工程复现

项目工程位于：https://github.com/Super-Radar/CTSAI-A100/tree/main/community/projects/a100-matlab-adc-signal-processing

进入：

community/projects/a100-matlab-adc-signal-processing/

使用 MATLAB 运行：

main

配置测试：

addpath(fullfile(pwd, 'tests'));
run_configuration_tests;

主要输出写入：

results/

包括 PNG、CSV 和 MAT 文件。

项目同时提供完整技术说明文档，覆盖：

• ADC 数据格式；
• LFMCW 信号模型；
• Range FFT；
• Doppler FFT；
• 杂波抑制；
• CFAR；
• 亚栅格细化；
• 阵列模型；
• Angle FFT；
• DML；
• MUSIC；
• OMP；
• MIMO / TDM / DDMA。

如果希望继续深入具体数学推导或者修改算法，可以直接进入项目源码与技术文档进一步研究。

当前边界

作为一套面向算法学习、验证和二次开发的工程，当前结果仍需要在明确边界下理解：

角度结果依赖真实阵列位置、固定相位误差、通道幅相一致性和有限 Snapshot；

MUSIC、DML 或 OMP 获得更加尖锐的谱峰，并不意味着雷达物理孔径被算法“放大”；

MTI 会衰减极低速目标；

SVD 删除过多主奇异分量，也可能将稳定目标或极低速目标作为背景削弱；

MIMO / DDMA 的物理解调必须建立在明确的 TX Chirp Schedule、Phase Increment 或等价 Doppler Offset 上；

亚栅格插值能够改善单峰的位置量化误差，但不能突破带宽与 CPI 决定的真实分辨能力；

当前工程定位是可复现的算法学习与二次开发基线，不将单组数据结果包装为 CTSAI A100 的产品级性能承诺。