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

    function [Sest,d] =OMPmod3(Phi,y,strue,K);

% Othogonal matching pursuit of Tropp et Gilbert
% modified  : finds at each iteration the vector
% that maximizes the loss in the residual
% and also makes sure the selected vectors have
% no linear dependencies
% Implementation probably computationnally very inefficient
% Finds greedy approx of s s.t. phi*s = y
% y : data
% phi : measurement matrix
% K :number of vectors used in the approx (typ. = S or M)
% 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
%
Sest=zeros(size(strue));
Pos=[]; % positions indexes of components of s
rtm1=y; % first residual
Phit=[]; % Matrix of the colums used to represent y 
t=1; %iteration number
d=0; % error flag
Mat=Phi;
    for t=1:K;
        for v=1:length(Mat(1,:))
            Phit2=[Phit Mat(:,v)]; % takes every new vector
            
            if rank(Phit2)~=length(Phit2(1,:))
                score(v)=1e20;
            else
                sest=Phit2\y;r=y-Phit2*sest;
            score(v)=norm(r,2);
            end
        end     
        [i,j]=min(score);
         Phit=[Phit Mat(:,j)]; % add selected vector to Phit
         % remove now the selected vector to Mat
       if j==1
    Mat=Mat(:,2:length(Mat(1,:)));
end
if j==length(Mat(1,:))
    Mat=Mat(:,1:length(Mat(1,:))-1);
end
if (j~=1) & (j~=(-1))
    Mat=[Mat(:,1:j-1) Mat(:,j+1:length(Mat(1,:)))];
end
    sest=Phit\y;
    rtm1=y-Phit*sest;
    % finds what is the position in Phi of the jth vector in Mat
    if t==1
     Pos=[Pos j];
    else
        l=length(find(Pos<=j));t=j+l;
        while length(find(Pos<=t))~=l
            t=t+1;l=l+1;
        end
        Pos=[Pos t];
    end
    Sest(Pos)=sest;
error=norm(abs(Sest-strue),'inf');
if error < 1e-3
    d=1;
    disp('modified OMP : success')
    break
end
%t=length(Phi(1,:));
end;
if d==0
    disp('OMP mod : failed')
end