www.gusucode.com > 压缩感知的几种重构算法比较的源代码 > 压缩感知的几种重构算法比较的源代码/MMA14/SPA.m

    function [Sest,d]=SPA(Phi,K,y,utrue)

% Subspace Pursuit Algorithm of Wai
% K : sparsity of Sest
% Phi : obs matrix
% y: obs vector
% Written by David Mary
% This script/program is released under the Commons Creative Licence
% with Attribution Non-commercial Share Alike (by-nc-sa)
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% Disclaimer: the short answer, this script is for educational purpose only.
% More Disclaimer:  http://igorcarron.googlepages.com/disclaimer

% Init.
[k,z]=sort(abs(Phi'*y));k=flipud(k);z=flipud(z);%把Phi'*y的值先升序排列 再按列上下翻转 Z是位置集
T=z(1:K);phit=Phi(:,T);
yr=y-phit*pinv(phit)*y;

 t=1; 
 while t < 101    
[k2,z2]=sort(abs(Phi'*yr));k2=flipud(k2);z2=flipud(z2); % 矩阵的上下翻转
T2=z2(1:K);
T3=sort(union(T,T2));phit3=Phi(:,T3);
xp=(pinv(phit3))*y;
[k3,z3]=sort(abs(xp));k3=flipud(k3);z3=flipud(z3);
T=T3(z3(1:K));phit4=Phi(:,T);
yr=y-phit4*(pinv(phit4))*y;%pinv 伪逆
S=(pinv(phit4))*y;Sest=zeros(size(utrue));Sest(T)=abs(S);
n=norm(yr,2);
d=0;n2=norm(abs(Sest-utrue),'inf');%NORM(X,inf) is the infinity norm of X, the largest row sum,= max(sum(abs(X'))).
if n2 <1e-3
    d=1;t=1e10;
    disp('SPA: success');
end
 t=t+1;
 end
 if d==0
     disp('SPA: failed')
 end