2023年7月29日发(作者：)

压缩感知重构算法——SP算法SP（subspace pursuit）算法是压缩感知中⼀种⾮常重要的贪婪算法，它有较快的计算速度和较好的重构概率，在实际中应⽤较多。本⽂给出了SP算法的matlab代码，以及相应的测试函数。参考⽂献：Dai W, Milenkovic O. Subspace pursuit for compressive sensing signal reconstruction[J]. Information Theory, IEEETransactions on, 2009, 55(5): 2230-2249.⽂献下载地址：Matlab代码：SP_paper.m

function x=SP_paper(Phi,y,K)%SP算法%获取Phi矩阵的⾏数和列数[M,N]=size(Phi);%初始化步骤%将Phi的每列与y做相关，得到⼀个N*1的矩阵(列向量)correlation=Phi'*y;%对correlation取绝对值后排序，按从⼤到⼩的顺序[var,pos] = sort(abs(correlation),'descend');%声明⼀个空集T，⽤于记录Phi的列数标值T=[];T=union(T,pos(1:K));y_r=resid_paper(y,Phi(:,T));%迭代%使⽤如下形式的do---while结构% while(1)% if(condition)% break;% end% endcount=1;while(1)%根据残差计算待增加的列数，得到T_addcorrelation=Phi'*y_r;[var,pos] = sort(abs(correlation),'descend');T_add=union([],pos(1:K));%合并已有的T和T_addT=union(T,T_add);%x_p=((Phi(:,T)'*Phi(:,T))eye(length(T)))*Phi(:,T)'*y;%proj_paper(y,Phi(:,T));%更新下标记录T[var,pos] = sort(abs(x_p),'descend');%取前K个最⼤值T=union([],T(pos(1:K)));%计算新的残差y_r_n=resid_paper(y,Phi(:,T));%判断是否退出循环,且置为最⼤迭代100次if(norm(y_r_n)>=norm(y_r) || count>100)break;end%若不退出循环,进⾏新⼀轮的迭代y_r=y_r_n;count=count+1;end%退出循环后，做最后的数据输出x=zeros(N,1);x(T)=((Phi(:,T)'*Phi(:,T))eye(length(T)))*Phi(:,T)'*y;end

function y_r=resid_paper(y,Phi)%计算y在Phi上的投影残差%获取矩阵Phi的⾏数和列数,M没有⽤[M,N]=size(Phi);%判断矩阵(Phi'*Phi)是否可逆if(rank(Phi'*Phi)~=N)error('矩阵不可逆');endy_p=Phi*((Phi'*Phi)eye(N))*Phi'*y;y_r=y-y_p;end

%%%%%%%%%%%%%%%%%%%%%%%%ion [y,Phi,x]=dataGen(M,N,K)% 产⽣贪婪算法所需要的数据%⽣成-1/+1的原始信号xx = zeros(N,1);q = randperm(N); %y=randperm(n),是把1到n这些数随机打乱得到的⼀个数字序列。x(q(1:K)) = sign(randn(K,1));%⽣成测量矩阵PhiPhi = randn(M,N);Phi = orth(Phi')'; %求矩阵正交基,B = orth(A),B的列向量是正交向量%⽣成观测向量y = Phi*x;end

%%%%%%%%%%%%%%%%%%%%%%%%SP_;clc;%信号长度N = 256;%信号满⾜K-稀疏K = 30;%观测向量的长度M = 128;[y,Phi,x] = dataGen(M,N,K);xp = SP_paper(Phi,y,K);% [xp,support,iterCount]=SP_res(Phi,y,K);figuresubplot(3,1,1),stem(x,'.'),title('原始'),axis([0,N,-1,1])subplot(3,1,2),stem(xp,'.'),title('重构'),axis([0,N,-1,1])subplot(3,1,3),stem(abs(xp-x),'.'),title('残差'),axis([0,N,0,max(abs(xp-x))])