多快拍下的OMP算法均方误差分析代码

%% OMPMMV（多快拍下的OMP）算法
clc
clear all
close all

%% 阵列参数
M = 10;%%阵元数
N = 10;%快拍数
theta = [20 40];%%入射信源角度
K = length(theta);%信源数
d = 0.5;%阵元间距半波长
dd = 0:M-1;%阵列流形序号
Grid = -90:90;
L = length(Grid);

%% 构造接收信号
monte_num=200;%蒙特卡罗实验次数
SNR=-10:5:30;
SNR_num=length(SNR);
data=zeros(SNR_num,monte_num,K);
for snrindex=1:SNR_num
    snr=SNR(snrindex);
    for monindex=1:monte_num
        x = zeros(L,N);
        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) = norm(a(:,i)'*r(:,:,j),2);
            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));
        data(snrindex,monindex,:) = doa;
    end
end

%% RMSE均方误差
MSE=zeros(SNR_num,K);
%k建立几个维度，优先级从高到低依次为：doa→重复蒙特卡罗实验→SNR
for snrindex=1:SNR_num
    for monindex=1:monte_num
        for doaindex=1:K
            MSE(snrindex,doaindex)=MSE(snrindex,doaindex)+(data(snrindex,monindex,doaindex)-theta(doaindex))^2;
        end
    end
end
%%取各个信源估计误差的平均值
MSE_OMPMMV=MSE(:,1)+MSE(:,2);
RMSE_OMPMMV=sqrt(MSE_OMPMMV/(monte_num*K));

%% 绘图
figure()
semilogy(SNR,RMSE_OMPMMV,'-r^',LineWidth=2)
grid on
xlabel('SNR（dB）')
ylabel('RMSE')
legend('OMPMMV');