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