对于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 3for t = 1:s% 初始化原子投影proj=[];% 4for j=1:Nproj = [proj,norm(r'*PHI(:,j),1)];end% 选取最大投影对应的原子索引[~, index] = max(proj);Omega = [Omega, index];% 重构X的每行元素for k=1:KX_hat(Omega,k)=pinv(PHI(:,Omega))*Y(:,k);end%更新测量向量投影Y_tplus1 = PHI(:,Omega)*pinv(PHI(:,Omega))*Y;%更新残差r = Y - Y_tplus1;endend

本代码复现自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 boundx = zeros(N,L); % Initializing sparse vector to be recovereduset = randperm(N,k); x(uset,:) = sqrt(0.5)*(randn(k,L) + 1i*randn(k,L)); % Sparse vector initializednoise = sqrt(1/2)*(randn(M,L) + 1i*randn(M,L)); % zero mean, unit covariance complex noise vectorpower=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 boundx = zeros(N,L); % Initializing sparse vector to be recovereduset = randperm(N,sp); x(uset,:) = sqrt(0.5)*(randn(sp,L) + 1i*randn(sp,L)); % Sparse vector initializednoise = sqrt(1/2)*(randn(M,L) + 1i*randn(M,L)); % zero mean, unit covariance complex noise vectorpower=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，效果就不是很好了。