comparison solvers/SMALL_MP.m @ 1:7750624e0c73 version0.5

(none)
author idamnjanovic
date Thu, 05 Nov 2009 16:36:01 +0000
parents
children 257c289bf11d
comparison
equal deleted inserted replaced
0:5181bee80bc1 1:7750624e0c73
1 function [A]=SMALL_MP(Dict,X, m, maxNumCoef, errorGoal, varargin)
2 %%
3 %=============================================
4 % Sparse coding of a group of signals based on a given
5 % dictionary and specified number of atoms to use.
6 % input arguments: Dict - the dictionary
7 % X - the signals to represent
8 % m - number of atoms in Dictionary
9 % errorGoal - the maximal allowed representation error for
10 % each signal.
11 %
12 % optional: if Dict is function handle then Transpose Dictionary
13 % handle needs to be specified.
14 %
15 % based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
16 % with pinv() is changed with matching pursuit)
17 % Ivan Damnjanovic 2009
18 %=============================================
19 %%
20 %%
21 % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers
22 explicitD=0;
23 if isa(Dict,'float')
24 D =@(z) Dict*z;
25 Dt =@(z) Dict'*z;
26 explicitD=1;
27 elseif isobject(Dict)
28 D =@(z) Dict*z;
29 Dt =@(z) Dict'*z;
30 elseif isa(Dict,'function_handle')
31 try
32 DictT=varargin{1};
33 if isa(DictT,'function_handle');
34 D=Dict;
35 Dt=DictT;
36 else
37 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
38 end
39 catch
40 error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.');
41 end
42 else
43 error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
44 end
45 %%
46 [n,P]=size(X);
47
48 E2 = errorGoal^2;
49
50
51 A = sparse(m,size(X,2));
52 for k=1:1:P,
53
54 x = X(:,k);
55 residual=x;
56 indx = [];
57 j=0;
58
59 currResNorm = norm(residual);
60 errorGoal=errorGoal*currResNorm;
61
62 a = zeros(m,1);
63
64
65 while currResNorm>errorGoal & j < maxNumCoef,
66
67 j = j+1;
68 dir=Dt(residual);
69
70 [tmp__, pos]=max(abs(dir));
71 indx(j)=pos;
72 a(pos)=dir(pos);
73 e = zeros(m,1);
74 e(pos) = 1;
75 residual=residual-D(e)*a(pos);
76
77 currResNorm = norm(residual);
78 end;
79 if (~isempty(indx))
80 A(indx,k)=a(indx);
81 end
82 end;
83 return;
84
85