changeset 224:fd0b5d36f6ad danieleb

Updated the contents of this branch with the contents of the default branch.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Thu Apr 12 13:52:28 2012 +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/.hgtags	Wed Apr 11 16:28:46 2012 +0100
+++ b/.hgtags	Thu Apr 12 13:52:28 2012 +0100
@@ -3,5 +3,6 @@
 9ff69e8e049f936804d0e5876cd4d367be9f3c4a backup 14032011
 b14e1f6ee4bea90a8a894a52c1114a72aa818071 ver_1.1
 19e0af5709140e163faaf9d8cf4b83a664be1edc release_1.5
-30872eb3d1606b89f8d071cfdb7cddb874fab8fd release_1.51
-4205744092e6a16951e48acef27fb41e32851758 release_1.9
+a4d0977d45956b96579a3ed83a4e6a1869ee6055 danieleb
+a4d0977d45956b96579a3ed83a4e6a1869ee6055 danieleb
+0000000000000000000000000000000000000000 danieleb
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/dl_ramirez.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,189 @@
+function DL = dl_ramirez(Problem,DL)
+%% Dictionary learning with incoherent dictionary
+%
+% REFERENCE
+% I. Ramirez, F. Lecumberry and G. Sapiro, Sparse modeling with universal
+% priors and learned incoherent dictionaries.
+
+%%
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2011 Daniele Barchiesi.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   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.
+
+%% Test function
+if ~nargin, testdl_ramirez; return; end
+
+%% Parameters & Defaults
+X = Problem.b;					%matrix of observed signals
+
+% determine dictionary size %
+if (isfield(DL.param,'initdict')) %if the dictionary has been initialised
+	if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
+		dictSize = length(DL.param.initdict);
+	else
+		dictSize = size(DL.param.initdict,2);
+	end
+end
+if (isfield(DL.param,'dictsize'))
+	dictSize = DL.param.dictsize;
+end
+
+if (size(X,2) < dictSize)
+	error('Number of training signals is smaller than number of atoms to train');
+end
+
+
+% initialize the dictionary %
+if (isfield(DL.param,'initdict')) && ~isempty(DL.param.initdict);
+	if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
+		D = X(:,DL.param.initdict(1:dictSize));
+	else
+		if (size(DL.param.initdict,1)~=size(X,1) || size(DL.param.initdict,2)<dictSize)
+			error('Invalid initial dictionary');
+		end
+		D = DL.param.initdict(:,1:dictSize);
+	end
+else
+	data_ids = find(colnorms_squared(X) > 1e-6);   % ensure no zero data elements are chosen
+	perm = randperm(length(data_ids));
+	D = X(:,data_ids(perm(1:dictSize)));
+end
+
+
+% coherence penalty factor
+if isfield(DL.param,'zeta')
+	zeta = DL.param.zeta;
+else
+	zeta = 0.1; 
+end
+
+% atoms norm penalty factor
+if isfield(DL.param,'eta')
+	eta = DL.param.eta;
+else
+	eta = 0.1; 
+end
+
+% number of iterations (default is 40) %
+if isfield(DL.param,'iternum')
+    iternum = DL.param.iternum;
+else
+    iternum = 40;
+end
+
+% show dictonary every specified number of iterations
+if isfield(DL.param,'show_dict')
+    show_dictionary=1;
+    show_iter=DL.param.show_dict;
+else
+    show_dictionary=0;
+    show_iter=0;
+end
+
+tmpTraining = Problem.b1;
+Problem.b1 = X;
+if isfield(Problem,'reconstruct')
+    Problem = rmfield(Problem, 'reconstruct');
+end
+
+
+%% Main Algorithm
+Dprev = D;						%initial dictionary
+Aprev = D\X;				    %set initial solution as pseudoinverse
+for i = 1:iternum
+	%Sparse Coding by 
+	A = sparsecoding(X,D,Aprev);
+	%Dictionary Update
+	D = dictionaryupdate(X,A,Dprev,zeta,eta);
+	
+	Dprev = D;
+	Aprev = A;
+   if ((show_dictionary)&&(mod(i,show_iter)==0))
+       dictimg = SMALL_showdict(dico,[8 8],...
+            round(sqrt(size(dico,2))),round(sqrt(size(dico,2))),'lines','highcontrast');  
+       figure(2); imagesc(dictimg);colormap(gray);axis off; axis image;
+       pause(0.02);
+   end
+end
+
+Problem.b1 = tmpTraining;
+DL.D = D;
+
+end
+
+function A = sparsecoding(X,D,Aprev)
+%Sparse coding using a mixture of laplacians (MOL) as universal prior.
+
+%parameters
+K = size(D,2);	%number of atoms
+M = size(X,2);	%number of signals
+
+mu1 = mean(abs(Aprev(:)));			%first moment of distribution of Aprev
+mu2 = (norm(Aprev(:))^2)/numel(Aprev);%second moment of distribution of Aprev
+kappa = 2*(mu2-mu1^2)/(mu2-2*mu2^2); %parameter kappa of the MOL distribution
+beta = (kappa-1)*mu1;				%parameter beta of the MOL distribution
+
+E = X-D*Aprev;						%error term
+sigmasq = mean(var(E));				%error variance
+tau = 2*sigmasq*(kappa+1);			%sparsity factor
+
+%solve a succession of subproblems to approximate the non-convex cost
+%function
+nIter = 10;							%number of iterations of surrogate subproblem
+Psi   = zeros(K,M);					%initialise solution of subproblem
+for iIter=1:nIter
+	Reg = 1./(abs(Psi) + beta);
+	Psi = solvel1(X,D,tau,Reg);
+end
+A = Psi;
+end
+
+function Psi = solvel1(X,D,tau,A)
+	[K M] = size(A);
+	Psi = zeros(K,M);
+	for m=1:M
+		cvx_begin quiet
+			variable v(K)
+			minimise (norm(X(:,m)-D*v) + tau*norm(A(:,m).*v,1));
+		cvx_end
+		Psi(:,m) = v;
+	end
+end
+
+function D = dictionaryupdate(X,A,Dprev,zeta,eta)
+	D = (X*A' + 2*(zeta + eta)*Dprev)/(A*A' + 2*zeta*(Dprev'*Dprev) + 2*eta*diag(diag(Dprev'*Dprev)));
+end
+
+
+
+function Y = colnorms_squared(X)
+% compute in blocks to conserve memory
+Y = zeros(1,size(X,2));
+blocksize = 2000;
+for i = 1:blocksize:size(X,2)
+	blockids = i : min(i+blocksize-1,size(X,2));
+	Y(blockids) = sum(X(:,blockids).^2);
+end
+end
+
+function testdl_ramirez
+	clc
+	N = 10;		%ambient dimension
+	K = 20;		%number of atoms
+	M = 30;		%number of observed signals
+	X = randn(N,M);		%observed signals
+	D = normcol(randn(N,K)); %initial dictionary
+	Problem.b = X;		%sparse representation problem
+	Problem.b1 = X;
+	DL = SMALL_init_DL('dl_ramirez');
+	DL.param.initdict = D;
+	DL.param = struct('initdict',D,...
+					  'zeta',0.5,...
+					  'eta',0.5);
+	DL = SMALL_learn(Problem,DL);
+end

--- a/DL/two-step DL/dico_decorr.m	Wed Apr 11 16:28:46 2012 +0100
+++ b/DL/two-step DL/dico_decorr.m	Thu Apr 12 13:52:28 2012 +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	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,208 @@
+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)
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1))); %target coherence (based on coherence lower bound)
+minCoherence = [0.1:0.1:1];
+
+% 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
+name = dicUpdate{1}; %use ksvd update
+SMALL.DL(1:36) = 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
+%save('SMALL','SMALL');
+
+% learn with random initialisation and mailhe decorrelation
+for iMu=1:10
+	SMALL.DL(1+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu)); %parameters for the dictionary learning
+	SMALL.DL(1+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(1+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
+
+% 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
+for iMu=1:10
+	SMALL.DL(13+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu),'initdict',dctDict); %parameters for the dictionary learning
+	SMALL.DL(13+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(13+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
+
+% learn with dct initialisation and tropp decorrelation
+SMALL.DL(24).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(24) = SMALL_learn(SMALL.Problem,SMALL.DL(24)); %learn dictionary
+
+
+% Learn with gabor initialisation and no decorrelation
+SMALL.DL(25).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(25) = SMALL_learn(SMALL.Problem,SMALL.DL(25)); %learn dictionary
+
+% learn with gabor initialisation and mailhe decorrelation
+for iMu=1:10
+	SMALL.DL(25+iMu).param = struct('data',SMALL.Problem.b,'Tdata',nActiveAtoms,...
+		'dictsize',dictSize,'iternum',iterNum,'memusage','high','solver',solver,...
+		'decFcn','mailhe','coherence',minCoherence(iMu),'initdict',gaborDict); %parameters for the dictionary learning
+	SMALL.DL(25+iMu) = SMALL_learn(SMALL.Problem,SMALL.DL(25+iMu)); %learn dictionary
+	%save('SMALL','SMALL');
+end
+
+% learn with gabor initialisation and tropp decorrelation
+SMALL.DL(36).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(36) = SMALL_learn(SMALL.Problem,SMALL.DL(36)); %learn dictionary
+
+%% Evaluate coherence and snr of representation for the various methods
+sigNoiseRatio = zeros(36,1);
+mu = zeros(36,1);
+for i=1:36
+	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)));
+
+figure,
+%subplot(3,1,1)
+hold on, grid on
+title('Data Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(1),sigNoiseRatio(1),'ks');
+plot(mu(12),sigNoiseRatio(12),'kd');
+plot(mu(2:11),sigNoiseRatio(2:11),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
+
+figure
+%subplot(3,1,2)
+hold on, grid on
+title('DCT Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(13),sigNoiseRatio(13),'ks');
+plot(mu(24),sigNoiseRatio(24),'kd');
+plot(mu(14:23),sigNoiseRatio(14:23),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
+
+figure
+%subplot(3,1,3)
+hold on, grid on
+title('Gabor Initialisation')
+plot([1 1],[0 25],'k-');
+plot(mu(25),sigNoiseRatio(25),'ks');
+plot(mu(36),sigNoiseRatio(36),'kd');
+plot(mu(26:35),sigNoiseRatio(26:35),'k*-');
+plot([minMu minMu],[0 25],'k--')
+set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+legend({'\mu_{max}','K-SVD','Grassmannian','INK-SVD','\mu_{min}'});
+xlabel('\mu');
+ylabel('SNR (dB)');
+
+% minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));
+% maxSNR = max(sigNoiseRatio);
+%
+% figure, subplot(2,2,1)
+% snrMat = buffer(sigNoiseRatio(1:9),3);
+% bar(snrMat');
+% title('Signal to noise ratio')
+% xlabel('Initial dictionary')
+% ylabel('SNR (dB)')
+% set(gca,'XTickLabel',{'data','dct','gabor'});
+% 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')
+% 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')

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/SMALL_test_coherence2.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,128 @@
+%
+%
+%
+
+clc, clear, close all
+
+%% Parameteres
+nTrials   = 10;										%number of trials of the experiment
+
+% Dictionary learning parameters
+toolbox   = 'TwoStepDL';							%dictionary learning toolbox
+dicUpdate = 'ksvd';									%dictionary learning updates
+dicDecorr = {'iterproj','ink-svd','shrinkgram'};	%dictionary decorrelation methods
+minCoherence = linspace(0.1,1,10);					%coherence levels
+
+iterNum   = 20;				%number of iterations
+epsilon   = 1e-6;			%tolerance level
+dictSize  = 512;			%number of atoms in the dictionary
+percActiveAtoms = 5;		%percentage of active atoms
+
+% Test signal parameters
+signal    = audio('music03_16kHz.wav'); %audio signal
+blockSize = 256;						%size of audio frames
+overlap   = 0.5;						%overlap between consecutive frames
+
+% Dependent parameters
+nActiveAtoms = fix(blockSize/100*percActiveAtoms); %number of active atoms
+
+% Initial dictionaries
+gaborParam = struct('N',blockSize,'redundancyFactor',2,'wd',@rectwin);
+gaborDict = Gabor_Dictionary(gaborParam);
+initDicts = {[],gaborDict};
+
+%% Generate audio approximation problem
+signal			 = buffer(signal,blockSize,blockSize*overlap,@rectwin);	%buffer frames of audio into columns of the matrix S
+SMALL.Problem.b  = signal.S;
+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
+nDecorrAlgs = length(dicDecorr);		%number of decorrelation algorithms
+nCorrLevels = length(minCoherence);		%number of coherence levels
+nInitDicts  = length(initDicts);		%number of initial dictionaries
+
+SMALL.DL(nTrials,nInitDicts,nCorrLevels,nDecorrAlgs) = SMALL_init_DL(toolbox); %create dictionary learning structures
+for iTrial=1:nTrials
+	for iInitDicts=1:nInitDicts
+		for iCorrLevels=1:nCorrLevels
+			for iDecorrAlgs=1:nDecorrAlgs
+				SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).toolbox = toolbox;
+				SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).name = dicUpdate;
+				SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).profile = true;
+				SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).param = ...
+					struct( 'data',SMALL.Problem.b,...
+					'Tdata',nActiveAtoms,...
+					'dictsize',dictSize,...
+					'iternum',iterNum,...
+					'memusage','high',...
+					'solver',solver,...
+					'decFcn',dicDecorr{iDecorrAlgs},...
+					'coherence',minCoherence(iCorrLevels),...
+					'initdict',initDicts(iInitDicts));
+				
+				SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs) = ...
+					SMALL_learn(SMALL.Problem,SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs));
+				save('SMALL_DL','SMALL');
+			end
+		end
+	end
+end
+
+%% Evaluate coherence and snr of representation for the various methods
+sr = zeros(size(SMALL.DL));				%signal to noise ratio
+mu1 = zeros(nTrials,nInitDicts,nCorrLevels,nDecorrAlgs,blockSize);	%cumulative coherence
+mu2 = zeros(nTrials,nInitDicts,nCorrLevels,nDecorrAlgs,blockSize);	%cumulative coherence
+dic(size(SMALL.DL)) = dictionary;		%initialise dictionary objects
+for iTrial=1:nTrials
+	for iInitDicts=1:nInitDicts
+		for iCorrLevels=1:nCorrLevels
+			for iDecorrAlgs=1:nDecorrAlgs
+				%Sparse representation
+				SMALL.Problem.A = SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).D;
+				tempSolver = SMALL_solve(SMALL.Problem,solver);
+				%calculate snr
+				sr(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs) = ...
+					snr(SMALL.Problem.b,SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).D*tempSolver.solution);
+				%calculate mu
+				dic(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs) = ...
+					dictionary(SMALL.DL(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).D);
+				mu1(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs,:) = ...
+					dic(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).cumcoherence;
+				mu2(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs,:) = ...
+					dic(iTrial,iInitDicts,iCorrLevels,iDecorrAlgs).cumcoherence(2);
+			end
+		end
+	end
+end
+
+%% Plot results
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));	%lowe bound on coherence
+initDictsNames = {'Data','Gabor'};
+dicDecorrNames = {'IPR','INK-SVD'};
+lineStyles     = {'k.-','r*-','b+-'};
+for iInitDict=1:nInitDicts
+	figure, hold on, grid on
+	title([initDictsNames{iInitDict} ' Initialisation']);
+	plot([1 1],[0 25],'k-');
+	for iDecorrAlgs=1:nDecorrAlgs-1
+		coherenceLevels = squeeze(mean(mu(:,iInitDict,:,iDecorrAlgs,1),1));
+		meanSNRs		= squeeze(mean(sr(:,iInitDict,:,iDecorrAlgs),1));
+		stdSNRs			= squeeze(std(sr(:,iInitDict,:,iDecorrAlgs),0,1));
+		errorbar(coherenceLevels,meanSNRs,stdSNRs,lineStyles{iDecorrAlgs});
+	end
+	plot([minMu minMu],[0 25],'k--')
+	
+	set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+	legend([{'\mu_{max}'},dicDecorrNames,{'\mu_{min}'}]);
+	xlabel('\mu');
+	ylabel('SNR (dB)');
+end
+
+%% 
+mu2 = squeeze(mean(mu2,1));
+mu = squeeze(mean(mu,1));

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/SMALL_test_mocod.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,167 @@
+%% SMALL_test_mocod
+% Script that tests the twostep dictionary learning algorithm with MOCOD
+% dictionary update.
+%
+% REFERENCES
+% D. Barchiesi and M. D. Plumbely, Learning incoherenct dictionaries for 
+% sparse approximation using iterative projections and rotations.
+%% Clear and close
+clc, clear, close all
+
+%% Parameteres
+nTrials   = 10;					%number of trials of the experiment
+
+% Dictionary learning parameters
+toolbox   = 'TwoStepDL';		%dictionary learning toolbox
+dicUpdate = 'mocod';			%dictionary learning updates
+zeta	  = logspace(-2,2,10);	%range of values for the incoherence term				
+eta  	  = logspace(-2,2,10);	%range of values for the unit norm term
+iterNum   = 20;					%number of iterations
+epsilon   = 1e-6;				%tolerance level
+dictSize  = 512;				%number of atoms in the dictionary
+percActiveAtoms = 5;			%percentage of active atoms
+
+% Test signal parameters
+signal    = audio('music03_16kHz.wav'); %audio signal
+blockSize = 256;						%size of audio frames
+overlap   = 0.5;						%overlap between consecutive frames
+
+% Dependent parameters
+nActiveAtoms = fix(blockSize/100*percActiveAtoms); %number of active atoms
+
+% Initial dictionaries
+gaborParam = struct('N',blockSize,'redundancyFactor',2,'wd',@rectwin);
+gaborDict = Gabor_Dictionary(gaborParam);
+initDicts = {[],gaborDict};				%cell containing initial dictionaries
+
+%% Generate audio approximation problem
+%buffer frames of audio into columns of the matrix S
+signal			 = buffer(signal,blockSize,blockSize*overlap,@rectwin);	
+SMALL.Problem.b  = signal.S;
+
+%copy signals from training set b to test set b1 (needed for later functions)
+SMALL.Problem.b1 = SMALL.Problem.b;
+
+%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
+nInitDicts  = length(initDicts);%number of initial dictionaries
+nZetas = length(zeta);			%number of incoherence penalty parameters
+nEtas  = length(eta);			%number of unit norm penalty parameters
+
+%create dictionary learning structures
+SMALL.DL(nTrials,nInitDicts,nZetas,nEtas) = SMALL_init_DL(toolbox); 
+for iTrial=1:nTrials
+	for iInitDicts=1:nInitDicts
+		for iZetas=1:nZetas
+			for iEtas=1:nEtas
+				SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).toolbox = toolbox;
+				SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).name = dicUpdate;
+				SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).profile = true;
+				SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).param = ...
+					struct('data',SMALL.Problem.b,...	%observed data
+					'Tdata',nActiveAtoms,...			%active atoms
+					'dictsize',dictSize,...				%number of atoms
+					'iternum',iterNum,...				%number of iterations
+					'memusage','high',...				%memory usage
+					'solver',solver,...					%sparse approx solver
+					'initdict',initDicts(iInitDicts),...%initial dictionary
+					'zeta',zeta(iZetas),...				%incoherence penalty factor
+					'eta',eta(iEtas));					%unit norm penalty factor
+				%learn dictionary
+				SMALL.DL(iTrial,iInitDicts,iZetas,iEtas) = ...
+					SMALL_learn(SMALL.Problem,SMALL.DL(iTrial,iInitDicts,iZetas,iEtas));
+			end
+		end
+	end
+end
+
+%% Evaluate coherence and snr of representation
+sr = zeros(size(SMALL.DL));					%signal to noise ratio
+mu = zeros(iTrial,iInitDicts,iZetas,iEtas);	%coherence
+dic(size(SMALL.DL)) = dictionary;			%initialise dictionary objects
+for iTrial=1:nTrials
+	for iInitDicts=1:nInitDicts
+		for iZetas=1:nZetas
+			for iEtas=1:nEtas
+				%Sparse representation
+				SMALL.Problem.A = SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).D;
+				tempSolver = SMALL_solve(SMALL.Problem,solver);
+				%calculate snr
+				sr(iTrial,iInitDicts,iZetas,iEtas) = ...
+					snr(SMALL.Problem.b,...
+					SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).D*tempSolver.solution);
+				%calculate mu
+				dic(iTrial,iInitDicts,iZetas,iEtas) = ...
+					dictionary(SMALL.DL(iTrial,iInitDicts,iZetas,iEtas).D);
+				mu(iTrial,iInitDicts,iZetas,iEtas) = ...
+					dic(iTrial,iInitDicts,iZetas,iEtas).coherence;
+			end
+		end
+	end
+end
+
+%% Plot results
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));%lower bound on coherence
+initDictsNames = {'Data','Gabor'};							%names of initial dictionaries
+lineStyles     = {'k.-','r*-','b+-'};						
+for iInitDict=1:nInitDicts
+	figure, hold on, grid on
+	%print initial dictionary as figure title
+	DisplayFigureTitle([initDictsNames{iInitDict} ' Initialisation']);
+% 	set(gcf,'Units','Normalized');
+% 	txh = annotation(gcf,'textbox',[0.4,0.95,0.2,0.05]);
+% 	set(txh,'String',[initDictsNames{iInitDict} ' Initialisation'],...
+% 		'LineStyle','none','HorizontalAlignment','center');
+	%calculate mean coherence levels and SNRs over trials
+	coherenceLevels = squeeze(mean(mu(:,iInitDict,:,:),1));
+	meanSNRs		= squeeze(mean(sr(:,iInitDict,:,:),1));
+%	stdSNRs		= squeeze(std(sr(:,iInitDict,iZetas,iEtas),0,1));
+	%plot coherence levels
+	subplot(2,2,1)
+	surf(eta,zeta,coherenceLevels);
+	set(gca,'Xscale','log','Yscale','log','ZLim',[0 1.4]);
+	view(gca,130,20)
+	xlabel('\eta');
+	ylabel('\zeta');
+	zlabel('\mu');
+	title('Coherence')
+	
+	%plot SNRs
+	subplot(2,2,2)
+	surf(eta,zeta,meanSNRs);
+	set(gca,'Xscale','log','Yscale','log','ZLim',[0 25]);
+	view(gca,130,20)
+	xlabel('\eta');
+	ylabel('\zeta');
+	zlabel('SNR (dB)');
+	title('Reconstruction Error')
+	
+	%plot mu/SNR scatter
+	subplot(2,2,[3 4])
+	mus = mu(:,iInitDict,:,:);
+	mus = mus(:);
+	SNRs = sr(:,iInitDict,:,:);
+	SNRs = SNRs(:);
+	[un idx] = sort(mus);
+	plot([1 1],[0 25],'k')
+	hold on, grid on
+	scatter(mus(idx),SNRs(idx),'k+');
+	plot([minMu minMu],[0 25],'k--')
+	set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+	xlabel('\mu');
+	ylabel('SNR (dB)');
+	legend([{'\mu_{max}'},'MOCOD',{'\mu_{min}'}]);
+	title('Coherence-Reconstruction Error Tradeoff')
+	
+% 	plot([minMu minMu],[0 25],'k--')
+% 	
+% 	set(gca,'YLim',[0 25],'XLim',[0 1.4]);
+% 	legend([{'\mu_{max}'},dicDecorrNames,{'\mu_{min}'}]);
+% 	xlabel('\mu');
+% 	ylabel('SNR (dB)');
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/DisplayFigureTitle.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,5 @@
+function DisplayFigureTitle(title)
+
+set(gcf,'Units','Normalized');
+txh = annotation(gcf,'textbox',[0.4,0.95,0.2,0.05]);
+set(txh,'String',title,'LineStyle','none','HorizontalAlignment','center');

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/SMALL_ImgDeNoiseResult.m.orig	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,50 @@
+function SMALL_ImgDeNoiseResult(SMALL)   
+%%  Represents the results of Dictionary Learning for Image denoising
+%
+%   Function gets as input SMALL structure and plots  Image Denoise
+%   results: Original Image, Noisy Image and for learned dictionaries and 
+%   denoised images 
+%
+
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2010 Ivan Damnjanovic.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   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.
+%%
+
+
+figure('Name', sprintf('Image %s (training set size- %d, sigma - %d)',SMALL.Problem.name, SMALL.Problem.n, SMALL.Problem.sigma));
+
+m=size(SMALL.solver,2)+1;
+maxval=SMALL.Problem.maxval;
+im=SMALL.Problem.Original;
+imnoise=SMALL.Problem.Noisy;
+
+subplot(2, m, 1); imagesc(im/maxval);colormap(gray);axis off; axis image;        % Set aspect ratio to obtain square pixels
+title('Original image');
+
+subplot(2,m,m+1); imagesc(imnoise/maxval);axis off; axis image; 
+title(sprintf('Noisy image, PSNR = %.2fdB', SMALL.Problem.noisy_psnr ));
+
+for i=2:m
+    
+    subplot(2, m, i); imagesc(SMALL.solver(i-1).reconstructed.Image/maxval);axis off; axis image; 
+    title(sprintf('%s Denoised image, PSNR: %.2f dB in %.2f s',...
+        SMALL.DL(i-1).name, SMALL.solver(i-1).reconstructed.psnr, SMALL.solver(i-1).time ),'Interpreter','none');
+    if strcmpi(SMALL.DL(i-1).name,'ksvds')
+        D = kron(SMALL.Problem.basedict{2},SMALL.Problem.basedict{1})*SMALL.DL(i-1).D;
+    else
+        D = SMALL.DL(i-1).D;
+    end
+    dictimg = SMALL_showdict(D,SMALL.Problem.blocksize,...
+        round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast');
+    
+    subplot(2,m,m+i);imagesc(dictimg);axis off; axis image; 
+    title(sprintf('%s dictionary in %.2f s',...
+        SMALL.DL(i-1).name, SMALL.DL(i-1).time),'Interpreter','none');
+    
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/audio.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,77 @@
+%% AUDIO OBJECT CLASS
+% Class designed to analyse and process audio signals
+%%
+classdef audio
+	properties (SetAccess = protected)
+		s				%vector containing the audio signal
+		fs				%sampling frequency
+		nBits			%number of bits per sample
+		name			%string containing the name of the audio file
+		format			%string containing the format of the audio file
+		bufferOperator	%struct containing the parameters of the buffer operator
+		S				%matrix containing frames of audio
+	end
+	
+	methods
+		%% Constructor
+		function obj = audio(varargin)
+			%%% obj = audio(varargin)
+			% Audio object constructor.
+			% INPUT: either a path to an audio file, or the following
+			% arguments.
+			% - s: vector containing the audio samples
+			% - fs: sampling frequency
+			% - nBits: number of bits per sample
+			% - name: name of the audio object
+			% - format: format of the audio object
+			%
+			% if no arguments are specified, prompt for the choice of an
+			% audio file
+			if ~nargin
+				[fileName,pathname] = uigetfile({'*.wav; *.aiff;'},'Select an audio file');
+				varargin{1} = strcat(pathname,filesep,fileName);
+			end
+			% if a file is specified, read it from disk
+			if ischar(varargin{1})
+				[~, obj.name obj.format] = fileparts(varargin{1});
+				switch obj.format
+					case '.wav'
+						[obj.s obj.fs obj.nBits] = wavread(varargin{1});
+					otherwise
+						error('Unsupported audio format')
+				end
+			% if properties are specified, set them to input values
+			else
+				obj.s = varargin{1};
+				if nargin>1, obj.fs = varargin{2}; else obj.fs = []; end
+				if nargin>2, obj.nBits = varargin{3}; else obj.nBits = []; end
+				if nargin>3, obj.name = varargin{4}; else obj.name = []; end
+				if nargin>4, obj.format = varargin{5}; else obj.format = []; end
+			end
+			obj.S = [];
+			obj.bufferOperator = [];
+		end
+		
+		%% Playback functions
+		function player = play(obj, player)
+			if ~exist('player','var') || isempty(player)
+				player = audioplayer(obj.s,obj.fs);
+			end
+			play(player);
+		end
+		
+		function player = stop(obj, player)
+			if ~exist('player','var') || isempty(player)
+				player = audioplayer(obj.s,obj.fs);
+			end
+			stop(player)
+		end
+		
+		function player = pause(obj, player)
+			if ~exist('player','var') || isempty(player)
+				player = audioplayer(obj.s,obj.fs);
+			end
+			pause(player)
+		end
+	end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/buffer.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,33 @@
+%% Buffer function
+% Buffers the samples of the audio object into the columns of the matrix S
+% based on the input parameters
+%%
+function obj = buffer(obj,wLength,overlap,window,method)
+
+%% Check inputs & Defaults
+error(nargchk(2, 5, nargin, 'struct'));
+if rem(length(obj.s),wLength)
+%	error('The wLength must be an integer divisor of the signal length!'); 
+end
+if ~exist('overlap','var') || isempty(overlap), overlap = 0; end
+if ~exist('method','var') || isempty(method), method = 'standard'; end
+
+%% Buffer audio
+switch lower(method)
+	case 'standard'
+		if ~exist('window','var') || isempty(window), window = @rectwin; end
+		validWindows = {'hanning','hamming','triang','rectwin'};
+		if ~sum(strcmpi(validWindows,func2str(window)));
+			error('The window chosen is not valid because it cannot be inverted!');
+		end
+		obj.S = diag(window(wLength))*buffer(obj.s,wLength,overlap,'nodelay');
+% 	case 'lot'
+% 		if ~exist('window','var') || isempty(window), window = 'sin2'; end
+% 		s_lot = lot(obj.s,wLength,'id',overlap,window);
+% 		obj.S = buffer(s_lot,wLength);
+	otherwise
+		error('Please specify a valid buffer method');
+end
+
+obj.bufferOperator = struct('wLength',wLength,'overlap',...
+	overlap,'window',window,'method',method);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/plot.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,92 @@
+function plot(obj)
+
+figure,
+playbackPanelH = uipanel(gcf,'Units','Normalized','Position',[.3 0 .4 .1]);
+
+buttWidth = 1/6;
+centers = linspace(0,1,7)-buttWidth/2;
+rewButtonH = uicontrol(playbackPanelH,'Style','pushbutton','String','<<','Units',...
+	'Normalized','Position',[centers(2) 0.2 buttWidth 0.6]);
+ffButtonH = uicontrol(playbackPanelH,'Style','togglebutton','String','>>','Units',...
+	'Normalized','Position',[centers(3) 0.2 buttWidth 0.6]);
+stopButtonH = uicontrol(playbackPanelH,'Style','pushbutton','String','stop','Units',...
+	'Normalized','Position',[centers(4) 0.2 buttWidth 0.6]);
+playButtonH = uicontrol(playbackPanelH,'Style','togglebutton','String','play','Units',...
+	'Normalized','Position',[centers(5) 0.2 buttWidth 0.6]);
+pauseButtonH = uicontrol(playbackPanelH,'Style','togglebutton','String','||','Units',...
+	'Normalized','Position',[centers(6) 0.2 buttWidth 0.6]);
+
+%% Plot the time domain signal
+s = obj.s;
+fs = obj.fs;
+plot((1:length(s))/fs,s);
+title('Audio signal')
+xlabel('time (s)');
+axis tight
+
+player = audioplayer(s,fs);
+set(player,'TimerPeriod',0.1);
+set(player,'StartFcn',@plotTransportBar);
+set(player,'TimerFcn',@updateTransportBar);
+set(player,'StopFcn',@deleteTransportBar);
+
+%% Add playback controls
+set(playButtonH,'Callback',@play_callback);
+set(stopButtonH,'Callback',@stop_callback);
+set(pauseButtonH,'Callback',@pause_callback);
+set(rewButtonH,'Callback',@rew_callback);
+set(ffButtonH,'Callback',@ff_callback);
+
+	function play_callback(~,~)
+		set(player,'SampleRate',fs);
+		play(player,player.CurrentSample);
+		set(pauseButtonH,'Value',0);
+		set(ffButtonH,'Value',0);
+	end
+
+	function pause_callback(~,~)
+		pause(player);
+		set(playButtonH,'Value',0);
+		set(ffButtonH,'Value',0);
+	end
+
+	function stop_callback(~,~)
+		stop(player);
+		set(playButtonH,'Value',0);
+		set(pauseButtonH,'Value',0);
+		set(ffButtonH,'Value',0);
+	end
+
+	function ff_callback(~,~)
+		set(player,'SampleRate',1.5*fs);
+		set(pauseButtonH,'Value',0);
+		set(playButtonH,'Value',0);
+	end
+
+	function rew_callback(~,~)
+		stop(player);
+		play(player);
+		set(pauseButtonH,'Value',0);
+		set(playButtonH,'Value',1);
+	end
+
+%% Transport Bar functions
+	function plotTransportBar(~,~)
+		global tbH
+		xLim = get(gca,'Xlim');
+		yLim = get(gca,'YLim');
+		tbH = line([xLim(1) xLim(1)],yLim,'Color','k');
+	end
+
+	function updateTransportBar(hObject,~)
+		global tbH
+		currentSample = hObject.CurrentSample;
+		pos = currentSample/fs;
+		set(tbH,'XData',[pos pos]);
+	end
+
+	function deleteTransportBar(~,~)
+		global tbH
+		delete(tbH);
+	end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/unbuffer.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,27 @@
+function s = unbuffer(obj)
+
+%% Check inputs and Defaults
+if ~isprop(obj,'bufferOperator') || isempty(obj.bufferOperator)
+	error('You must buffer a signal before unbuffer it, come on!');
+end
+
+switch lower(obj.bufferOperator.method)
+	%Unbuffer using overlap-add method
+	case 'standard'
+		w = obj.bufferOperator.window(obj.bufferOperator.wLength);
+		S = diag(1./w)*obj.S;
+		%Non overlapping case
+		if obj.bufferOperator.overlap == 0
+			s = S(:);
+		%Overlapping case
+		else
+			Stemp = S(1:obj.bufferOperator.wLength-obj.bufferOperator.overlap,1:end);
+			s = [Stemp(:); S(obj.bufferOperator.wLength-obj.bufferOperator.overlap+1:end,end)];
+		end
+	%Unbuffer using inverse lot with identity local transform
+	case 'lot'
+		s = ilot(obj.S(:),obj.bufferOperator.wLength,'id',...
+			obj.bufferOperator.overlap,obj.bufferOperator.window);
+	otherwise
+		error('Please specify a valid buffer method');
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/cumcoherence.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,30 @@
+function mu = cumcoherence(obj,p)
+% Calculates the p-cumulative coherence of the dictionary, defined as
+% \mu_p(k) = max_|I|=k{max_j\notin I{(\sum_i \in I}|<\phi_i,\phi_j>|^p)^1/p}
+%
+% INPUT
+% obj: dictionary object
+% p: power (default 1)
+%
+% OUTPUT
+% mu: p-cumulative coherence
+if ~exist('p','var') || isempty(p), p = 1; end
+
+obj = normalize(obj);
+[M N] = size(obj.phi);
+mu = zeros(M,1);
+for m=1:M
+	c = zeros(N);
+	for i=1:N
+		c(:,i) = abs(obj.phi'*obj.phi(:,i)).^p;
+		c(i,i) = 0;
+	end
+	c = sort(c,'descend');
+	c = c(1:m,:);
+	if m==1
+		mu(m) = max(c.^(1/p));
+	else
+		mu(m) = max(sum(c).^(1/p));
+	end
+end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/dictionary.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,111 @@
+classdef dictionary
+    %% Dictionary for sparse representation
+    properties
+        phi         %Matrix containing the dictionary
+        len         %Length of basis functions
+        nAtoms      %Number of basis function
+        name        %String containing the matrix ensemble from which the dictionary is drawn
+    end
+    properties (Dependent = true)
+        redundancy      %Redundancy of the dictionary: nAtoms/len
+        coherence       %Maximum inner product of different basis
+        isNormalised    %True if the atoms have unit norm
+        rank            %rank of the dictionary
+    end
+    
+    methods
+        %% Constructor
+        function obj = dictionary(phi,len,nAtoms)
+            % obj = dictionary(phi,len,nAtoms)
+            % INPUTS:
+            % - phi: either a string specifying a matrix ensamble or a
+            % matrix defining an explicit dictionary
+            % - len: length of the atoms (only for implicit dictionaries)
+            % - nAtoms: number of atoms (only for implicit dictionaries)
+            if nargin
+                if ~ischar(phi)
+                    [obj.len obj.nAtoms] = size(phi);
+                    obj.phi              = phi;
+                    obj.name             = 'explicit';
+                else
+                    switch lower(phi)
+                        case 'dct'
+                            obj.phi = dctmatrix(len,nAtoms);
+                        case 'grassmannian'
+                            obj.phi = grassmanian(len,nAtoms);
+                        otherwise
+                            obj.phi = MatrixEnsemble(len,nAtoms,phi);
+                    end
+                    obj.len    = len;
+                    obj.nAtoms = nAtoms;
+                    obj.name   = lower(phi);
+                end
+            end
+		end
+        %% Dependent properties
+        function redundancy = get.redundancy(obj)
+            redundancy = obj.nAtoms/obj.len;
+        end
+        function coherence = get.coherence(obj)
+            obj.phi = normcol(obj.phi);
+            G = obj.phi'*obj.phi;
+            G = G - eye(size(G));
+            coherence = max(abs(G(:)));
+        end
+        function isNormalised = get.isNormalised(obj)
+            isNormalised = norm(sum(conj(obj.phi).*obj.phi) - ...
+                ones(1,obj.nAtoms))<1e-9;
+        end
+        function r = get.rank(obj)
+            r = rank(obj.phi);
+        end
+        %% Operations
+        function obj = normalize(obj)
+            obj.phi = normcol(obj.phi);
+		end
+        %% Visualization
+        function image(obj)
+            %Image of the dictionary
+            if isreal(obj.phi)
+                imagesc(obj.phi);
+                title('Dictionary');
+                xlabel('Atom number');
+            else
+                subplot(2,1,1)
+                imagesc(real(obj.phi));
+                title('Real');
+                xlabel('Atom number');
+                subplot(2,1,2)
+                imagesc(imag(obj.phi));
+                title('Imaginary');
+                xlabel('Atom number');
+            end
+        end
+        function imagegram(obj)
+            G = obj.phi'*obj.phi;
+            imagesc(G);
+            title('Gram Matrix')
+        end
+        function plot(obj,n)
+            %Plot of the n-th basis
+            if isreal(obj.phi)
+                plot(obj.phi(:,n));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2))]);
+            else
+                subplot(2,1,1)
+                plot(real(obj.phi(:,n)));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2)) ' - Real']);
+                subplot(2,1,2)
+                plot(imag(obj.phi(:,n)));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2)) ' - Imaginary']);
+            end
+		end
+        function movie(obj)
+            %Movie of the basis
+            for i=1:size(obj.phi,2)
+                obj.plot(i);
+                pause(1/25);
+            end
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/mtimes.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,11 @@
+function C = mtimes(A,B)
+isAdic = strcmp(class(A),'dictionary');
+isBdic = strcmp(class(B),'dictionary');
+if isAdic && ~isBdic
+    C = A.phi*B;
+elseif isBdic && ~isAdic
+    C = A*B.phi;
+elseif isAdic && isBdic
+    C = A.phi*B.phi;
+end
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/plotcumcoherencebounds.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,65 @@
+function plotcumcoherencebounds(obj)
+
+[d N] = size(obj.phi);
+
+mu1 = cumcoherence(obj,1);
+mu2 = cumcoherence(obj,2);
+mu1_min = (1:d)*sqrt((N-d)/(d*(N-1)));
+
+figure, subplot(2,1,1)
+hold on, grid on
+plot(1:d,mu1,'k-')
+plot(1:d,mu2,'k--')
+plot(1:d,mu1_min,'b-.')
+set(gca,'XScale','log','YScale','log');
+axis tight
+ylabel('p-cumulative coherence')
+xlabel('k')
+legend('\mu_1(k)','\mu_2(k)','\mu_{1_{min}}(k)')
+
+% temp = conv(mu1,[1 1]);
+% kMax = find(temp<1,1,'last');
+% title(['Worst case bound: k_{max} = ' num2str(kMax)])
+% 
+% % subplot(2,1,2)
+% p_BP_fails = exp(-1./(8*mu1(1)^2*(1:d)));
+% p_TH_fails = exp(-1./(128*mu2.^2));
+% plot(1:d,p_BP_fails,'k-+');
+% plot(1:d,p_TH_fails,'k-o');
+% legend('P(BP fails)','P(TH fails)');
+
+% % Plot cumulative coherence with lower and upper bounds
+% mumin = (1:d)*sqrt((N-d)/(d*(N-1)));
+% mumax = (1:d);
+% figure,
+% subplot(1,6,1:2)
+% hold on, grid on
+% plot(1,mumax(1),'k-s');
+% plot(1,mu(1),'ko');
+% plot(1,mumin(1),'k-d')
+% 
+% 
+% subplot(1,6,3:6)
+% hold on, grid on
+% plot(2:d,mumax(2:end),'k-s');
+% plot(2:d,mu(2:end),'k-o');
+% plot(2:d,mumin(2:end),'k-d');
+% set(gca,'XScale','log','YScale','log');
+% axis tight
+% xlabel('k');
+% ylabel('\mu(k)');
+% ylim([mumin(1) 10])
+% 
+% % Plot recovery bounds
+% plot(2:d,1-mu(1:d-1),'r-o')
+% plot(2:d,1-mumin(1:d-1)','r-d')
+% plot([2 d],[1/3 1/3],'b');
+% legend('\mu_{max}(k)','\mu(k)','\mu_{min}(k)','Exact-Sparse \mu','Exact-Sparse \mu_{min}','Sparse');
+
+
+% v = conv(mu,[1 1]);
+% ind = find(v<1, 1, 'last');
+% 
+% line([ind ind], [min(mu) max(mu)],'Color','red');
+% title(['Minimum allowed sparsity (Tanner):' num2str(ind/obj.len)]);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/spark.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,35 @@
+function varargout = Spark(obj)
+% Calculates the minimum number of linearly dependent atoms in the matrix A
+% WARNING: this function computes a combinatorial function, use only if the
+% size of the problem is small (i.e. <20)
+if nargout <= 1
+    A = obj.phi;
+    k = size(A,2);
+    if k>20
+        warning('This function computes a combinatorial function, use only if thesize of the problem is small (i.e. <20).');
+        choice = input('The calculation of spark will take a long time... do you wish to continue anyway (y/n)','s');
+        if strcmpi( choice,'n')
+            return
+        end
+    end
+    sigma = 2;
+    while true
+        P = nchoosek(1:k,sigma);
+        for i=1:size(P,1)
+            r = rank(A(:,P(i,:)));
+            if r<sigma
+                varargout{1} = sigma;
+                return
+            end
+        end
+        sigma = sigma + 1;
+        if sigma==k
+            varargout{1} = inf;
+            return
+        end
+    end
+else
+    %% TODO: calculate lower and upper bounds on the spark
+    varargout{1} = 2;
+    varargout{2} = inf;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/dctmatrix.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,24 @@
+function D = dctmatrix(N,K,type)
+
+error(nargchk(1,3,nargin,'struct'));
+if ~exist('type','var') || isempty(type), type='II'; end
+if ~exist('K','var') || isempty(K), K=N; end
+
+[c r] = meshgrid(0:K-1,0:N-1);
+switch type
+    case 'I'
+        D = cos(pi*c.*r/(K-1));
+        D(1,:) = D(1,:)/sqrt(2);
+        D(end,:) = D(end,:)/sqrt(2);
+    case 'II'
+        D = cos(pi*(2*c+1).*r/(2*K));
+        D(1,:) = D(1,:)/sqrt(2);
+    case 'III'
+        D = cos(pi*(2*r+1).*c/(2*K));
+        D(:,1) = D(:,1)/sqrt(2);
+    case 'IV'
+        D = cos(pi*(2*r+1+2*c+4*c.*r)/(4*K));
+    otherwise
+        error('unsupported dct type');
+end
+D = normcol(D);
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/grassmannian.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,47 @@
+function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
+% grassmanian attempts to create an n by m matrix with minimal mutual
+% coherence using an iterative projection method.
+%
+% [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
+%
+% REFERENCE
+% M. Elad, Sparse and Redundant Representations, Springer 2010.
+
+%% Parameters and Defaults
+error(nargchk(2,7,nargin));
+
+if ~exist('verb','var')  || isempty(verb),  verb = false; end %verbose output
+if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
+if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.9; end %shrinking factor
+if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
+if ~exist('nIter','var') || isempty(nIter), nIter = 10; end %number of iterations
+
+%% Main algo
+A = normc(initA); %normalise columns
+[Uinit Sigma] = svd(A);
+G = A'*A; %gram matrix
+
+muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence (equiangular tight frame)
+res = zeros(nIter,1);
+if verb
+	fprintf(1,'Iter		mu_min		mu \n');
+end
+
+% optimise gram matrix
+for iIter = 1:nIter
+	gg  = sort(abs(G(:))); %sort inner products from less to most correlated
+	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
+	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
+	[U S V] = svd(G);	%compute new SVD of gram matrix
+	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
+	G = U*S*V';			%update gram matrix
+	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
+	if verb
+		Geye = G - eye(size(G));
+		fprintf(1,'%6i    %12.8f  %12.8f  \n',iIter,muMin,max(abs(Geye(:))));
+	end
+end
+
+[V_gram Sigma_gram] = svd(G);				%calculate svd decomposition of gramian
+Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma); %calculate singular values of dictionary
+A = Uinit*Sigma_new*V_gram';				%update dictionary

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/iterativeprojections.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,97 @@
+function [dic mus srs] = iterativeprojections(dic,mu,Y,X,params)
+% grassmanian attempts to create an n by m matrix with minimal mutual
+% coherence using an iterative projection method.
+%
+% REFERENCE
+%
+
+%% Parameters and Defaults
+if ~nargin, testiterativeprojections; return; end
+
+if ~exist('params','var') || isempty(param), params = struct; end
+if ~isfield(params,'nIter'), params.nIter = 10; end %number of iterations
+if ~isfield(params,'eps'),	 params.eps = 1e-9;	 end %tolerance level
+[n m] = size(dic);
+
+SNR = @(dic) snr(Y,dic*X);									 %SNR function
+MU  = @(dic) max(max(abs((dic'*dic)-diag(diag(dic'*dic))))); %coherence function
+
+%% Main algorithm
+dic = normc(dic);	%normalise columns
+alpha = m/n;		%ratio between number of atoms and ambient dimension
+
+mus = zeros(params.nIter,1);		%coherence at each iteration
+srs = zeros(params.nIter,1);		%signal to noise ratio at each iteration
+iIter = 1;
+while iIter<=params.nIter && MU(dic)>mu
+	fprintf(1,'Iteration number %u\n', iIter);
+	% calculate snr and coherence
+	mus(iIter) = MU(dic);
+	srs(iIter) = SNR(dic);
+	
+	% calculate gram matrix
+	G = dic'*dic;
+	
+	% project onto the structural constraint set
+	H = zeros(size(G));				%initialise matrix
+	ind1 = find(abs(G(:))<=mu);		%find elements smaller than mu
+	ind2 = find(abs(G(:))>mu);		%find elements bigger than mu
+	H(ind1) = G(ind1);				%copy elements belonging to ind1
+	H(ind2) = mu*sign(G(ind2));		%threshold elements belonging to ind2
+	H(1:m+1:end) = 1;				%set diagonal to one
+	
+	% project into spectral constraint set
+	[~ , S, V] = svd(H);
+	%G = alpha*(V(:,1:n)*V(:,1:n)');
+	G = V(:,1:n)*S(1:n,1:n)*V(:,1:n)';
+	
+	% calculate dictionary
+	[~, S V] = svd(G);
+	dic = sqrt(S(1:n,:))*V';
+	
+	% rotate dictionary
+	options = struct('nIter',100,'step',0.001);
+	[~, ~, W] = rotatematrix(Y,dic*X,'conjgradlie',options);
+	dic = W*dic;
+	
+	iIter = iIter+1;
+end
+if iIter<params.nIter
+	mus(iIter:end) = mus(iIter);
+	srs(iIter:end) = srs(iIter);
+end
+% Test function
+function testiterativeprojections
+clc
+%define parameters
+n = 256;							%ambient dimension
+m = 512;							%number of atoms
+N = 1024;							%number of signals
+mu_min = sqrt((m-n)/(n*(m-1)));		%minimum coherence
+
+%initialise data
+X = sprandn(m,N,1);					%matrix of coefficients
+phi = normc(randn(n,m));			%dictionary
+temp = randn(n);
+W = expm(0.5*(temp-temp'));			%rotation matrix
+Y = W*phi*X;						%observed signals
+
+%optimise dictionary
+[~, mus srs] = iterativeprojections(phi,0.2,Y,X);
+
+%plot results
+nIter = length(mus);
+
+figure, 
+subplot(2,1,1)
+plot(1:nIter,srs,'kd-');
+xlabel('nIter');
+ylabel('snr (dB)');
+grid on
+
+subplot(2,1,2), hold on
+plot(1:nIter,mus,'ko-');
+plot([1 nIter],[mu_min mu_min],'k')
+grid on
+legend('\mu','\mu_{min}');
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/rotatematrix.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,114 @@
+function [Dhat cost W] = rotatematrix(D,Phi,method,param)
+% 
+%
+%
+% REFERENCE
+% M.D. Plumbley, Geometrical Methods for Non-Negative ICA: Manifolds, Lie
+% Groups and Toral Subalgebra, Neurocomputing
+
+%% Parse inputs and set defaults
+if ~nargin, testrotatematrix; return, end
+
+if ~exist('param','var')  || isempty(param),  param  = struct; end
+if ~exist('method','var') || isempty(method), method = 'conjgradLie'; end
+if ~isfield(param,'nIter'), param.nIter = 100; end	%number of iterations
+if ~isfield(param,'eps'), param.eps = 1e-9; end		%tolerance level
+if ~isfield(param,'step'), param.step = 0.01; end
+
+J = @(W) 0.5*norm(D-W*Phi,'fro');		%cost function
+
+% Initialise variables
+cost = zeros(param.nIter,1);			%cost at each iteration
+W	 = eye(size(Phi,1));				%rotation matrix
+grad = ones(size(W));					%gradient
+t	 = param.step;						%step size
+Gprev = 0;								%previous gradient
+Hprev = 0;								%previous Lie search direction
+iIter = 1;								%iteration counter
+
+%% Main algorithm
+while iIter<=param.nIter && norm(grad,'fro')>eps
+	cost(iIter) = J(W);				%calculate cost
+	grad = (W*Phi-D)*Phi';			%calculate gradient
+	switch method
+		case 'unconstrained'		% gradient descent
+			eta = param.step;
+			W = W - eta*grad;		% update W by steepest descent
+		case 'tangent'				% self correcting tangent
+			eta = param.step;
+			W = W - 0.5*eta*(grad - W*grad'*W);
+			[U , ~, V] = svd(W);
+			W = U*V';
+		case 'steepestlie'			%steepest descent in Lie algebra
+			eta = param.step;
+			B = 2*skew(grad*W');	% calculate gradient in Lie algebra
+			W = expm(-eta*B)*W;		% update W by steepest descent
+		case 'linesearchlie'		% line search in Lie algebra
+			B = 2*skew(grad*W');	% calculate gradient in Lie algebra
+			H = -B;					% calculate direction as negative gradient
+			t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra
+			W = expm(t*H)*W;		% update W by line search
+		case 'conjgradlie'			% conjugate gradient in Lie algebra
+			G = 2*skew(grad*W');	% calculate gradient in Lie algebra
+			H = -G + polakRibiere(G,Gprev)*Hprev; %calculate conjugate gradient direction
+			t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra
+			W = expm(t*H)*W;		% update W by line search
+			Hprev = H;				% save search direction
+			Gprev = G;				% save gradient
+	end
+	iIter = iIter+1;				% update iteration counter
+end
+Dhat = W*Phi;						%rotate matrix
+cost(iIter:end) = cost(iIter-1);	%zero-pad cost vector
+end
+
+%% Support functions
+function gamma = polakRibiere(G,Gprev)
+%Polak-Ribiere rule for conjugate direction calculation
+gamma = G(:)'*(G(:)-Gprev(:))/(norm(Gprev(:))^2);
+if isnan(gamma) || isinf(gamma)
+	gamma = 0;
+end
+end
+
+function t = searchline(J,H,W,t)
+%Line search in one-parameter Lie subalgebra
+t = fminsearch(@(x) J(expm(x*H)*W),t);
+end
+
+function B = skew(A)
+%Skew-symmetric matrix
+B = 0.5*(A - A');
+end
+
+
+%% Test function
+function testrotatematrix
+clear, clc, close all
+n = 256;							%ambient dimension
+m = 512;							%number of atoms
+param.nIter = 300;				%number of iterations
+param.step  = 0.001;			%step size
+param.mu    = 0.01;			%regularization factor (for tangent method)
+methods = {'tangent','linesearchlie','conjgradlie'};
+
+Phi = randn(n,m);				%initial dictionary
+Qtrue = expm(skew(randn(n)));	%rotation matrix
+D = Qtrue*Phi;					%target dictionary
+
+cost = zeros(param.nIter,length(methods));
+times = zeros(param.nIter,length(methods));
+for iIter=1:length(methods)
+	tic
+	[~, cost(:,iIter)] = rotatematrix(D,Phi,methods{iIter},param);
+	times(:,iIter) = linspace(0,toc,param.nIter);
+	sprintf('Method %s completed in %f seconds \n',methods{iIter},toc)
+end
+
+figure, plot(times,cost)
+set(gca,'XScale','lin','Yscale','log')
+legend(methods)
+grid on
+xlabel('time (sec)')
+ylabel('J(W)')
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/shrinkgram.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,67 @@
+function [dic mus] = shrinkgram(dic,mu,dd1,dd2,params)
+% grassmanian attempts to create an n by m matrix with minimal mutual
+% coherence using an iterative projection method.
+%
+% [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
+%
+% REFERENCE
+% M. Elad, Sparse and Redundant Representations, Springer 2010.
+
+%% Parameters and Defaults
+if ~nargin, testshrinkgram; return; end
+
+if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.9; end %shrinking factor
+if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
+if ~exist('params','var') || isempty(params), params = struct; end
+if ~isfield(params,'nIter'), params.nIter = 100; end
+
+%% Main algo
+dic = normc(dic);	%normalise columns
+G = dic'*dic;		%gram matrix
+[n m] = size(dic);
+
+MU  = @(G) max(max(abs(G-diag(diag(G))))); %coherence function
+
+mus   = ones(params.nIter,1);
+iIter = 1;
+% optimise gram matrix
+while iIter<=params.nIter && MU(G)>mu
+	mus(iIter) = MU(G);		%calculate coherence
+	gg  = sort(abs(G(:)));	%sort inner products from less to most correlated
+	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
+	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
+	[U S V] = svd(G);		%compute new SVD of gram matrix
+	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
+	G = U*S*V';				%update gram matrix
+	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
+	iIter = iIter+1;
+end
+%if iIter<params.nIter
+%	mus(iIter:end) = mus(iIter-1);
+%end
+
+[V_gram Sigma_gram] = svd(G);				%calculate svd decomposition of gramian
+dic = sqrt(Sigma_gram(1:n,:))*V_gram';		%update dictionary
+
+function testshrinkgram
+clc
+%define parameters
+n = 256;							%ambient dimension
+m = 512;							%number of atoms
+N = 1024;							%number of signals
+mu_min = sqrt((m-n)/(n*(m-1)));		%minimum coherence
+
+%initialise data
+phi = normc(randn(n,m));			%dictionary
+
+%optimise dictionary
+[~, mus] = shrinkgram(phi,0.2);
+
+%plot results
+nIter = length(mus);
+
+figure, hold on
+plot(1:nIter,mus,'ko-');
+plot([1 nIter],[mu_min mu_min],'k')
+grid on
+legend('\mu','\mu_{min}');

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/tests/test_buffer.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,48 @@
+%% Buffer Test
+% Test script for the function audio->buffer which takes an audio object as
+% an input and buffers it into frames. The test assesses wether all the
+% possible buffer methods are invertible.
+%%
+
+N = 500;					% number of audio samples
+TOL = 200;					% tolerance (in decibels)
+verb = true;				% verbose
+
+obj = audio(randn(N,1));	% audio object
+% wLengts: one window, two windows, maximum number of windows
+temp = factor(N);
+wLengths = [N; fix(N/2); fix(N/temp(end))];
+% overlaps: zero, half window, max
+overlapNames = {'zero','half-window','maximum'};
+overlaps = {@(n)0, @(n)n/2, @(n)n-1};
+% windows: valid windows
+windows = {@hanning,@hamming,@triang,@rectwin};
+% methods: valid methods
+methods = {'standard'};
+
+
+nLen  = length(wLengths);
+nOver = length(overlaps);
+nWin  = length(windows);
+nMet  = length(methods);
+count = 1;
+for iLen=1:nLen
+	for iOver=1:nOver
+		for iWin=1:nWin
+			for iMet=1:nMet
+				if verb
+				printf('\n buffer test %d/%d - %d window length, %s overlap, %s window and %s method ... ',...
+					count,nMet*nWin*nOver*nLen,wLengths(iLen),overlapNames{iOver},func2str(windows{iWin}),methods{iMet});
+				end
+				obj = buffer(obj,wLengths(iLen),overlaps{iOver}(wLengths(iLen)),windows{iWin},methods{iMet});
+				s_rec = obj.unbuffer;
+				if snr(obj.s,s_rec) > TOL
+					if verb, printf('Passed'); end
+				else
+					error('Test failed!');
+				end
+				count = count+1;
+			end
+		end
+	end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/snr.m	Thu Apr 12 13:52:28 2012 +0100
@@ -0,0 +1,6 @@
+function x = snr(s,r)
+% SNR calculates the signal-to-noise ratio between the signal s and its
+% representation r defined as:
+% SNR = 20*log10(||s||_2/||s-r||_2)
+err = s-r;
+x = mag2db(norm(s(:))/norm(err(:)));
\ No newline at end of file