view solvers/SMALL_MP.m @ 9:28f2b5fe3483

(none)

function [A]=SMALL_MP(Dict,X, m, maxNumCoef, errorGoal, varargin) 
%%
%=============================================
% Sparse coding of a group of signals based on a given 
% dictionary and specified number of atoms to use. 
%  input arguments: Dict - the dictionary
%                   X - the signals to represent
%                   m - number of atoms in Dictionary
%                   errorGoal - the maximal allowed representation error for
%                   each signal.
%
% optional:         if Dict is function handle then Transpose Dictionary
%                   handle needs to be specified.
%
% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
% with pinv() is changed with matching pursuit)
% Ivan Damnjanovic 2009
%=============================================
%%
%%
% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers 
explicitD=0;
if     isa(Dict,'float')      
            D =@(z) Dict*z;  
            Dt =@(z) Dict'*z;
            explicitD=1;
elseif isobject(Dict) 
            D =@(z) Dict*z;  
            Dt =@(z) Dict'*z;
elseif isa(Dict,'function_handle') 
    try
        DictT=varargin{1};
        if isa(DictT,'function_handle'); 
                D=Dict; 
                Dt=DictT;
        else
                error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
        end
    catch
        error('If Dictionary is a function handle, Transpose Dictionary  needs to be specified. Exiting.');
    end
else
    error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
end
%%
[n,P]=size(X);

%E2 = errorGoal^2;


A = sparse(m,size(X,2));
for k=1:1:P,
   
    x = X(:,k);
    residual=x;
	indx = [];
	j=0;
    
	currResNorm = norm(residual);
    errGoal=errorGoal*currResNorm;
	
    a = zeros(m,1);


    while currResNorm>errGoal && j < maxNumCoef,
		
        j = j+1;
        dir=Dt(residual);
        
        [tmp__, pos]=max(abs(dir));
        indx(j)=pos;
        a(pos)=dir(pos);
        e = zeros(m,1);
        e(pos) = 1;
        residual=residual-D(e)*a(pos);
    
		currResNorm = norm(residual);
   end;
   if (~isempty(indx))
       A(indx,k)=a(indx);
   end
end;
return;