OMP算法估计DOA

%% 经典OMP算法
clc
clear all
close all

%% 阵列参数
M = 10;%%阵元数
theta = [20 30 51];%%入射信源角度
K = length(theta);%信源数
d = 0.5;%阵元间距半波长
dd = 0:M-1;%阵列流形序号
SNR = 20;%信噪比

%% 构造接收信号
grid = -90:90;
L = length(grid);
x = zeros(L,1);
x(theta+91) = 1;
a = exp(1i*2*pi*d*dd'*sind(grid));%阵列流形矩阵
y = a*x;%信号接收矩阵
y = awgn(y,SNR,'measured');%添加高斯白噪声 （measured 表示使用实际测量的信噪比）

%% 初始化
r(:,1) = y;%初始残差
label_set = [];%初始标签集
%% 迭代
iter = K;%迭代次数
for j = 1:iter
    for i = 1:L
        correlated_atom(i) = abs(r(:,j)'*a(:,i));
    end
    [~,index] = max(abs(correlated_atom));%辨识
    label_set = [label_set,index];
    x_est{:,j} = pinv(a(:,label_set)'*a(:,label_set))*a(:,label_set)'*y;%估计
    r(:,j+1) = y-a(:,label_set)*x_est{:,j};%更新残差 
end
doa = sort(grid(label_set));
disp(doa)