OMP算法和MP算法的RMSE对比

以下代码是两个算法的均方根误差（RMSE）对比，

一、main函数代码 

%% (SNR-RMSE对比)——MP,OMP
clc
clear all
close all

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

%% 蒙特卡罗实验
monte_num=1000;%蒙特卡罗实验次数
SNR=-10:5:30;
snrnum=length(SNR);
Narg=2;%2种不同的方法
data=zeros(snrnum,monte_num,Narg,K);
for snrindex=1:snrnum
    snr=SNR(snrindex);
    for monindex=1:monte_num
        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 表示使用实际测量的信噪比）
        data(snrindex,monindex,1,:) = MP(y,K,L,a,grid);
        data(snrindex,monindex,2,:) = OMP(y,K,L,a,grid);
    end
end

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

figure()
semilogy(SNR,RMSE_MP,'-mo',LineWidth=2)
hold on
semilogy(SNR,RMSE_OMP,'-gp',LineWidth=2)
grid on
xlabel('SNR（dB）')
ylabel('RMSE')
legend('MP','OMP');

二、MP函数代码 

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

三、OMP函数代码 

function doa_OMP = OMP(y,K,L,a,grid)
%经典OMP算法
%% 初始化
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_OMP=sort(grid(label_set));
end