一统江湖：毫米波雷达开发手册之大话线谱估计

严肃声明

​ 严格禁止未经本人允许转载本文工作成果（包括引言、证明、代码分析等），一经发现被用于学术或商业用途，将予以法律警告。欢迎和笔者深度合作，探讨学术话题。
本人发现部分网站（如https://it.cha138.com/等野鸡网站）未经本人允许私自将本人博客内容转发公开，本人严禁杜绝此类行为，一经发现将投诉举报相关行为并追究违法责任！！！

写在前面

​ 深知新手在接触毫米波雷达板硬件时需要花费的沉没成本，因此在行将告别毫米波雷达之际，总结这两年以来在毫米波雷达上的一些经验和教训。

​ 本文档用于为实现基于AWR1243BOOST等单板毫米波雷达开发提供参考指南与解决方案，主要包括硬件配置、基础参数、信号模型、应用DEMO开发以及可深入研究方向思考等；为更好地匹配后续级联雷达应用的学习路线，在本手册中会尽可能同化单板雷达和级联雷达中的相关表述。

​ 本指南作者信息：Xl，联系方式：xxxxx@zju.edu.cn。未经本人允许，请勿用于商业和学术用途。

​ 希望后者在使用本指南时可以考虑引用作者在毫米波雷达旅途中的相关工作，如本文参考文献[1].
本章节为可深入研究方向思考章节之线谱估计，主要讨论线谱估计有关的几个问题。
欢迎各位读者通过邮件形式与笔者交流讨论，本章节完整程序请私信笔者，希望使用本代码时能够提供一份引用和Star，以表示对笔者工作的尊重，谢谢！在后续将定时维护更新。
https://github.com/DingdongD/TDMA-MIMO

往期内容：
炉火纯青：毫米波雷达开发手册之大话空间谱估计
登堂入室：毫米波雷达开发手册之信号模型
初出茅庐：毫米波雷达开发手册之基础参数
扬帆起航：毫米波雷达开发手册之硬件配置
眼观四海：自动驾驶&4D成像毫米波雷达 如今几何？

讨论概述

在正式开始本文讨论前，不妨思考下列问题：
频率估计是什么？为什么要频率估计？
空间谱估计的方法是否可以用来做频率估计？
如果可以，是否直接套用到频率信号上即可？

讨论1

频率估计通常通常可以分为非参数化估计和参数化估计，在参数化估计中可以细分为连续谱估计和线谱估计。连续谱估计的常用方法是ARMA模型，线谱估计则是认为目标信号在某个域上可以用有限的原子/基信号去表征，其思想与子空间方法类似（但并不局限于子空间方法）。
在下面的讨论中，我们讨论的是线谱估计，线谱估计在毫米波雷达的距离、多普勒速度、空间信号估计以及呼吸、心跳信号检测中是贯穿全线的，如果要有一个词来一统雷达信号处理，笔者认为线谱估计能胜其任。
说到这里，频率估计的意义想必很明确了，目的就是为了从含噪信号中能够将信号映射到频域、空域分别解析出在该维度下的信息，总的来说可以如下总结（如有不合理，欢迎指正）：
快时间维度---映射至频域---距离信息
慢时间维度---映射至频域---多普勒速度信息
相位维度 ---映射至频域---呼吸/心跳信息
阵列维度 ---映射至空域---方位/俯仰信息

讨论2

空间谱估计方法能否直接用到频域上？不可以，但也可以。
很多人可能像@秋姐姐一样，自以为自己写的空间谱估计等价于频域信号估计，这里要注意的区别是：
1、空间谱估计，也称为波达角（DOA，Direction-of-Arrival）估计，这一估计通常是有角度频率限制的，例如，常用的限制区间为[0:180°]或[0:pi]；但是频率信号估计通常是没有约束条件的，我们可以对整个频带内的型号进行估计；
2、空间谱估计是在空域信号上的估计，空域信号对应的是阵列流形空间的表征模式，其通常可以分为单快拍和多快拍下信号，在单快拍下信号维度为\(\mathbb R^{M\times 1}\)，在多快拍下信号维度为\(\mathbb R^{M\times N}\)，在这里\(M\)表示阵元数目，\(N\)表示快拍数目（这里的快拍通常是能够表征同一距离采样点下的统计分布的累积，通常在毫米波雷达中对应慢时间维度的快拍）。而我们的频域信号维度通常是\(\mathbb R^{L}\)，\(L\)表示快拍数目（这里的快拍可以是快时间维度或慢时间维度的快拍，取决于需要恢复的信息所在维度）。
OK，不难发现：空间谱估计和频域信号估计的信号维度实际上是不统一的，前者为单快拍矢量或多快拍矩阵，而后者为矢量。但这并不影响，通常频域信号估计是对某个确定维度（距离、多普勒等）的估计，而我们所得到的数据不出意外为时域信号，时域信号实质上是对该维度信息的累积观测/统计表征，因此我们可以通过对这个矢量信号重构为类似于空域信号的模式表征，那么我们就可以基于新的模式表征实现空间谱估计（此空间谱估计不是真的空间谱估计，而是使用了空间谱估计方法追求了超分辨效果罢了）

那么，要怎么构建这个新的类似空间谱估计的模式表征呢？其实这个思想和子空间方法中的旋转子空间方法（ESPRIT）也是类似的，通过将大矩阵滑动分解为若干小矩阵来实现旋转不变性；这里也是通过滑动窗口模拟空间维度上的统计积累。

M = 16; % 虚拟阵元数目，这个可以抽象为空间谱估计中的阵元
N = length(data); % data为你的时域输入信号
    for i = 1 : N - M
        X(:, i) = data(i + M - 1 : -1 : i).'; % 构建新的模式表征
    end

讨论3

下面，我们希望给出一份实际案例来更好地帮助各位理解。
下面的程序案例中给出了实信号和复信号在FFT/MUSIC算法下的频率谱测试，如果想尝试其他空间谱估计的超分辨算法（旋转子空间方法/压缩感知OMP、L1NORM、IAA、OFFGRID、ANM等），可以阅读笔者的炉火纯青：毫米波雷达开发手册之大话空间谱估计
到这里为此，我们某种程度上已经统一了毫米波雷达信号处理的思想，并且实际上也已经为各种应用（距离维信号重构、多普勒维信号重构、生命体征重构等）开篇。

clc;clear;close all;
%% 本文件用于实现多频率信号的估计
%% By Xuliang

fs = 40e3; % 采样频率
Tc = 4e-3; % 持续时间
n = 1 / fs : 1 / fs : Tc; %采样点
Npoint = 2^nextpow2(length(n));

snr = 10; % 信噪比
% sig = cos(2*pi*3259*n + pi /3); % 单目标
% sig = 2 * cos(2 * pi * 2182 * n + pi * 4 /3) + 1.6 * cos(2 * pi * 2550 * n+ pi / 8) + 1 * cos(2 * pi * 4100 * n + pi / 4); % 多目标

sig = 3 * exp(1j * 2 * pi * 2182 * n) + 1.5 * exp(1j * 2 * pi * 2550 * n) + exp(1j * 2 * pi * 4250 * n);
sig = awgn(sig, snr, 'measured'); % 加高斯白噪声

fftdata = (fft(sig, Npoint)); % FFT
freq_idx = fs / Npoint : fs / Npoint : fs / 2; % 频率间隔
figure(1);
plot(freq_idx, abs(fftdata(1:end/2))); % FFT可视化
title('FFT谱');

M = 32; % 信号子空间维度，可以通过修改这个值来检查算法效果
P = 3; % 信源数目 这里要考虑目标数目*2

if isreal(sig)
    P = P * 2; % 如果为实信号 需要考虑信源*2 这个可以去结合特征值进一步理解
else
    P = P; % 复信号 信源数目不变
end

music_spectrum = FREQ_MUSIC(sig, M, P, Npoint, fs); % MUSIC估计
figure(2);
plot(freq_idx, db(music_spectrum(1:Npoint/2))); 
title('MUSIC谱');

function [PoutMusic] = FREQ_MUSIC(data, M, P, Npoint, fs)
    % By Xuliang
	% data: 输入信号(呼吸信号，快时间信号，慢时间信号均可) 一维 N * 1
    % P: 信源个数
    % M: 自相关矩阵阶数
    % Npoint: 频率点数
	% fs : 采样频率
    
    searchGrids = fs / Npoint : fs / Npoint : fs; % 频率间隔
    N = length(data); % 输入信号数据长度
    for i = 1 : N - M
        X(:, i) = data(i + M - 1 : -1 : i).'; % 构建新的模式表征
    end
    
    snap = N - M; % 快拍数
    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]' * searchGrids(id) / fs); % 导向矢量
        PoutMusic(id) = (abs(1 / (atheta_vec' * EN * EN' * atheta_vec))) ; % 功率谱计算
    end
end

讨论4

突然间多了一个讨论4，也不知道是为什么？@秋同学在5.31日这一天成功地引起了“咱也不知道说什么”的新一轮烧脑。事情是什么个事情了？故事如下：
小秋买了一块毫米波雷达，并根据说明书采集了数据格式为256×100×86×1的原始雷达数据，256表示快时间维度，100表示慢时间维度，86表示空间维度（阵元），小秋对着手册把原始信号按传统的距离多普勒FFT-CFAR处理得到了每个检测目标对应的空域信号\(X\in \mathbb R^{86\times 1}\),然后小秋又用信号X进行了DOA估计，哎，那么很好奇且很刚的小秋就来了个问题：为什么我不直接用满足\(R^{86 \times 25600}\)的数据做成像？小秋坚决觉得这样做是可以的，小鱼表示不同意他的观点。两人的观点如下：
小秋：可以直接用\(R^{86 \times 25600}\)的数据来做DOA估计；
小鱼：不建议使用\(R^{86 \times 25600}\)的数据来做DOA估计。
两人就这样挣扎了一下午...
最后小鱼来这里贴了一下小鱼的看法。
DOA估计一般是对确定距离单元下的空域信号分析，其通常为关于阵列的矢量或矩阵信号，矢量体现为单快拍，矩阵体现为长时观测信号的积累即多快拍，这里的快拍一般指慢时间维度。小秋的想法是希望利用快时间维度的快拍信息来进行DOA估计，存在的问题在于：每个距离单元可能对应多个目标，他们的空间分布统计特性不一致，这种长时观测积累不利于一致性判断；此外，这样处理可能会引入过多的目标，超出阵列可检测的最大目标容限，以现有的DOA估计算法很难处理过多目标的空域信号。

参考文献

[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.