一統江湖:毫米波雷達開發手冊之大話線譜估計

2023-05-28 12:00:08

寫在前面

​ 深知新手在接觸毫米波雷達板硬體時需要花費的沉沒成本,因此在行將告別毫米波雷達之際,總結這兩年以來在毫米波雷達上的一些經驗和教訓。

​ 本檔案用於為實現基於AWR1243BOOST等單板毫米波雷達開發提供參考指南與解決方案,主要包括硬體設定基礎引數訊號模型應用DEMO開發以及可深入研究方向思考等;為更好地匹配後續級聯雷達應用的學習路線,在本手冊中會盡可能同化單板雷達和級聯雷達中的相關表述。

​ 本指南作者資訊:Xuliang,聯絡方式:[email protected]。未經本人允許,請勿用於商業和學術用途。

​ 希望後者在使用本指南時可以考慮參照作者在毫米波雷達旅途中的相關工作,如本文參考文獻[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.