写在前面
本文档用于为实现基于AWR1243BOOST等单板毫米波雷达开发提供参考指南与解决方案，主要包括硬件配置、基础参数、信号模型、应用DEMO开发以及可深入研究方向思考等；为更好地匹配后续级联雷达应用的学习路线，在本手册中会尽可能同化单板雷达和级联雷达中的相关表述。
Xuliang，联系方式：22134033@zju.edu.cn。未经本人允许，请勿用于商业和学术用途。
本章节为可深入研究方向思考章节之空间谱估计，主要解读子空间方法和压缩感知方法。
欢迎各位读者通过邮件形式与笔者交流讨论，本章节完整程序请私信笔者，希望使用本代码时能够提供一份引用和Star，以表示对笔者工作的尊重，谢谢！在后续将定时维护更新。
往期内容：
登堂入室：毫米波雷达开发手册之信号模型
眼观四海：自动驾驶&4D成像毫米波雷达 如今几何？
扬帆起航：毫米波雷达开发手册之硬件配置
初出茅庐：毫米波雷达开发手册之基础参数
空间谱估计算法
信号模型
目标空间通常是一个由信号源的参数与复杂环境参数张成的空间。观察空间是利用空间按一定方式排列的阵元来接收目标空间的辐射信号，接收数据中往往包含信号特征（方位、距离、极化等）和空间环境特征（噪声、杂波、干扰等），观察空间是一个多维空间。估计空间是利用空间谱估计技术从复杂的观察数据中提取信号的特征参数。对M阵列而言，基于窄带信号的假设，可以认为\(\frac{1}{c}\boldsymbol{r}^{T}\boldsymbol{a}\ll\frac{1}{B},\因此阵列信号以向量表示为\(\boldsymbol{s(t}=s(t[e^{-jr_1^Tk},e^{-jr_2^Tk},....,e^{-jr_M^Tk}]\，,若设第一个阵元为基准点且初始相位为0，那么可以得到导向矢量（方向矢量，\(M×1\维）为:\(\boldsymbol a(\boldsymbol\theta=\left[\boldsymbol1,,\boldsymbol e^{-jr_2^Tk},....,\boldsymbol e^{-jr_M^Tk}\right]^T=\left[\boldsymbol1,,\boldsymbol e^{jr_2^Tk},....,\boldsymbol e^{jr_M^Tk}\right]^H.\
M阵列-K信号源模型而言，
M个具有全向性的阵元按任意排列构成，并设有K个具有相同中心频率\(\omega_0\、波长为\(\lambda\的空间窄带平面波分别以角\(\Phi_1,\Phi_2,\dots\Phi_K\入射，\(\Phi_{i}=(\theta_{i},\phi_{i}\，阵列第\(m\个阵元的输出可以表示为：\(x_m(t=\Sigma_{i=1}^K s_i(te^{j\omega_0\tau_m(\Phi_i} + n_m(t\AWR1243BOOST等毫米波雷达中需要分析方位图、点云图，方位角和俯仰角是实现这些分析的基础。
\(X(t=[x_1(t,x_2(t,...,x_M(t]^T\
\(S(t=[s_1(t,s_2(t,...,s_M(t]^T\
观察空间可用噪声空间、目标空间和流形表示，即满足下式：\(\boldsymbol{X}(\boldsymbol{t}=\boldsymbol{A}(\boldsymbol{\Phi}S(\boldsymbol{t}+N(\boldsymbol{t}\overline{}\
名词解释
波数向量：\(|\boldsymbol{k}|=\frac{\omega}{c}=\frac{2\pi}{\lambda}\，空间距离变化
1m时相位的变化量延时时间：\(\frac{1}{c}\boldsymbol{r}^T\boldsymbol{a}\
滞后相位：\(r^Tk\
方向矢量/导向矢量：阵列相对基准阵列的相位向量\(a(\Phi\
阵列流形：流形表示由
K个信号源表示的方向向量组成的矩阵\(A(\Phi\快时间维&慢时间维：脉冲雷达往往伴随快时间和慢时间维这两个概念，毫米波
FMCW雷达本质上也是一种脉冲雷达。快时间维是针对单个脉冲而言，是对接收回波信号按行存储。快时间维的采样频率为采样率。慢时间维是针对多个脉冲而言，是对接收回波信号按列存储。每M个脉冲参与处理，脉冲间的时间间隔为脉冲重复间隔（Pulse Repetition Interval，PRI），脉冲重复频率（Pulse Repetition Frequency，PRF）是1/PRI。慢时间维的采样频率为PRF。暴力来说，在毫米波雷达里面快时间维可以理解为距离维，慢时间维可以理解为多普勒或速度维。毫米波雷达与空间信号关系
AWR1243BOOST等雷达采集数据通常由
DCA1000回传，回传数据文件需要经过解析得到雷达信号矩阵，该矩阵的维度为4D，分别为距离维ADC采样点数×虚拟天线通道数×Chirp Loop数目（单帧发送的Chirp数目）×帧数，但是我们这里研究的空间信号是2D，这要如何匹配呢？4D tensor数据实际在处理的时候需要逐帧处理，每帧的
3d tensor数据经过快时间维FFT和慢时间维FFT（或称2D-FFT）得到距离多普勒谱图，再在距离多普勒谱图上做恒虚警率检测得到对应距离和速度的目标信号，此时我们得到的目标信号维度是关于天线维度的向量，这个向量可以理解为单快拍下的空间信号矩阵\(X\；实际上，我们也可以通过仅对快时间FFT后的数据作恒虚警率检测得到特定距离的目标，这个时候我们得到的目标信号维度是和天线维度、Chirp数目相关的矩阵，那这个矩阵可以理解为多快拍下的空间信号矩阵\(X\，这里谈到的\(X\是和信号模型中提到的\(X\等价的。\(\boldsymbol{R}=\boldsymbol{E}(\bigl(X(\boldsymbol{t}-\boldsymbol{m}_x(\boldsymbol{t}\bigr\bigl(X(\boldsymbol{t}-\boldsymbol{m}_x(\boldsymbol{t}\bigr^H\
\(\textbf{R}=E\{\boldsymbol{X}(\boldsymbol{t}\boldsymbol{X}(\boldsymbol{t}^{\boldsymbol{H}}\}=E\left\{\big(\boldsymbol{A}(\boldsymbol{\Phi}\boldsymbol{S}(\boldsymbol{t}+\boldsymbol{N}(\boldsymbol{t}\big\big(\boldsymbol{A}(\boldsymbol{\Phi}\boldsymbol{S}(\boldsymbol{t}+\boldsymbol{N}(\boldsymbol{t}\big^{\boldsymbol{H}}\right\}=A(\PhiS(\boldsymbol{t}\boldsymbol{S}(\boldsymbol{t}^H A(\boldsymbol{\Phi}^H+\boldsymbol{\sigma}^2I=A(\boldsymbol{\Phi}\boldsymbol{R}_s A(\boldsymbol{\Phi}^H+\boldsymbol{R}_N\
目标子空间与噪声子空间的加形式，也可以通过特征分解表示为\(M\个特征值与特征矢量的积和形式。
子空间方法
Vanilla-MUSIC
这种将阵列信号的协方差矩阵进行特征分解，得到与信号分量对应的信号子空间和信号分量正交的噪声子空间，利用两个子空间的正交性来估计信号参数。根据信号子空间与噪声子空间的正交性，可以得到以下表述：
\(A(\PhiR_{s}A(\Phi^{H}U_{N}=0\
因为\(\boldsymbol{A\neq0},\boldsymbol{R_{s}\neq0},\故\(\boldsymbol{A(\Phi}^{H}\boldsymbol{U}_{N}=\boldsymbol{0}\，这表明在无噪声情况下的信号\(x(t=\Sigma_{i=1}^M s_i(ta_i(\theta\to x(t\in span\{a_1,a_2,...,a_K\}.\协方差矩阵大特征值对应的特征向量张成的空间与入射信号的导向矢量张成的空间是同一个空间。即\(\begin{matrix}U_{S}=[e_{1},e_{2},...,e_{K}],U_{N}=[e_{K+1},e_{K+2},...,e_{N}]\\ \textit{span}\{e_{1},e_{2},...,e_{K}\}=\textit{span}\{a_{1},a_{2},...,a_{K}\}\\ \end{matrix}\
根据信号子空间的导向矢量与噪声子空间正交的原理，信号入射方向上会出现极小值，因此空间谱函数可以表示为：\(\mathbf{P}_{\text{music}}(\theta=\frac{1}{\boldsymbol{a}^H(\PhiU_N U_N^H\boldsymbol{a}(\Phi}\
下面给出
Vanilla-MUSIC的单独求解方位角和联合求解角案例。function [PoutMusic] = DOA_MUSIC(X, P, searchGrids % By Xuliang % X: 输入信号 Channel * ChirpNum % P: 目标数目 % PoutMusic: 输出功率谱 M = size(X, 1; % 阵元数 snap = size(X, 2; % 快拍数 RX = X * X' / snap; % 协方差矩阵 [V, D] = eig(RX; % 特征值分解 eig_value = real(diag(D; % 提取特征值 [B, I] = sort(eig_value, 'descend'; % 排序特征值 EN = V(:, I(P+1:end; % 提取噪声子空间 PoutMusic = zeros(1, length(searchGrids; for id = 1 : length(searchGrids atheta_vec = exp(-1j * 2 * pi * [0:M-1]' * 1 / 2 * sind(searchGrids(id; % 导向矢量 PoutMusic(id = (abs(1 / (atheta_vec' * EN * EN' * atheta_vec ; % 功率谱计算 end end function [index1,index2,Pmusic] = MUSIC_2D(X,P %% MUSIC算法：用于方位角和俯仰角的联合估计 % By Xuliang % X: 输入信号 Channel * ChirpNum % P: 目标数目 M = size(X, 1; % 阵元数 snap = size(X, 2; % 快拍数 d2rad = pi / 180; % pi rad = 180 ° 1°= pi/180 rad lambda = physconst('lightspeed' / 77e9; % 波长 D = 0.5 * lambda; % 阵元间距 D1 = 0:D:(M-1*D; % TX0和TX2的线阵间距 D2 = 2 * D : D : 5 * D; % TX1的线阵间距 D_AZI = [D1,D2]; % TX0 TX2 和 TX1的水平间距 D_ELE = [zeros(1,8 D*ones(1,4]; % TX0/TX2和TX1的纵向间距 %% 计算协方差矩阵 RX = X1 * X1' / snap; [EV,D] = eig(RX; % 特征值分解 EVA = diag(D;% 特征值对角线提取并转为一行 [~,I] = sort(EVA;% 特征值排序从小到大，eig默认有排序 EV = fliplr(EV(:,I;% 特征矢量排序 EN = EV(:,P+1:M; % 噪声子空间 %% 遍历每个方位角和俯仰角 计算空间谱 for ele = 1:181 % 俯仰角遍历 phim(ele = ele - 91; phim1 = d2rad * phim(ele; for azi = 1:181 % 方位角遍历 theta(azi = azi - 91; theta1 = d2rad * theta(azi; a = exp(-1j * 2 * pi * (D_AZI * cos(phim1 * sin(theta1 + D_ELE * sin(phim1 / lambda.'; % 构建导向矢量 Pmusic(azi,ele = 1 / (a' * EN *EN' * a; % 空间谱函数 end end Pmusic = abs(Pmusic; [index1,index2] = find(Pmusic == max(max(Pmusic; % 找到空间谱谱峰并返回俯仰角和方位角 end% 如何对特征值可视化判断信源数？ [V, D] = eig(RX; % 特征值分解 SP = V(:, M-P+1:M; EN = V(:, 1:M-P; figure(1; % 特征值分布的可视化 plot(diag(D,'kd-'; xlabel('Number of Eigenvalues'; ylabel('Eigenvalues'; % 如何验证信号和噪声空间的正交性？ s_eigen_fft = fft(SP; n_eigen_fft = fft(EN; plot(abs(s_eigen_fft(:,1:2, 'ks-','LineWidth',1.5; hold on; plot(abs(n_eigen_fft(:,1:2, 'rd-','LineWidth',1.5; % legend('Signal Space', 'Noise Space';ESPRIT
Waiting .。
DML
压缩感知
网格模型
无网格模型
参考文献
[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。