对于SOMP算法的测试

刚开始只上传了SOMP算法的代码，并没有过多介绍。

所以本篇文章对SOMP算法用法进行一个介绍

SOMP算法代码

function [X_hat] = MMV_SOMP(Y, PHI, s)
    % SOMP:同时正交匹配追踪 simultaneous orthogonal matching pursuit
    %         论文：J. Determe, J. Louveaux, L. Jacques and
    %         F. Horlin, \On the Noise Robustness of Simultaneous Orthogonal Matching Pursuit,"
    %         in IEEE Transactions on Signal Processing, vol. 65, no. 4, pp. 864-875, Feb. 15
    %         2017. doi: 10.1109/TSP.2016.2626244. http://ieeexplore.ieee.org/document/7738592/
    % write:zts
    % date：23,11,17
    % version:1.0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 参数：
    % Y: 观测矩阵，每列是一个观测向量
    % PHI: 字典矩阵
    % s: 稀疏度
    
    % 获取矩阵大小 
  
    [~, N] = size(PHI);
    [~, K] = size(Y);
	% 1 初始化残差和索引
    r = Y;
	% 初始化索引集
    Omega = [];
    % 初始化重构矩阵
    X_hat=zeros(N,K);
	% 2 3
    for t = 1:s
	
		% 初始化原子投影
		proj=[];
		% 4
		for j=1:N
			proj = [proj,norm(r'*PHI(:,j),1)];
		end
        % 选取最大投影对应的原子索引
        [~, index] = max(proj);
        Omega = [Omega, index];
        % 重构X的每行元素
        for k=1:K
             X_hat(Omega,k)=pinv(PHI(:,Omega))*Y(:,k);
        end
		%更新测量向量投影
		Y_tplus1 = PHI(:,Omega)*pinv(PHI(:,Omega))*Y;
        %更新残差
        r = Y - Y_tplus1;
 
    end

end

本代码复现自SOMP论文。

SOMP算法解决的是具有公共稀疏集的压缩感知内稀疏信号恢复问题。

要求信号是行稀疏。最初的代码是要求已知稀疏度的，后面在各个场景里应用的代码是通过门限来判断的，一般的门限适合噪声方差有关的。

测试信噪比-NMSE

如果观测矩阵不是独立同分布的，可能算法就不会工作，当然很多算法对这点的要求都很高。

第一个先测试信噪比对算法影响

clc;
clear all;
N = 128; % length of vector to be recovered
M = 64; % number of measurement
L=16;% number of layer
k =13; % Sparsity level
niter = 50; % number of iteration
SNR_dB=0:5:25;
NMSE_dB = zeros(1,length(SNR_dB));
for n=1:length(SNR_dB)
    SNR_dB_now = SNR_dB(n);
    A = sqrt(1/2)*(randn(M,N) + 1i*randn(M,N)); % Sensing matrix construction for theroetical bound
    x = zeros(N,L); % Initializing sparse vector to be recovered
    
    uset = randperm(N,k); 
    x(uset,:) = sqrt(0.5)*(randn(k,L) + 1i*randn(k,L)); % Sparse vector initialized
    noise = sqrt(1/2)*(randn(M,L) + 1i*randn(M,L)); % zero mean, unit covariance complex noise vector
    power=sum(abs(A*x).^2,'all')/M/L;
    var = power/(10^(0.1*SNR_dB_now));
    noise = sqrt(var)*noise;
    y = A*x + noise; % create measurement
    %% 算法测试
    [X_hat] = MMV_SOMP(y, A, k);
     NMSE_dB(n) = 10*log10( norm(X_hat-x)^2/(norm(x)^2));
end
plot(SNR_dB,NMSE_dB,'r-o');
legend('SOMP');
xlabel('SNR\_dB'),ylabel('NMSE\_dB');
title('测试');

可以看出性能还不错，缺点就是复杂度过高。当然如果没有这个公共支撑集，我们看执行L次的单测量OMP。

OMP的算法在之前的博客里贴的也有，当然问GPT也能问出来，OMP应该是最容易复现的代码了，就不再次贴出来了。

可以看出OMP没有利用公共稀疏集，所以性能有些损失。当然这是单次实验，可能偶尔有重合的机会。

测试稀疏度-NMSE

clc;
clear all;
N = 128; % length of vector to be recovered
M = 64; % number of measurement
L=16;% number of layer
k =5:5:40; % Sparsity level
niter = 50; % number of iteration
SNR_dB=20;
NMSE_dB_SOMP = zeros(1,length(k));
NMSE_dB_OMP = zeros(1,length(k));
for n=1:length(k)
    sp = k(n);
    A = sqrt(1/2)*(randn(M,N) + 1i*randn(M,N)); % Sensing matrix construction for theroetical bound
    x = zeros(N,L); % Initializing sparse vector to be recovered
    
    uset = randperm(N,sp); 
    x(uset,:) = sqrt(0.5)*(randn(sp,L) + 1i*randn(sp,L)); % Sparse vector initialized
    noise = sqrt(1/2)*(randn(M,L) + 1i*randn(M,L)); % zero mean, unit covariance complex noise vector
    power=sum(abs(A*x).^2,'all')/M/L;
    var = power/(10^(0.1*SNR_dB));
    noise = sqrt(var)*noise;
    y = A*x + noise; % create measurement
    %% 算法测试
   [X_hat] = MMV_SOMP(y, A, sp);
   NMSE_dB_SOMP(n) = 10*log10( norm(X_hat-x)^2/(norm(x)^2));
    x_omp=	MMV_OMP(y,A,L,sp);
    NMSE_dB_OMP(n) = 10*log10( norm(x_omp-x)^2/(norm(x)^2));
end
plot(k,NMSE_dB_SOMP,'r-o',k,NMSE_dB_OMP,'g-+');
legend('SOMP','OMP');
xlabel('稀疏度'),ylabel('NMSE\_dB');
title('20dB测试');

可以看出来SOMP还是比较平滑的，最大稀疏度是40/128=0.3125.而到了OMP，效果就不是很好了。