一.基于压缩感知(CS)的DOA估计方法-OMP-CS算法

阅读须知：

1.本文为本人原创作品仅供学习参考，未经过本人同意禁止转载和抄袭。

2.要想无障碍阅读本文需要一定的压缩感知理论以及压缩感知信号重构算法基础。

3.话不多说，直接开搞。

1 基于压缩感知DOA估计方法原理

        假设有K个远场窄带信号入射到有M个天线的均匀线阵上，第k个信号的入射角度为𝜃𝑘。t时刻阵列接收的单快拍数据矢量可以表示为：

上式是经典的阵列接收数据模型，不管是圆阵还是均匀阵(管它什么牛马阵，都是这个表达式，变的只有阵列流行矩阵A)。其中n(t)表示阵列接收噪声。下面精彩的来了，好好看好好学。对阵列流行矩阵A进行扩展，形成完备的冗余字典G，使它包含了所有可能的方位角度，即 ：

其中，Q为方位角θk的等网格划分数目。我们可以看到，利用冗余字典式G可以将式(1)转化为稀疏表示的问题，即：

 其中，δ为Q维系数向量。即x基于字典G是K-稀疏的。非零元素的位置对应向量的对应角代表了θ的值，非零元素幅值即为信号在采样时刻幅度。观察式(3)可知，x是传感器阵列接收的数据为已知，G为字典矩阵已知，δ是信号的稀疏表示，它是未知的，不为零的元素出现的位置即表征了DOA信息。通过(3)将DOA估计转化为了信号重构的问题，因此只需要求解下式的优化问题，即可重构δ，即：

 考虑到噪声的影响，式(4)可以转化为带有不等式约束的优化问题：

综上所述，基于压缩感知DOA估计的方法归根到底就是一个压缩感知信号重构的问题，关(4)或(5)的信号重构问题有很多算法可以解决，比如我们熟知的贪婪类算法——正交匹配追踪算法(OMP)、广义正交匹配算法(GOMP)、子空间追踪算法(SP)等。接下来在本文中，我将使用OMP算法求解。

2 OMP算法原理

        OMP算法是比较早提出的贪婪类算法。贪婪类算法基于最小化𝑙0范数问题进行稀疏信号构。此类算法因其理论简单、计算量小等优点被广泛应用于实际问题中。该类算法首先从感知矩阵中选取合适的原子，然后将选好的原子组成支撑集，最后利用支撑集来估计信号。因此，该类算法主要包括两部分：原子选取和信号估计。OMP算法利用内积匹配准则进行原子识别，选出与信号残差最匹配的一个原子，用其索引更新支撑集，利用最小二乘法估计信号。

算法具体步骤如下所示：

需要注意的是：上图中所说的观测值为x,即阵列接收的数据。感知矩阵为G,即字典矩阵。稀疏信号x即为δ。

另外：如果对OMP算法仍不熟悉的，可以多看看相关文献，这里为了减少文章冗余，不再赘述。

3 实验结果

        设一个单快拍入射信号，接收阵元为16，入射信号数目为2，信号入射角为[-10°，20°]，得到的仿真结果如下图所示：

        为了探究信噪比与稀疏度对算法的影响，我们分别给出不同信噪比和不同稀疏度(对应目标信源数量)实验图如图2所示：

 根据实验结果得到以下结论：

1.基于压缩感知DOA估计精度与信号重构算法相关，从图中可以看出低信噪比情况下利用OMP算法求解信号精度感人，甚至出现错误估算的情况。因此基于OMP算法的压缩感知DOA估计不适用于低信噪比情况，而在高信噪比时，OMP算法估计精度较好。

2.基于压缩感知DOA估计方法与传统的MUSIC、ESPRIT子空间的方法相比，该方法对信源数(这里对应于稀疏度)的正确估计依赖不强，但需要注意只有在信源估计数量大于信源数量时依赖才不强，反之则可能直接丢失信号。

4 OMP-CS Matlab仿真代码

%Author:Poulen
%Data : 2023.6.12
%Algorithm simulation: OMP-CS Algorithm

%该算法是基于（正交匹配追踪算法）OMP压缩感知信号重构算法的DOA估计方法

clear 
close all;
clc;

%产生信号
%仿真初始值设定
M=16;        %阵元单元
c=3e8;       %光速
f0=77e9;     %初始频率
lambda=c/f0; %波长
slope=30e12; %调频斜率
time=60e-6;  %60us
d=0:lambda/2:(M-1)*lambda/2;%阵列天线
thita=[-10,20];%波达方向
K = length(thita); 

%产生原始信号
N=1;%选取单个信号快拍
t=linspace(0,time,N);

A=zeros(M,K);%导向向量空间 M*K
S=zeros(K,N);%信号空间

St = exp(1j*2*pi*(f0*t(:)+1/2*slope*t(:).^2));
f = 100+f0; %非相干信号的频率
for i = 1:K
    A(:,i) = exp(-1j*2*pi/lambda*d(:)*sind(thita(i)));
%     S(i,:) = St;
    S(i,:) = exp(1j*2*pi*f*t(:));
    f = 1000+f;
end
S = A*S; %产生阵列接收数据

%向数据添加白噪声
SNR = 30; %单位dB
S = S +(randn(size(S)).*std(S))/db2mag(SNR);

%算法步骤：
%%%%%%%%%%%%%%%%
%step 1: 生成冗余矩阵G,使得信号满足稀疏重建的,包含了DOA方向的所有可能的角
%step 2: 信号重构，关于压缩感知信号重构的方法多种多样，在本次demo中选取了正交匹配追踪算法重构稀疏信号
%%%%%%%%%%%%%%%%

%生成过完备冗余字典矩阵，即感知矩阵，其包含所有可能的DOA

scale = -60:59; %所有可能的角度
G = zeros(M,length(scale)); %感知矩阵
output = zeros(length(scale),1);
for i=1:length(scale)
    G(:,i) = exp(-1j*2*pi/lambda*d(:)*sind(scale(i)));
end

%设置初值
K_Sparse = 5; %稀疏度   注意：该稀疏度可随意更改，这里我直接选用为仿真的信源数
r = S;        %残差
R = [];       %原子集
umax_index = [];
%OPM算法信号重构过程
for t=1:K_Sparse
    u = abs(G' * r); 
    temp=find(u==max(u));
    umax_index = [umax_index,temp];
    R = [R,G(:,temp)];
    x = pinv(R)*S ;
    r = S - R * x;
end
output(umax_index) = x; %重构K-sprase信号

figure;
plot(scale,abs(output),'LineWidth',1.6,'Color',[0 0 0]);
xlabel('AngleRange(°)');
ylabel('Amplitude(dB)');
title('OMP-CS Algorithm');
legend([' K-Sparse = ' num2str(K_Sparse) ' ']);

%基于压缩感知DOA估计方法
%在低信噪比是，DOA估计误差较大
%该方法与传统的MUSIC、ESPRITzi子空间的方法相比，该方法对信源数的正确估计依赖不强，注意：只有在信源估计数量
%大于信源数量时依赖才不强，反之则可能直接丢失信号

结束语

      由于实力有限，文中难免出现不足之处。在后续的文章中，我将会更新一些有关DOA估计的相关方法，提供大家借鉴与学习，还望大家多多支持支持。