526互联

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

发布时间 2023-05-28 11:38:18作者: 信海

写在前面

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

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

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

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

往期内容：
炉火纯青：毫米波雷达开发手册之大话空间谱估计
登堂入室：毫米波雷达开发手册之信号模型
初出茅庐：毫米波雷达开发手册之基础参数
扬帆起航：毫米波雷达开发手册之硬件配置
眼观四海：自动驾驶&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

参考文献

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