Mercurial > hg > smallbox

changeset 156:a4d0977d4595 danieleb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First branch commit, danieleb
author danieleb
date Tue, 30 Aug 2011 11:12:31 +0100
parents af307f247ac7
children 00a8473e4b85
files .hgignore DL/two-step DL/SMALL_two_step_DL.m DL/two-step DL/dico_decorr.m data/audio/wav/oboe.mf.c4b4.wav examples/SMALL_test_coherence.m util/SMALL_solve.m
diffstat 6 files changed, 299 insertions(+), 53 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Tue Aug 30 11:12:31 2011 +0100
@@ -0,0 +1,20 @@
+toolboxes/CVX
+toolboxes/GPSR
+toolboxes/KSVD
+toolboxes/KSVDS
+toolboxes/SPARCO
+toolboxes/SparseLab
+toolboxes/Sparsify
+toolboxes/SPGL1
+solvers/SMALL_ompGabor/omp2mexGabor\.mexmaci64
+solvers/SMALL_ompGabor/ompmexGabor\.mexmaci64
+util/ksvd utils/addtocols\.mexmaci64
+util/ksvd utils/col2imstep\.mexmaci64
+util/ksvd utils/collincomb\.mexmaci64
+util/ksvd utils/im2colstep\.mexmaci64
+util/ksvd utils/rowlincomb\.mexmaci64
+util/ksvd utils/sprow\.mexmaci64
+util/Rice Wavelet Toolbox/mdwt\.mexmaci64
+util/Rice Wavelet Toolbox/midwt\.mexmaci64
+util/Rice Wavelet Toolbox/mirdwt\.mexmaci64
+util/Rice Wavelet Toolbox/mrdwt\.mexmaci64
--- a/DL/two-step DL/SMALL_two_step_DL.m	Fri Jul 29 12:35:52 2011 +0100
+++ b/DL/two-step DL/SMALL_two_step_DL.m	Tue Aug 30 11:12:31 2011 +0100
@@ -76,7 +76,7 @@
 end
 % determine if we should do decorrelation in every iteration  %
 
-if isfield(DL.param,'coherence')
+if isfield(DL.param,'coherence') && isscalar(DL.param.coherence)
     decorrelate = 1;
     mu = DL.param.coherence;
 else
@@ -108,13 +108,20 @@
 % main loop %
 
 for i = 1:iternum
+    disp([num2str(i) '/' num2str(iternum)]);
     Problem.A = dico;
     solver = SMALL_solve(Problem, solver);
     [dico, solver.solution] = dico_update(dico, sig, solver.solution, ...
         typeUpdate, flow, learningRate);
-    if (decorrelate)
-        dico = dico_decorr(dico, mu, solver.solution);
-    end
+    dico = normcols(dico);
+        switch DL.param.decFcn
+            case 'mailhe'
+                dico = dico_decorr(dico, mu, solver.solution);
+            case 'tropp'
+                [n m] = size(dico);
+                dico = grassmanian(n,m,[],[],[],dico,true);
+            otherwise
+        end
     
    if ((show_dictionary)&&(mod(i,show_iter)==0))
        dictimg = SMALL_showdict(dico,[8 8],...
@@ -139,4 +146,4 @@
   Y(blockids) = sum(X(:,blockids).^2);
 end
 
-end
\ No newline at end of file
+end
--- a/DL/two-step DL/dico_decorr.m	Fri Jul 29 12:35:52 2011 +0100
+++ b/DL/two-step DL/dico_decorr.m	Tue Aug 30 11:12:31 2011 +0100
@@ -9,6 +9,8 @@
     %   Result:
     %   dico: a dictionary close to the input one with coherence mu.
     
+    eps = 1e-6; % define tolerance for normalisation term alpha
+    
     % compute atom weights
     if nargin > 2
         rank = sum(amp.*amp, 2);
@@ -20,7 +22,7 @@
     % coherence mu. niter can be adjusted to needs.
     niter = 1;
     while niter < 5 && ...
-            max(max(abs(dico'*dico -eye(length(dico))))) > mu + 10^-6
+            max(max(abs(dico'*dico -eye(length(dico))))) > mu + eps
         % find pairs of high correlation atoms
         colors = dico_color(dico, mu);
         
@@ -36,7 +38,7 @@
                 
                 % update the atom
                 corr = dico(:,index(1))'*dico(:,index(2));
-                alpha = sqrt((1-mu*mu)/(1-corr*corr));
+                alpha = sqrt((1-mu*mu)/(1-corr^2+eps));
                 beta = corr*alpha-mu*sign(corr);
                 dico(:,index(2)) = alpha*dico(:,index(2))...
                     -beta*dico(:,index(1));
Binary file data/audio/wav/oboe.mf.c4b4.wav has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/SMALL_test_coherence.m	Tue Aug 30 11:12:31 2011 +0100
@@ -0,0 +1,210 @@
+clear
+
+%% Parameteres
+
+% Dictionary learning parameters
+toolbox   = 'TwoStepDL'; %dictionary learning toolbox
+dicUpdate = {'ksvd','mailhe'}; %dictionary updates
+iterNum   = 20; %number of iterations
+
+% Test signal parameters
+signal    = audio('music03_16kHz.wav'); %audio signal
+blockSize = 256; %size of audio frames
+dictSize  = 512; %number of atoms in the dictionary
+overlap   = 0.5; %overlap between consecutive frames
+sigma     = 1e6; %snr of noise (set to be negligible so that the problem becomes approximation rather than denoising)
+percActiveAtoms = 5; %percentage of active atoms
+
+% Dependent parameters
+nActiveAtoms = fix(blockSize/100*percActiveAtoms); %number of active atoms
+epsilon      = 1/sigma; %error constraint for sparse representation step (corresponds to noise applied to signals)
+minCoherence = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1))); %target coherence (based on coherence lower bound)
+
+% Initial dictionaries
+dctDict = dictionary('dct',blockSize,dictSize);
+dctDict = dctDict.phi;
+gaborParam = struct('N',blockSize,'redundancyFactor',2,'wd',@rectwin);
+gaborDict = Gabor_Dictionary(gaborParam);
+
+%% Generate audio denoising problem with low noise (audio representation)
+SMALL.Problem = generateAudioDenoiseProblem(signal.s,[],blockSize,...
+    dictSize,overlap,sigma); % generate representation problem
+SMALL.Problem.b1 = SMALL.Problem.b; % copy signals from training set b to test set b1 (needed for later functions)
+
+% omp2 sparse representation solver
+ompParam = struct('X',SMALL.Problem.b,'epsilon',epsilon,'maxatoms',nActiveAtoms); %parameters
+solver = SMALL_init_solver('ompbox','omp2',ompParam,false); %solver structure
+
+
+%% Test ksvd dictionary update
+name = dicUpdate{1}; %use ksvd update
+SMALL.DL(1:9) = SMALL_init_DL(toolbox,name); %create dictionary learning structures
+
+% learn with random initialisation and no decorrelation
+SMALL.DL(1).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none'); %parameters for the dictionary learning
+SMALL.DL(1) = SMALL_learn(SMALL.Problem,SMALL.DL(1)); %learn dictionary
+
+% learn with random initialisation and mailhe decorrelation
+SMALL.DL(2).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence); %parameters for the dictionary learning
+SMALL.DL(2) = SMALL_learn(SMALL.Problem,SMALL.DL(2)); %learn dictionary
+
+% learn with random initialisation and tropp decorrelation
+SMALL.DL(3).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence); %parameters for the dictionary learning
+SMALL.DL(3) = SMALL_learn(SMALL.Problem,SMALL.DL(3)); %learn dictionary
+
+% Learn with dct initialisation and no decorrelation
+SMALL.DL(4).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none','initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(4) = SMALL_learn(SMALL.Problem,SMALL.DL(4)); %learn dictionary
+
+% learn with dct initialisation and mailhe decorrelation
+SMALL.DL(5).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(5) = SMALL_learn(SMALL.Problem,SMALL.DL(5)); %learn dictionary
+
+% learn with dct initialisation and tropp decorrelation
+SMALL.DL(6).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(6) = SMALL_learn(SMALL.Problem,SMALL.DL(6)); %learn dictionary
+
+% Learn with gabor initialisation and no decorrelation
+SMALL.DL(7).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none','initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(7) = SMALL_learn(SMALL.Problem,SMALL.DL(7)); %learn dictionary
+
+% learn with gabor initialisation and mailhe decorrelation
+SMALL.DL(8).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(8) = SMALL_learn(SMALL.Problem,SMALL.DL(8)); %learn dictionary
+
+% learn with gabor initialisation and tropp decorrelation
+SMALL.DL(9).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(9) = SMALL_learn(SMALL.Problem,SMALL.DL(9)); %learn dictionary
+
+%% Test mailhe dictionary update
+name = dicUpdate{2}; %use mailhe update
+SMALL.DL(10:18) = SMALL_init_DL(toolbox,name); %create dictionary learning structure
+
+% learn with random initialisation and no decorrelation
+SMALL.DL(10).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none'); %parameters for the dictionary learning
+SMALL.DL(10) = SMALL_learn(SMALL.Problem,SMALL.DL(10)); %learn dictionary
+
+% learn with random initialisation and mailhe decorrelation
+SMALL.DL(11).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence); %parameters for the dictionary learning
+SMALL.DL(11) = SMALL_learn(SMALL.Problem,SMALL.DL(11)); %learn dictionary
+
+% learn with random initialisation and tropp decorrelation
+SMALL.DL(12).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence); %parameters for the dictionary learning
+SMALL.DL(12) = SMALL_learn(SMALL.Problem,SMALL.DL(12)); %learn dictionary
+
+% Learn with dct initialisation and no decorrelation
+SMALL.DL(13).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none','initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(13) = SMALL_learn(SMALL.Problem,SMALL.DL(13)); %learn dictionary
+
+% learn with dct initialisation and mailhe decorrelation
+SMALL.DL(14).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(14) = SMALL_learn(SMALL.Problem,SMALL.DL(14)); %learn dictionary
+
+% learn with dct initialisation and tropp decorrelation
+SMALL.DL(15).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence,'initdict',dctDict); %parameters for the dictionary learning
+SMALL.DL(15) = SMALL_learn(SMALL.Problem,SMALL.DL(15)); %learn dictionary
+
+% Learn with gabor initialisation and no decorrelation
+SMALL.DL(16).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','none','initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(16) = SMALL_learn(SMALL.Problem,SMALL.DL(16)); %learn dictionary
+
+% learn with gabor initialisation and mailhe decorrelation
+SMALL.DL(17).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','mailhe','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(17) = SMALL_learn(SMALL.Problem,SMALL.DL(17)); %learn dictionary
+
+% learn with gabor initialisation and tropp decorrelation
+SMALL.DL(18).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+    'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+    'decFcn','tropp','coherence',minCoherence,'initdict',gaborDict); %parameters for the dictionary learning
+SMALL.DL(18) = SMALL_learn(SMALL.Problem,SMALL.DL(18)); %learn dictionary
+
+%% Evaluate coherence and snr of representation for the various methods
+sigNoiseRatio = zeros(18,1);
+mu = zeros(18,1);
+for i=1:18
+    SMALL.Problem.A = SMALL.DL(i).D;
+    tempSolver = SMALL_solve(SMALL.Problem,solver);
+    sigNoiseRatio(i) = snr(SMALL.Problem.b,SMALL.DL(i).D*tempSolver.solution);
+    dic(i) = dictionary(SMALL.DL(i).D);
+    mu(i) = dic(i).coherence;
+end
+
+%% Plot results
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));
+maxSNR = max(sigNoiseRatio);
+
+figure, subplot(2,2,1)
+snrMat = buffer(sigNoiseRatio(1:9),3);
+bar(snrMat');
+title('SNR - KSVD Update')
+xlabel('Initial dictionary')
+ylabel('SNR (dB)')
+set(gca,'XTickLabel',{'data','dct','gabor'},'YLim',[0 maxSNR+1]);
+legend('none','Mailhe','Tropp')
+grid on
+
+subplot(2,2,2), grid on
+snrMat = buffer(sigNoiseRatio(10:18),3);
+bar(snrMat');
+title('SNR - Mailhe Update')
+xlabel('Initial dictionary')
+ylabel('SNR (dB)')
+set(gca,'XTickLabel',{'data','dct','gabor'},'YLim',[0 maxSNR+1]);
+legend('none','Mailhe','Tropp')
+grid on
+
+subplot(2,2,3), hold on, grid on
+title('Coherence - KSVD Update')
+muMat = buffer(mu(1:9),3);
+line([0.5 3.5],[1 1],'Color','r');
+bar(muMat');
+line([0.5 3.5],[minMu minMu],'Color','k');
+set(gca,'XTick',1:3,'XTickLabel',{'data','dct','gabor'},'YLim',[0 1.05])
+legend('\mu_{max}','none','Mailhe','Tropp','\mu_{min}')
+ylabel('\mu')
+xlabel('Initial dictionary')
+
+subplot(2,2,4), hold on, grid on
+title('Coherence - Mailhe Update')
+muMat = buffer(mu(10:18),3);
+line([0.5 3.5],[1 1],'Color','r');
+bar(muMat');
+line([0.5 3.5],[minMu minMu],'Color','k');
+set(gca,'XTick',1:3,'XTickLabel',{'data','dct','gabor'},'YLim',[0 1.05])
+legend('\mu_{max}','none','Mailhe','Tropp','\mu_{min}')
+ylabel('\mu')
+xlabel('Initial dictionary')
--- a/util/SMALL_solve.m	Fri Jul 29 12:35:52 2011 +0100
+++ b/util/SMALL_solve.m	Tue Aug 30 11:12:31 2011 +0100
@@ -15,7 +15,7 @@
 %   published by the Free Software Foundation; either version 2 of the
 %   License, or (at your option) any later version.  See the file
 %   COPYING included with this distribution for more information.
-%   
+%
 %%
 
 if isa(Problem.A,'float')
@@ -45,51 +45,58 @@
 
 start=cputime;
 tStart=tic;
-if strcmpi(solver.toolbox,'sparselab')
-    y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']);
-elseif strcmpi(solver.toolbox,'sparsify')
-    y = eval([solver.name,'(b, A, n, ''P_trans'', AT,',solver.param,');']);
-elseif (strcmpi(solver.toolbox,'spgl1')||strcmpi(solver.toolbox,'gpsr'))
-    y = eval([solver.name,'(b, A,',solver.param,');']);
-elseif (strcmpi(solver.toolbox,'SPAMS'))
-    y = eval([solver.name,'(b, A, solver.param);']);
-elseif (strcmpi(solver.toolbox,'SMALL'))
-    if isa(Problem.A,'float')
-        y = eval([solver.name,'(A, b, n,',solver.param,');']);
-    else
-        y = eval([solver.name,'(A, b, n,',solver.param,',AT);']);
-    end
-elseif (strcmpi(solver.toolbox, 'ompbox'))
-    G=A'*A;
-    epsilon=solver.param.epsilon;
-    maxatoms=solver.param.maxatoms;
-    y = eval([solver.name,'(A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']);
-elseif (strcmpi(solver.toolbox, 'ompsbox'))
-    basedict = Problem.basedict;
-     if issparse(Problem.A)
-        A = Problem.A;
-      else
-        A = sparse(Problem.A);
-      end
-    G = dicttsep(basedict,A,dictsep(basedict,A,speye(size(A,2))));
-    epsilon=solver.param.epsilon;
-    maxatoms=solver.param.maxatoms;
-    y = eval([solver.name,'(basedict, A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']);
-    Problem.sparse=1;
-%   To introduce new sparse representation algorithm put the files in
-%   your Matlab path. Next, unique name <TolboxID> for your toolbox and
-%   prefferd API <Preffered_API> needs to be defined.
-%
-% elseif strcmpi(solver.toolbox,'<ToolboxID>')
-%
-%     % - Evaluate the function (solver.name - defined in the main) with
-%     %   parameters given above
-%
-%     y = eval([solver.name,'(<Preffered_API>);']);
-
-else
-    printf('\nToolbox has not been registered. Please change SMALL_solver file.\n');
-    return
+switch solver.toolbox
+    case 'sparselab'
+        y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']);
+    case 'sparsify'
+        y = eval([solver.name,'(b, A, n, ''P_trans'', AT,',solver.param,');']);
+    case 'spgl1'
+        y = eval([solver.name,'(b, A,',solver.param,');']);
+    case 'gpsr'
+        y = eval([solver.name,'(b, A,',solver.param,');']);
+    case 'SPAMS'
+        y = eval([solver.name,'(b, A, solver.param);']);
+    case 'SMALL'
+        if isa(Problem.A,'float')
+            y = eval([solver.name,'(A, b, n,',solver.param,');']);
+        else
+            y = eval([solver.name,'(A, b, n,',solver.param,',AT);']);
+        end
+    case 'ompbox'
+        G = A'*A;
+        maxatoms=solver.param.maxatoms;
+        switch solver.name
+            case 'omp'
+                y = omp(A,b,G,maxatoms,'checkdict','off');
+            case 'omp2'
+                epsilon=solver.param.epsilon;
+                y = omp2(A,b,G,epsilon,'maxatoms',maxatoms,'checkdict','off');
+        end
+    case 'ompsbox'
+        basedict = Problem.basedict;
+        if issparse(Problem.A)
+            A = Problem.A;
+        else
+            A = sparse(Problem.A);
+        end
+        G = dicttsep(basedict,A,dictsep(basedict,A,speye(size(A,2))));
+        epsilon=solver.param.epsilon;
+        maxatoms=solver.param.maxatoms;
+        y = eval([solver.name,'(basedict, A, b, G,epsilon,''maxatoms'',maxatoms,''checkdict'',''off'');']);
+        Problem.sparse=1;
+        %   To introduce new sparse representation algorithm put the files in
+        %   your Matlab path. Next, unique name <TolboxID> for your toolbox and
+        %   prefferd API <Preffered_API> needs to be defined.
+        %
+        % elseif strcmpi(solver.toolbox,'<ToolboxID>')
+        %
+        %     % - Evaluate the function (solver.name - defined in the main) with
+        %     %   parameters given above
+        %
+        %     y = eval([solver.name,'(<Preffered_API>);']);
+    otherwise
+        printf('\nToolbox has not been registered. Please change SMALL_solver file.\n');
+        return
 end
 
 %%