压缩感知重构算法——SP算法

SP（subspace pursuit）算法是压缩感知中一种非常重要的贪婪算法，它有较快的计算速度和较好的重构概率，在实际中应用较多。本文给出了SP算法的matlab代码，以及相应的测试函数。

参考文献：Dai W, Milenkovic O. Subspace pursuit for compressive sensing signal reconstruction[J]. Information Theory, IEEE Transactions on, 2009, 55(5): 2230-2249.

文献下载地址：http://arxiv.org/pdf/0803.0811.pdf

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_add    correlation=Phi'*y_r;    [var,pos] = sort(abs(correlation),'descend');    T_add=union([],pos(1:K));        %合并已有的T和T_add    T=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;endfunction 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

测试代码：

dataGen.m

function [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_main.m

clear;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))])