炉火纯青：毫米波雷达开发手册之大话空间谱估计

严肃声明

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

写在前面

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

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

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

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

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

空间谱估计算法

信号模型

空间谱估计是利用空间阵列实现空间信号的参数估计的技术，空间谱估计系统由空间信号入射、空间阵列接收以及参数估计等三部分组成，相应地可分为三个空间即目标空间、观察空间以及估计空间。

目标空间通常是一个由信号源的参数与复杂环境参数张成的空间。观察空间是利用空间按一定方式排列的阵元来接收目标空间的辐射信号，接收数据中往往包含信号特征（方位、距离、极化等）和空间环境特征（噪声、杂波、干扰等），观察空间是一个多维空间。估计空间是利用空间谱估计技术从复杂的观察数据中提取信号的特征参数。

阵列接收信号为$s(t)$，目标空间源信号载波为$\exp(j\omega t)$，信号在空间沿波束向量$k$的方向传播，基准点处信号为$s(t)\exp(j\omega t)$，则距离基准点处的阵元接收信号为$s_{r}(t)=s(t-\frac{1}{c}\boldsymbol{r}^{T}/t)e^{j(\omega t-r^{T}k)}$.

对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(t)e^{j\omega_0\tau_m(\Phi_i)} + n_m(t)$

其中，$s_i(t)$表示入射到阵列的第$i$个源信号，$n_m(t)$表示第个阵元的加性噪声，$\tau_m(\Phi_i)$为来自$\Phi_i$方向的源信号投射至第$m$个阵元时相对选定参考点的时延，$\theta_i$表示俯仰角，$\phi_i$表示方位角。在AWR1243BOOST等毫米波雷达中需要分析方位图、点云图，方位角和俯仰角是实现这些分析的基础。

根据第$m$个阵元的输出，可以推广至空间的表示，记观察空间为$X(t)$,噪声空间为$N(t)$，流形为$A(\Phi)$, $S(t)$为源信号。那么有如下：

$X(t)=[x_1(t),x_2(t),...,x_M(t)]^T$

$N(t)=[n_1(t),n_2(t),...,n_M(t)]^T$

$S(t)=[s_1(t),s_2(t),...,s_M(t)]^T$

$A(\Phi)=[a(\Phi_1),a(\Phi_2),...,a(\Phi_K)]$

观察空间可用噪声空间、目标空间和流形表示，即满足下式：$\boldsymbol{X}(\boldsymbol{t})=\boldsymbol{A}(\boldsymbol{\Phi})S(\boldsymbol{t})+N(\boldsymbol{t})\overline{}$

流形$\boldsymbol{A}(\boldsymbol{\Phi})$与阵列的形状、信号源的来向有关，通常在实际应用中天线阵的形状一旦固定则不会发生改变，因此流形任意一列总是和目标空间源信号的来向有关。阵列形状通常有均匀线阵、均匀圆阵、L型线阵、平面阵列、任意阵列等。

名词解释

波数向量：$|\boldsymbol{k}|=\frac{\omega}{c}=\frac{2\pi}{\lambda}$，空间距离变化1m时相位的变化量

信号相对于基准点的延时时间：$\frac{1}{c}\boldsymbol{r}^T\boldsymbol{a}$

电磁波传播至离基准点$r$处的阵元相对于电磁波传播到基准点的滞后相位：$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$等价的。

无论是相对运动还是静止，固定帧和距离单元后，由目标空间入射的源信号方向可以认为不发生变化或保持平稳随机，统计特性不随时间变化，因此定义阵列的协方差矩阵$R$为：

$\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$

${\boldsymbol m}_{{\boldsymbol x}}({\boldsymbol t})={\boldsymbol E}[{\boldsymbol X}({\boldsymbol t})],{\boldsymbol m}_{{\boldsymbol x}}({\boldsymbol t})={\boldsymbol0}$

$\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(\Phi)S(\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$

阵列信号$X(t)$的协方差可以表示为目标子空间与噪声子空间的加形式，也可以通过特征分解表示为$M$个特征值与特征矢量的积和形式。

$\boldsymbol\Sigma=\mathbf{diag}(\boldsymbol\lambda_1,\boldsymbol\lambda_2,...,\boldsymbol\lambda_M),$通过特征值排序，$\boldsymbol{\lambda}_{1}\geq\boldsymbol{\lambda}_{2}\geq\cdots\geq\boldsymbol{\lambda}_{K}\succ\boldsymbol{\lambda}_{K+1}=\cdots=\boldsymbol{\lambda}_{M}=\boldsymbol{\sigma}^{2},$前$K$个个特征值可以认为是与目标子空间信号相关的，由其对应的特征向量可以表示目标信号子空间$U_s$，后$M-K$个特征值完全取决于噪声，数值等于$\sigma^2$，，由其对应的特征向量构成噪声子空间$U_N$。因此，阵列信号$X(t)$的协方差可以进一步表示为：$\boldsymbol{R}=\boldsymbol{U\Sigma}U^H=\boldsymbol{U_S}\boldsymbol{\Sigma_S}U_S^H+\boldsymbol{U_N}\boldsymbol{\Sigma_N}U_N^H$

在下面，我们将对子空间方法和压缩感知方法进行相对全面的推导、分析和实现，考虑到部分方法的内容较多，笔者采用了贴笔记的方式代替MARKDOWN重新手打，以减少笔者的负担。

子空间方法

MUSIC

这种将阵列信号的协方差矩阵进行特征分解，得到与信号分量对应的信号子空间和信号分量正交的噪声子空间，利用两个子空间的正交性来估计信号参数。根据信号子空间与噪声子空间的正交性，可以得到以下表述：

$\boldsymbol{R}U_N=[U_S,U_N]\boldsymbol{\Sigma}\begin{bmatrix}U_S^H\\U_N^H\end{bmatrix}\boldsymbol{U}_N=[U_S,U_N]\boldsymbol{\Sigma}\begin{bmatrix}\boldsymbol{O}\\I\end{bmatrix}=[U_S,U_N]\begin{bmatrix}\boldsymbol{\sigma}_S^2&\cdot\\.&\boldsymbol{\sigma}_N^2\end{bmatrix}\begin{bmatrix}\boldsymbol{O}\\I\end{bmatrix}=\sigma_N^2U_N$

$A(\Phi)R_{s}A(\Phi)^{H}U_{N}=0$

$U_N^H A(\Phi)R_s A(\Phi)^H U_N=(A(\Phi)^H U_N)^H R_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(t)a_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}$

因信号子空间和噪声子空间相互正交，可知信号子空间的导向矢量$A(\Phi)$与噪声子空间正交，但实际由于干扰$\boldsymbol{a}^{H}(\boldsymbol{\Phi})U_{N}=\boldsymbol{0}$不完全成立，即不完全满足正交性。故采用最小优化搜索（零谱搜索）来估计波达方向：$\Phi_{\text{music}}=\boldsymbol{arg_\Phi min}a^H(\boldsymbol\Phi)U_N U_N^H a(\boldsymbol\Phi)$

根据信号子空间的导向矢量与噪声子空间正交的原理，信号入射方向上会出现极小值，因此空间谱函数可以表示为：$\mathbf{P}_{\text{music}}(\theta)=\frac{1}{\boldsymbol{a}^H(\Phi)U_N U_N^H\boldsymbol{a}(\Phi)}$

那么通过遍历$\theta$，当$\mathbf{P}_{\text{music}}(\theta)$由极大值时，对应角为估计的角度。

下面给出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

区别于MUSIC算法利用自协方差矩阵中噪声子空间和信号子空间正交的特性，ESPRIT利用了自协方差矩阵信号子空间的旋转不变特性。
主要原理：利用ULA中相邻子阵间存在固定间距，固定间距反映了各相邻子阵间的一个固定关系（子阵间的旋转不变性）,下面给出了基于最小二乘思想的谱推导,其中$X_1$和$X_2$为两个子阵，满足$\Phi$的相位旋转特性。

$$\boldsymbol{X}=\left[ \begin{array}{c} \boldsymbol{X}_1\\ \boldsymbol{X}_2\\ \end{array} \right] =\left[ \begin{array}{c} \boldsymbol{A}\\ \boldsymbol{A}\Phi\\ \end{array} \right] \boldsymbol{S}+\boldsymbol{N}=\bar{A}\boldsymbol{S}+\boldsymbol{N} \\ \boldsymbol{R}=\boldsymbol{E}\left[ \boldsymbol{XX}^H \right] =\bar{A}\boldsymbol{S}\bar{A}^H+\boldsymbol{R}_N \\ \hat{R}=\boldsymbol{U}_S\Sigma _S\boldsymbol{U}_{S}^{H}+\boldsymbol{U}_{\boldsymbol{N}}\Sigma _{\boldsymbol{N}}\boldsymbol{U}_{\boldsymbol{N}}^{H} \\ \boldsymbol{U}_S=\left[ \begin{array}{c} \boldsymbol{U}_{S1}\\ \boldsymbol{U}_{S2}\\ \end{array} \right] =\left[ \begin{array}{c} \boldsymbol{AT}\\ \boldsymbol{A}\Phi \boldsymbol{T}\\ \end{array} \right] \,\, \boldsymbol{U}_{\boldsymbol{S}2}=\boldsymbol{U}_{\boldsymbol{S}1}\boldsymbol{T}^{-1}\Phi \boldsymbol{T}=\boldsymbol{U}_{\boldsymbol{S}1}\Psi \\ $$

为更好地理解ESPRIT，我们给出了演示代码，并和总体最小二乘约束下的TLS-ESPRIT进行了对比。

%% 本程序对比了ESPRIT算法的不同形式
%% By Xuliang
clc;clear;close all;
% 参数初始化
M = 12; % 阵元数目
a = [0:M-1]'; % 导向序号
 
f0 = 77e6; % 频率
c = 3e8; % 光速
lambda = c / f0; % 波长
d = lambda / 2; % 阵元间距
 
snap = 128; % 快拍数
fs = 1000; % 采样频率
t = 1 / fs * (0:snap-1);  % 时间
 
% 线性调频信号生成
P = 5; % 信号源数目.
thetas = [-30 5 10 30 50 40 60 70 20]; % 信号源方向
s = (randn(P, snap) + 1j * randn(P, snap)) / sqrt(2);  % 可以通过随机数来控制相位
% s = exp(1j * randn(P, snap) * 2 * pi);
a_theta = exp(1j * 2 * pi * d / lambda * a * sind(thetas)); % 信号导向矢量
 
X0 = a_theta(:,1:P) * s(1:P,:); 
 
theta_step = 0.1; % 遍历网格步长
theta_grids = -90 : theta_step : 90; % 遍历网格区间
snr = -10 : 1 : 30; % 信噪比
 
 
% ESPRIT参数设计
ESP_NUM = 4; % 每个子阵选取阵元数目为M-ESP_NUM
 
esprit_angle1 = zeros(length(snr), length(theta_grids));
esprit_angle2 = zeros(length(snr), length(theta_grids));
monte_num = 200;
Pd_esp1 = zeros(1, length(snr));
Pd_esp2 = zeros(1, length(snr));
Pe_esp1 = zeros(1, length(snr));
Pe_esp2 = zeros(1, length(snr));
Pv_esp1 = zeros(1, length(snr));
Pv_esp2 = zeros(1, length(snr));
 
tic
for i = 1 : length(snr)
   esp1_count_num = zeros(1, monte_num);
   esp1_estimated_bias = zeros(1, monte_num);
   esp2_count_num = zeros(1, monte_num);
   esp2_estimated_bias = zeros(1, monte_num);
   
   for monte_idx = 1 : monte_num
       X = awgn(X0, snr(i), 'measured'); % 完整基带信号
   
       % ESPRIT方法
       % LS-ESPRIT
       ESP_X1 = X(1:M-ESP_NUM, :);
       ESP_X2 = X(2:M-ESP_NUM+1, :);
       ESP_XX = [ESP_X1; ESP_X2]; % 维度变成 2 * (M-ESP_NUM)
       ESP_R1 = ESP_XX * ESP_XX' / snap; % 两个子阵的自协方差矩阵
       [ESP_V1, ~] = eig(ESP_R1);
       ESP_US = ESP_V1(:, 2*(M-ESP_NUM)-P+1 : 2*(M-ESP_NUM));
       ESP_US1 = ESP_US(1:(M-ESP_NUM),:);
       ESP_US2 = ESP_US(1+(M-ESP_NUM):2*(M-ESP_NUM),:);
       ESP_ANS1 = inv(ESP_US1' * ESP_US1) * ESP_US1' * ESP_US2; % LS-ESPRIT
       [ESP_V1, ESP_D1] = eig(ESP_ANS1);
       ESP_D1 = diag(ESP_D1);
       ESP_A1 = angle(ESP_D1);
       ESP_A1 = sort(asin(ESP_A1/pi)/pi*180);
       
       % TLS-ESPRIT
       ESP_US12 = [ESP_US1 ESP_US2];
       [ESP_V2, ~] = eig(ESP_US12' * ESP_US12);
       ESP_EN2 = ESP_V2(:,1:P); % 获取噪声子空间
       ESP_ANS2 = -ESP_EN2(1:P,:) * inv(ESP_EN2(P+1:2*P,:)); % TLS-ESPRIT
       [ESP_V2, ESP_D2] = eig(ESP_ANS2);
       ESP_D2 = diag(ESP_D2);
       ESP_A2= angle(ESP_D2);
       ESP_A2 = sort(asin(ESP_A2/pi)/pi*180);
 
       for target_idx = 1 : P
           for detected_idx = 1 : length(ESP_A1)
               if abs(ESP_A1(detected_idx) - thetas(target_idx)) <= 2.0
                   esp1_count_num(1, monte_idx) = esp1_count_num(1, monte_idx) + 1;
                   esp1_estimated_bias(1, monte_idx) = abs(ESP_A1(detected_idx) - thetas(target_idx));
               end
           end
           
           for detected_idx = 1 : length(ESP_A2)
               if abs(ESP_A2(detected_idx) - thetas(target_idx)) <= 2.0
                   esp2_count_num(1, monte_idx) = esp2_count_num(1, monte_idx) + 1;
                   esp2_estimated_bias(1, monte_idx) = abs(ESP_A2(detected_idx) - thetas(target_idx));
               end
           end
       end
       
   end
   Pd_esp1(1,i) = sum(esp1_count_num) / monte_num / P; % 计算准确率
   Pe_esp1(1,i) = sqrt(mean(esp1_estimated_bias.^2)); % 计算均方根误差
   Pv_esp1(1,i) = var(esp1_estimated_bias); % 计算估计方差
   Pd_esp2(1,i) = sum(esp2_count_num) / monte_num / P; % 计算准确率
   Pe_esp2(1,i) = sqrt(mean(esp2_estimated_bias.^2)); % 计算均方根误差
   Pv_esp2(1,i) = var(esp2_estimated_bias); % 计算估计方差
end
toc
 
figure(1);
axes('Fontsize', 14);
plot(snr, Pd_esp1, 'Color','#006400','LineWidth',1.5,'Marker','+');hold on;
plot(snr, Pd_esp2, 'Color','#8B0000','LineWidth',1.5,'Marker','d');hold on;
xlabel('\fontname{Times New Roman}SNR (dB)');ylabel('\fontname{Times New Roman}Detected Probability');
legend('LS-ESPRIT','TLS-ESPRIT');
 
figure(2);
axes('Fontsize', 14);
plot(snr, Pe_esp1, 'Color','#006400','LineWidth',1.5,'Marker','+');hold on;
plot(snr, Pe_esp2, 'Color','#8B0000','LineWidth',1.5,'Marker','d');hold on;
xlabel('\fontname{Times New Roman}SNR(dB)');ylabel('\fontname{Times New Roman}Estimated Error');
legend('LS-ESPRIT','TLS-ESPRIT');
 
figure(3);
axes('Fontsize', 14);
plot(snr, Pv_esp1, 'Color','#006400','LineWidth',1.5,'Marker','+');hold on;
plot(snr, Pv_esp2, 'Color','#8B0000','LineWidth',1.5,'Marker','d');hold on;
xlabel('\fontname{Times New Roman}SNR(dB)');ylabel('\fontname{Times New Roman}Estimated Variance');
legend('LS-ESPRIT','TLS-ESPRIT');
 

DML

确定性最大似然估计方法DML和随机最大似然估计方法SML均为子空间拟合方法的典型代表。主要思想是将阵列流型矩阵与接收数据的子空间进行拟合，通过将观测信号的似然函数定义为含有未知参数的条件概率密度函数，使选定未知的参数使似然函数尽可能大。
本处，主要介绍DML方法，其将噪声过程建模为平稳高斯随机白噪声过程，并假设信号波形是确定性信号（但输入波形是需估计的参数），下面给出具体的推导。

function PML = DOA_ML(M, P, X)
    % DML/SML的一维扫描空间谱实现[希望有大佬能教一下K维度的]
    % M：阵列阵元数目
    % P：信源数目
    % X：基带信号
    % By Xuliang
    
    f0 = 77e9; % 本振频率
    a = [0:M-1]'; % 阵列序号
    c = 3e8; % 光速
    lambda = c / f0; % 波长
    d = lambda / 2; % 阵元间距
 
    % 特征值分解
    RX = X * X' / size(X,2);
    [VV, DD] = eig(RX);
    DD = diag(DD);
    US = VV(:, M-P+1:M);
    UN = VV(:, 1:M-P);theta_step = 1; % 遍历网格步长
    DS = US(M-P+1:M, :);
    sigma2 = mean(DD(1:M-P)); % 通常采用噪声特征值的平均值表示估计
    WSOPT = (diag(DD(M-P+1:M)) - sigma2*eye(P)) * pinv(diag(DD(M-P+1:M))); % 最优权重
 
    theta_grids = -90 : theta_step : 90; % 遍历网格区间
 
    for j = 1 :length(theta_grids)
        theta_vecs = exp(-1j * 2 * pi * d / lambda * a * sind(theta_grids(j)));
         Pa = theta_vecs * inv(theta_vecs' * theta_vecs) * theta_vecs'; % 投影矩阵
         A_pos = inv(theta_vecs' * theta_vecs) * theta_vecs';
         WNOPT = A_pos * US * WSOPT * US' * A_pos';
         PML(j) = db(abs(trace(Pa * RX))); % DML和SML的最大化功率谱形式
    end
 
end
 

压缩感知

网格模型

OMP

function PoutOMP = DOA_OMP(y, A, P)
    % y ：测量信号
    % A ：测量矩阵
    % P ：信号数目
    % By Xuliang
    
    k = size(A, 2); % 原子数目
    n = size(A, 1); % 阵元数目
    X = zeros(1, k); % 待估计信号
    rn = y; % 残差初始化为y,n*1
    
    Sup_set = []; % 支持集
    Matched_S = []; % 匹配信号
    while true
        inner_product = A' * rn; % 求解内积 k * 1
        [~, pos] = max(inner_product); % 求解最大内积值
        Sup_set = union(Sup_set, pos); % 更新支持集索引
        Atom_Matrix = A(:, Sup_set); % 存放原子矩阵 n * k
%         A(:, pos) = zeros(n, 1); % 将测量矩阵pos列置0 防止再用
        Matched_S = pinv(Atom_Matrix' * Atom_Matrix) * Atom_Matrix' * y; % 基于LS准则更新
        rn = y - Atom_Matrix * Matched_S; % 更新残差项
        if length(Matched_S) == P
           break; 
        end
    end
    theta_all = 1 : k;
    PoutOMP = ones(1, k); % 将扫描角度均置1
    PoutOMP(setdiff(theta_all, Sup_set)) = 0; % 将A中存在B中不存在的数据置为0
    
end

L1-Norm

这个类似LASSO问题，同样可以参考先前写的范数问题，不过多赘述。

function PoutL1Norm = DOA_L1Norm(X, A)
    % 本程序为L1-Norm的函数实现文件
    % X ： 基带信号
    % A ：过完备基
    % By Xuliang
    
    snap = size(X, 2); % 快拍数目
    M = size(X, 1); % 阵元数目
    thetaNum = size(A, 2); % 原子数目
    cvx_begin quiet
        variables p q
        variable SSV1(thetaNum, snap) complex
        variables r(thetaNum)
        expressions Rn(M, snap) complex
 
        minimize(p + 2 * q); 
        subject to
            Rn = X - A * SSV1; % 求残差
            Rvec = vec(Rn); % 矩阵转向量
            norm(Rvec) <= p; % 第一个不等式约束
            sum(r) <= q; % 第二个不等式约束
            for i = 1 : thetaNum
                norm(SSV1(i, :)) <= r(i);
            end
    cvx_end
    PoutL1Norm = abs(SSV1(:, :)  / max(SSV1(:, :))); % 输出功率谱
end
 

L1-SVD

function PoutSVD = DOA_L1SVD(X, A, P)
    % 本程序为L1-SVD的函数实现文件
    % X ：基带信号
    % A ：过完备基
    % P ： 信源数目
    % By Xuliang
    
    [M, snap] = size(X); % M 阵元 snap 快拍
    thetaNum = size(A, 2); % 扫描网格点数
    DK1 = eye(P);
    DK2 = zeros(P, snap-P);
    DK = [DK1, DK2].'; % SNAP * P的选择矩阵
    [U, Sigm, V] = svd(X); % 奇异值分解
    Xsv = X * V * DK; % 获取新的接收矩阵
    
    % 低信噪比测试方案
    cvx_begin quiet
        variables p q
        variables r(thetaNum)
        variable SSV1(thetaNum, P) complex
        expressions Rn(M, P) complex
 
        minimize(p + 2.7 * q); % 优化目标
        subject to
            Rn = Xsv - A * SSV1; % 求残差
            Rvec = vec(Rn); % 矩阵转换为向量
            norm(Rvec) <= p; % 第一个不等式约束
            sum(r) <= q; % 第二个不等式约束
            for i = 1 : thetaNum % 第三个不等式约束
                norm(SSV1(i, :)) <= r(i);
            end
    cvx_end
    
    % 高信噪比测试方案
    % confidence_interval = 0.9; % 置信值
    % noise = X - X0; % 噪声估计
    % noise_var = var(noise(:)); % 噪声方差估计
    % regulari_param = compute_regulariParam(confidence_interval, noise_var, M, P); % 根据卡方分布反演门限值
    % cvx_begin quiet
    %     variable SSV1(length(theta_grids), P) complex
    %     minimize(sum(norms(SSV1, 2, 2)))
    %     subject to
    %         norm(Xsv - A * SSV1, 'fro') <= regulari_param 
    % cvx_end
    PoutSVD = abs(SSV1(:, :)  / max(SSV1(:, :))); % 求解功率谱
    
end

function PoutSVD = DOA_L1SVD(X, A, P)
    % 本程序为L1-SVD的函数实现文件
    % X ：基带信号
    % A ：过完备基
    % P ： 信源数目
    
    [M, snap] = size(X); % M 阵元 snap 快拍
    thetaNum = size(A, 2); % 扫描网格点数
    DK1 = eye(P);
    DK2 = zeros(P, snap-P);
    DK = [DK1, DK2].'; % SNAP * P的选择矩阵
    [U, Sigm, V] = svd(X); % 奇异值分解
    Xsv = X * V * DK; % 获取新的接收矩阵
    
    % 低信噪比测试方案
    cvx_begin quiet
        variables p q
        variables r(thetaNum)
        variable SSV1(thetaNum, P) complex
        expressions Rn(M, P) complex
 
        minimize(p + 2.7 * q); % 优化目标
        subject to
            Rn = Xsv - A * SSV1; % 求残差
            Rvec = vec(Rn); % 矩阵转换为向量
            norm(Rvec) <= p; % 第一个不等式约束
            sum(r) <= q; % 第二个不等式约束
            for i = 1 : thetaNum % 第三个不等式约束
                norm(SSV1(i, :)) <= r(i);
            end
    cvx_end
    
    % 高信噪比测试方案
    % confidence_interval = 0.9; % 置信值
    % noise = X - X0; % 噪声估计
    % noise_var = var(noise(:)); % 噪声方差估计
    % regulari_param = compute_regulariParam(confidence_interval, noise_var, M, P); % 根据卡方分布反演门限值
    % cvx_begin quiet
    %     variable SSV1(length(theta_grids), P) complex
    %     minimize(sum(norms(SSV1, 2, 2)))
    %     subject to
    %         norm(Xsv - A * SSV1, 'fro') <= regulari_param 
    % cvx_end
    PoutSVD = abs(SSV1(:, :)  / max(SSV1(:, :))); % 求解功率谱
    
end

OGSBI

function [PoutOGSBI,theta_rec] = DOA_OGSBI(X, theta_grids, P, svd_flag)
    % 本程序为OGSBI的函数实现文件
    % X ： 基带信号
    % A ：过完备基
    % P ： 信源数目
    % svd_flag ： 是否进行SVD分解
    % By Xuliang
    
    f0 = 77e9; % 频率
    cspeed = 3e8; % 光速
    lambda = cspeed / f0; % 波长
    d = lambda / 2; % 阵元间距
    [M, snap] = size(X);
    thetaNum = length(theta_grids);
    theta_step = theta_grids(2) - theta_grids(1);
    for m_id = 0 : M - 1
        for n_id = 1 : thetaNum
            temp = exp(-1j * 2 * pi * d / lambda * m_id * sind(theta_grids(n_id)));  % 构造扫描矩阵 M * thetaNum
            A(m_id + 1, n_id) = temp;
            B(m_id + 1, n_id) = (-1j * 2 * pi * d / lambda * m_id * cosd(theta_grids(n_id))) * temp; % 扫描矩阵的导数 逼近误差
        end
    end
    BHB = B' * B;
    
    if svd_flag % 对多快拍MVM信号进行SVD降维
        [U, Sigm, V] = svd(X, 'econ');
        X = X * V(:, 1:P);
    end
    snap = size(X,2); % 更新SVD后的快拍
    
    tolerance = 1e-3; % 收敛阈值
    iter_num = 1e3; % 迭代次数
    sigma2 = mean(var(X)) / 100; % 噪声功率初始化
    alpha0 = 1 / sigma2; % 噪声功率的导数
    alpha0_vec = zeros(iter_num, 1); % 存放每一次迭代的噪声功率参数
    alpha = mean(abs(A' * X), 2); % 初始化信号功率 thetaNum * P
    beta = zeros(thetaNum, 1); % 网格误差
    converged = false; % 收敛判断标记
    iter_beta = 1; % beta更新次数
    iter_id = 0; % 收敛次数
    idx = []; % 用于存放目标谱峰的索引
    % 信号功率参数的伽马分布先验
    rho = 1e-2 ;
    % 噪声功率参数的伽马分布先验
    c = 1e-4; 
    d = 1e-4;
    while ~converged
       iter_id = iter_id + 1; % 迭代序号
 
       Phi = A; % 修正后的扫描矩阵
       Phi(:, idx) = A(:, idx) + B(:, idx) * diag(beta(idx));
       alpha_old = alpha; % 保留上一次的信号功率向量
       C = 1 / alpha0 * eye(M) + Phi * diag(alpha) * Phi'; % 更新证据分布的协方差矩阵
       Cinv = inv(C); % 求逆
       Sigma = diag(alpha) - diag(alpha) * Phi' * Cinv * Phi * diag(alpha); % 获取后验分布的协方差矩阵
       mu = alpha0 * Sigma * Phi' * X; % 后验分布的均值 thetaNum * P
       
       % 更新信号功率参量alpha
       % rho足够小,可用泰勒展开逼近 alpha = E{X^2} = [E(X)]^2+Var(X)
       alpha = mean(abs(mu).^2,2) + real(diag(Sigma)); % 形状和尺度参数通常为实值
       if rho ~= 0
           alpha = (-snap + sqrt(snap^2 + 4 * snap * rho * alpha)) / (2 * rho);
       end
 
       % 更新噪声参量alpha0
       gammas = 1 - real(diag(Sigma)) ./ (alpha + eps); % E{|Y-Phi*X|}的方差项
       alpha0 = (snap * M + c - 1) / (norm(X - Phi * mu, 'fro')^2 + snap / alpha0 * sum(gammas) + d); 
       alpha0_vec(iter_id) = alpha0; % 存入alpha0_vec
 
       % 收敛条件判断
       if norm(alpha - alpha_old) / norm(alpha_old) < tolerance || iter_id >= iter_num
           converged = true;
           iter_beta = 5;
       end
 
       [temp, idx] = sort(alpha, 'descend'); % 降序排序
       idx = idx(1 : P); % 找到前P个最大峰值对应的entry
       beta = zeros(thetaNum, 1); % 谱存在稀疏性,因此将其他位置置为0,
       beta(idx) = temp(idx);
 
       % 更新beta因子
       Ps = real(conj(BHB(idx, idx)) .* (mu(idx, :) * mu(idx, :)' + snap * Sigma(idx, idx)));
       v = zeros(length(idx), 1); 
       for t = 1 : snap
           v = v + real(conj(mu(idx, t)) .* B(:, idx)' * (X(:, t) - A * mu(:, t)));
       end
       v = v - snap * real(diag(B(:, idx)' * A * Sigma(:, idx)));
       optim_beta = P \ v; % 基于LS准则求解满足最优条件的beta
       if any(abs(optim_beta) > theta_step / 2) || any(diag(Ps) == 0)
           for iter_i = 1 : iter_beta
              for iter_n = 1 : P
                  temp_beta = beta(idx); 
                  temp_beta(iter_n) = 0; % 去除了第n个条目的向量
                  beta(idx(iter_n)) = (v(iter_n) - Ps(iter_n, :) * temp_beta) / Ps(iter_n, iter_n);
 
                  if beta(idx(iter_n)) > theta_step / 2 % 约束角度上限
                      beta(idx(iter_n)) = theta_step / 2;
                  end
 
                  if beta(idx(iter_n)) < - theta_step / 2 % 约束角度下限
                      beta(idx(iter_n)) = - theta_step / 2;
                  end
 
                  if Ps(iter_n, iter_n) == 0
                      beta(idx(iter_n)) = 0;
                  end
 
              end
           end
       else
           beta = zeros(thetaNum, 1);
           beta(idx) = optim_beta; % 如果optim_beta介于-theta_step/2和theta_step/2之间
       end 
    end
    
    theta_rec = theta_grids' + beta;
    if svd_flag % 使用svd则需要复原
        x_rec = mu * V(:,1:size(mu,2))';
        xpower_rec = mean(abs(x_rec).^2,2) + real(diag(Sigma)) * P / snap;
    else
        xpower_rec = mean(abs(mu).^2,2) + real(diag(Sigma));
    end
    PoutOGSBI = abs(xpower_rec / max(xpower_rec));
end

function PoutRootSBI = DOA_RootSBI(X, SearchArea, etc)
    % 本程序为ROOTSI的函数实现文件
    % X ： 基带信号
    % SearchArea ：搜索网格
    % etc ： 控制稀疏向量的超参
    % By Xuliang,20230226
    
    [M, snap] = size(X); % 基带信号维度 阵元数 * 快拍数
    K_hat = length(SearchArea); % 搜索区域的长度/过完备基长度
    resolution = SearchArea(2) - SearchArea(1); % 网格分辨率
    SearchMidleft = SearchArea - resolution / 2; % 左边界
    SearchMidright = SearchArea + resolution / 2; % 右边界
    f0 = 77e9; % 频率
    c = 3e8; % 光速
    lambda = c / f0; % 波长
    d = lambda / 2; % 阵元间距
 
    A_theta = exp(-1j * 2 * pi * (0:M-1)' * d / lambda * sind(SearchArea)); % 过完备基
    a = 1e-4; b = 1e-4; % a b为控制功率参量伽马分布的双因子
    rho = 1e-2; % 控制稀疏向量\delta满足伽马分布的因子 Delta=diag(delta)
    max_iter = 500; % 最大迭代次数
    tolerance = 1e-4; % 收敛阈值
    
    % 参数初始化
    sigma2 = mean(var(X)) / 100; % 噪声功率初始化
    beta = 1 / sigma2; % 噪声功率的倒数，置信度
    delta = mean(abs(A_theta' * X), 2); % 初始化信号功率 thetaNum * P
 
    converged = false; % 收敛判断条件
    iter = 0; % 迭代次数
    
    while ~converged
       iter = iter + 1; 
       delta_last = delta;
       
       % 计算最大后验估计得到的信号S的均值和协方差
       Phi = A_theta;
       Sigmat = 1 / beta * eye(M) + Phi * diag(delta) * Phi'; % 信号X的协方差
       inv_Sigmat = inv(Sigmat);
       Sigma = diag(delta) - diag(delta) * Phi' * inv_Sigmat * Phi * diag(delta); % 协方差的维度为网格数*网格数
       mu = beta * Sigma * Phi' * X; % 均值的维度为网格数*快拍数
       diagSigma = 1 - real(diag(Sigma)) ./ delta; % X-Phi*S的协方差展开项
       
       % 更新控制稀疏向量功率参数的delta
       tempT = sum(mu .* conj(mu), 2) + snap * real(diag(Sigma));
       delta = (-snap + sqrt(snap^2 + 4 * rho * real(tempT))) / (2 * rho);
       
       % 更新功率置信度beta
       resiE = X - Phi * mu; % M * T维度值
       beta = (snap * M + a - 1) / (b + norm(resiE, 'fro')^2 + snap / beta * sum(diagSigma));
       
       % 判断是否收敛
       errors = norm(delta - delta_last) / norm(delta_last);
       if errors < tolerance || iter >= max_iter
           converged = true;
       end
       
       % 根细化
       f = sqrt(sum(mu .* conj(mu),2));
       [~, sort_ind] = sort(f); % 升序排序
       index_amp = sort_ind(end:-1:end-etc+1); % 选取排序后的顺序 最大值-->次大值
       
       for j = 1 : length(index_amp)
          idx = index_amp(j);
          mut = mu(idx, :); % 获取某个网格值对应的均值-快拍矢量
          gammat = Sigma(:, idx); % 获取某个网格值对应的协方差矢量
          phik = mut * mut' + snap * gammat(idx); % 第k个网格对应的phi值
          
          tempind = [1 : K_hat];
          tempind(idx) = []; % 满足j不等于i的条件 去除与i相等项目
          Xti = X - Phi(:, tempind) * mu(tempind, :);
          psik = snap * Phi(:, tempind) * gammat(tempind) - Xti * mut'; % 维度为M*1
      
          z1 = [1:M-1]'; % 根系数
          cvec = zeros(M, 1); % 初始化根矩阵
          cvec(1) = M * (M - 1) / 2 * phik;
          cvec(2:end) = z1 .* psik(2:end);
          
          ro = roots(cvec); % 求根
          abs_root = abs(ro);
          [~, indmin] = min((abs(abs_root-1))); % 最靠近单位圆上的为根
          angle_vec = asind(-angle(ro(indmin)) * lambda / (2 * pi * d));
          
          if angle_vec <= SearchMidright(idx) && angle_vec >= SearchMidleft(idx)
              SearchArea(idx) = angle_vec;
              A_theta(:, idx) = exp(-1j * 2 * pi * (0:M-1)' * d / lambda * sind(angle_vec));
          end
       end   
    end
    PoutRootSBI = mean(mu .* conj(mu), 2);
    PoutRootSBI = abs(PoutRootSBI / max(PoutRootSBI));
end
 

L1-SARCV

根据线性代数理论，阵列协方差矩阵$\boldsymbol{R}\in\mathbb{R}^{M\times M}$的每一列可以用$M$维复向量空间中的完备基任意表示。

过完备基$\{a\big(\theta_q\big)\}_{q=1}^Q(Q\gg M),\theta_q$为空间域采样的指定方向(以$1°$间隔在$-90°$至$90°$采样）。基于此，重新表示协方差矩阵$R$的第$m$列：

$$\tau_m=E\Big(x\big(t)x_m^*(t)\Big)=A\big(\theta\big)b_m+\sigma^2e_m,m=1\text{,}2,...,M$$

$A(\theta) = [a(\theta_1),a(\theta_2),...,a(\theta_Q)]$是$M\times Q$的阵列流形矩阵，$b_m$是过完备基对应的稀疏表示系数矢量。如果$\left\{\theta_{q}\right\}_{q=1}^{Q}$足够密集，$\left\{a{\left(\theta_{q}\right)}\right\}_{q=1}^{Q}$可以被认为和$\{a(\theta_k)\}_{k=1}^K$非常接近，理想的$b_m$应该是除了与基向量相关的$K$个元素外其他元素均为0，这意味$b_m$具有和信号DOA相关的稀疏结构（非零模式）。

将$r_m$对应的模型扩展到矩阵形式(SRACV模型)：

$$\boldsymbol{R}=A(\boldsymbol{\theta})B+\sigma^2I_M$$

$\boldsymbol{B}=[\boldsymbol{b}_{1},\boldsymbol{b}_{2},...,\boldsymbol{b}_{M}],\{b_{m}\}_{m=1}^M]$共享相同的稀疏结构(每个$b_m$的非零元素出现在$B$的相同行）；然后，通过引入$\boldsymbol{b}^{\circ}=\Bigl[\boldsymbol{b}_{1}^{\circ},\boldsymbol{b}_{2}^{\circ},...,\boldsymbol{b}_{Q}\Bigr]^{T},$其中$\left.b_{q}^{\circ}=\left\|(B)_{q^{c}}\right.\right\|_{2}$（第$q$个稀疏稀疏满足$B$的第$q$行矢量的$L_2$范数），这意味着$\{b_m\}_{m=1}^M$的稀疏结构可以用相同稀疏结构的$\boldsymbol{b}^{\circ}$来描述，等价于$\boldsymbol{b}^{\circ}$能够表征过完备基里$\{r_{m}\}_{m=1}^{M}$的系数特性。

在SRACV算法中， $$\boldsymbol{b}^{\circ}$$ 定义为阵列协方差矢量的系数表示。在忽略误差项$\sigma^2 I_M$的前提下，DOA估计问题等价于找到充分稀疏的$\boldsymbol{b}^{\circ}$就可以使$\{b_m\}_{m=1}^M$尽可能稀疏且拟合${r_m}_{m=1}^M$。

评估$\boldsymbol{b}^{\circ}$稀疏性的标准是引入L1范数解决SRACV问题。

Remark：SRACV模型不涉及信号协方差矩阵$R_s$，的秩先验，因此，无需去相关就可用于任意阵列下不相关和相关信号源的估计；估计精度受限于网格的分辨率($\theta$的步长）。

L1-SRACV问题求解可以表征为：

$$min_{B}\left\|b^*\right\|_1,\text{subject to}~R=A(\boldsymbol\theta)B+\sigma^2I_M$$

实际中，$R$未知，由于阵列接收数据长度有限，可以用数据协方差矩阵的最大似然估计来确定$\hat{R}=\frac{\sum_{t=1}^N[x(t)x^H(t)]}{N}=R+\Delta R$来估计，$\Delta R=\widehat{\boldsymbol{R}}-R$表示估计误差，通过沿列方向堆栈获得栈向量，且符合渐近高斯分布。

$$\text{rec}(\Delta R)\sim AsN(0,\frac{1}{N}R^T\otimes R)$$

据上，DOA估计问题可以描述为在拟合数据$vec(\Delta R)$至其数据模型时的最稀疏向量$\boldsymbol{b}^{\circ}$，问题可以进步表示为：

$$\min\hat{re}\left[\hat{R}-A(\theta)B-\sigma^2I_M\right]^HW^{-1}vec\left[\hat{R}-A(\theta)B-\sigma^2I_M\right]+\lambda\left\|\hat{b}^*\right\|_1$$

为了将$vec(\Delta R)$变成标准正态分布，$W^{-\frac12}=\sqrt{N}R^{-\frac12}\otimes R^{-\frac12}$.

那么则有$\|Wvec{R}-A(\theta)B-\sigma^2I_M\bigg]\|_2^2$服从卡方分布，引入参数$\beta$使得下式以很高的概率P成立，P是接近1的数(一般取0.999)

$$\left\|Wvec\left[\hat R-A(\theta)B-\sigma^2I_M\right]\right\|^2_2\le\beta$$

通过上面的推导分析，问题可以进一步优化如下：

$$min_B.\left\|b^\circ\right\|_1,\text{s.t.}\left\|z-\Psi\text{vec(B)}\right\|_2^2\le\beta\\ \text{z}=W^{-\frac{1}{2}}vec[R-\sigma^2I_M],\Psi=W^{-\frac{1}{2}}(I_M\otimes A(\theta))$$

在上面，$\sigma^2$可通过最大似然估计得到（$\hat{R}$中M-K个最小特征值的平均，这意味着需要信号源数目先验；如果不知道先验，则采用最小的特征值折衷)

上述问题，往往需要二阶锥优化来求解。

约束$\left\|b^{\circ}\right\|_{1}\leq g,g$为辅助变量；等价为$\operatorname*{min}_{B,Y,g}g$

根据$b^{\circ}$的每个元素都是$B$的行的L2范数，那么有：

$$\mathbf{1}^T\gamma\leq g,\|(B)_{q}.\|_2\leq\gamma_{q},q=1,2,...Q~and~\|\text{z}-\Psi\text{vec}(\text{B})\|_2^2\leq\beta$$

其中,$\mathbf{1}^T$为$Q\times 1$的单位向量，$\gamma$为第$q$个元素为$\gamma_q$的矢量。

代码复现如下：

function PoutSRACV = DOA_SRACV(X, A, P)
    % 本程序为L1-SRACV的函数实现文件
    % X ：基带信号
    % A ：过完备基
    % P ： 信源数目
    % By Xuliang
    
    [M, snap] = size(X); % 阵元 快拍
    thetaNum = size(A, 2); % 网格数目
    
    R = X * X' / snap; % 信号协方差矩阵
    [V, D] = eig(R); % 特征值分解
    eig_value = real(diag(D)); % 提取特征值
    [B, I] = sort(eig_value,'descend'); % 排序特征值
    EN = V(:, I(P+1:end)); % 提取噪声子空间向量
    diagD = diag(D);
    Rn = EN * diag(diagD(I(P+1:end))) * EN'; % 噪声协方差矩阵
    
    identity_vec = ones(thetaNum, 1); % 创建单位矢量 控制稀疏度
    p = 0.001; % 正则参数
    mu = chi2inv(1-p, M^2); % 根据卡方分布求出mu值
    z = vec(sqrt(snap) * R^(-1/2) * (R - Rn) * R^(-1/2));
    Phi = sqrt(snap) * kron((R^(-1/2)).', (R^(-1/2)*A));
    
    cvx_begin quiet
        variables gammas(thetaNum);
        variables g; % 辅助变量
        variable BB(thetaNum, M) complex; % 稀疏矩阵
        expression bb(thetaNum, 1); % 稀疏向量
 
        minimize(g); % 优化目标
        subject to
        identity_vec' * gammas <= g; % 第一个约束
        bb = cvx(zeros(thetaNum, 1));
        for qid = 1 : thetaNum % 第二个约束
            bb(qid, :) = norm(BB(qid, :));
            bb(qid) <= gammas(qid);
        end
        norm(z - Phi * vec(BB)) <= sqrt(mu); % 第三个约束
    cvx_end
    PoutSRACV = abs(bb(:, 1) / max(bb(:, 1))); % 求解功率谱
    
end
 

IAA

在这里，笔者实现了IAA算法的3种变体形式，但不再对具体的理论部分进行赘述，读者可以通过代码中所给出的参考文献进行验证，如有错误，欢迎在评论区留言指正。

function PoutIAA = DOA_IAA(X, A, params)
    % 本程序用于实现IAA算法的3种变体形式
    % 参考文献:AReduced ComplexityApproach to IAA Beamforming for
    % Efficient DOA Estimation of Coherent Sources
    % By Xuliang
    % X : 输入基带信号
    % A : 字典矩阵
    % params : (mode, iter_num1, iter_num2, beta_thres)
    % mode: 选择算法运行的模式 IAA-APES/IAA-ML/IAA-RC
    % iter_num1 : 第一轮迭代次数  iter_num2 第二轮迭代次数
    % beta_res : 门限系数【针对RC/RCML】
    
    [M, snap] = size(X); % 阵列数目 * 快拍
    thetaNum = size(A, 2); % 原子数目
    
    mode = params.mode;
    threshold = params.threshold;
    if strcmp(mode, "APES")
        P_old = (sum(abs(A' * X / M), 2) / snap).^2;
        iter_num = params.iter_num;
        for iter_id = 1 : iter_num
            R = A * spdiags(P_old, 0, thetaNum, thetaNum) * A'; 
            invR = inv(R); 
 
            P = zeros(thetaNum, 1); 
            for snapIdx = 1 :snap
                x = X(:, snapIdx); 
                invRx = invR * x; % M * 1
                invRa = invR * A; % M * thetaNum
                ainvRa = sum(conj(A) .* invRa, 1).'; 
                coeff = A' * invRx ./ real(ainvRa); 
                P = P + abs(coeff).^2;
            end
            P = P / snap;
 
            if norm(P_old - P) < threshold 
                break;
            end
            P_old = P;
        end
 
    elseif strcmp(mode, "ML")
        Ru = X * X' / snap;
        for k = 1 : thetaNum
            P(k) = pinv(A(:, k)' * pinv(Ru) * A(:, k));
        end
        
        Pold = diag(P);
        Rold = pinv(A * Pold * A');
        for iter_id = 1 : iter_num
            [~, idx] = sort(P);
            for k = 1 : thetaNum
               P_prev(idx(k)) = P(idx(k));
               P(idx(k)) = max(0, P(idx(k)) + (A(:, idx(k))' * Rold * (Ru - pinv(Rold)) * Rold * A(:, idx(k))) / (A(:, idx(k))' * Rold * A(:, idx(k)))^2);
               Rold = Rold - ((P(idx(k)) - P_prev(idx(k))) * Rold * A(:, idx(k)) * A(:, idx(k))' * Rold) * pinv(1 + (P(idx(k)) - P_prev(idx(k))) * A(:, idx(k))' * Rold * A(:, idx(k)));
            end
            
            if norm(P - P_prev) < threshold 
                break;
            end
        end
   
    elseif strcmp(mode, "RC")
        iter_num1 = params.iter_num1;
        iter_num2 = params.iter_num2;
        beta_thres = params.beta_thres;
 
        P = (sum(abs(A' * X / M), 2) / snap).^2;
        Pold = diag(P);
        Ru = X * X' / snap;
        for iter_id = 1 : iter_num1
            Rold = (A * Pold * A'); 
            for k = 1 : thetaNum
               w(k, :) = A(:, k)' * pinv(Rold) / (A(:, k)' * pinv(Rold) * A(:, k));
               P(k) = real(w(k, :) * Ru * w(k, :)');
            end
            Pold = diag(P);
        end
 
        [Pval, idx] = sort(P); % 升序
        [Rrow] = find(Pval > beta_thres * P(idx(1))); % 易误点: 不等式右边为排序后最小元素值乘以门限 左边为排序元素值 目的是找到第一个满足阈值的元素列标
        Ridx = Rrow(1);
 
        Q_left = (A(:, idx(1:Ridx-1)) * diag(P(idx(1 : Ridx-1))) * A(:, idx(1:Ridx-1))');
        for iter_id = 1 : iter_num2
            Q_right = (A(:, idx(Ridx:end)) * diag(P(idx(Ridx:end))) * A(:, idx(Ridx:end))'); 
            RR = Q_left + Q_right;
            for k = 1 : length(Rrow)
                w(idx(Rrow(k)), :) = A(:, idx(Rrow(k)))' * pinv(RR) / (A(:, idx(Rrow(k)))' * pinv(RR) * A(:, idx(Rrow(k))));
                P(idx(Rrow(k))) = w(idx(Rrow(k)), :) * Ru * w(idx(Rrow(k)), :)';
            end
        end
        
    else
        disp(["Check the input mode again!"]);
    end
           
   PoutIAA = P;
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.
[2]J. Yin and T. Chen, "Direction-of-Arrival Estimation Using a Sparse Representation of Array Covariance Vectors," in IEEE Transactions on Signal Processing, vol. 59, no. 9, pp. 4489-4493, Sept. 2011, doi: 10.1109/TSP.2011.2158425.
[3]D. Malioutov, M. Cetin and A. S. Willsky, "A sparse signal reconstruction perspective for source localization with sensor arrays," in IEEE Transactions on Signal Processing, vol. 53, no. 8, pp. 3010-3022, Aug. 2005, doi: 10.1109/TSP.2005.850882.