changeset 225:f0441ef3ae8f danieleb

Updating this branch in order to reflect the default branch.

--- a/DL/dl_ramirez.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,189 +0,0 @@
-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/examples/SMALL_test_coherence.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,208 +0,0 @@
-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')

--- a/examples/SMALL_test_coherence2.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,128 +0,0 @@
-%
-%
-%
-
-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));

--- a/examples/SMALL_test_mocod.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,167 +0,0 @@
-%% 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

--- a/util/classes/@audio/audio.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-%% 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

--- a/util/classes/@audio/buffer.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-%% 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);

--- a/util/classes/@audio/plot.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,92 +0,0 @@
-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

--- a/util/classes/@audio/unbuffer.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,27 +0,0 @@
-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

--- a/util/classes/@dictionary/cumcoherence.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-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

--- a/util/classes/@dictionary/dictionary.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,111 +0,0 @@
-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

--- a/util/classes/@dictionary/mtimes.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,11 +0,0 @@
-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

--- a/util/classes/@dictionary/plotcumcoherencebounds.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,65 +0,0 @@
-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

--- a/util/classes/@dictionary/spark.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,35 +0,0 @@
-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

--- a/util/classes/dictionaryMatrices/dctmatrix.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-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

--- a/util/classes/dictionaryMatrices/grassmannian.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-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

--- a/util/classes/dictionaryMatrices/iterativeprojections.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,97 +0,0 @@
-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}');
-

--- a/util/classes/dictionaryMatrices/rotatematrix.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,114 +0,0 @@
-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

--- a/util/classes/dictionaryMatrices/shrinkgram.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,67 +0,0 @@
-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}');

--- a/util/classes/tests/test_buffer.m	Thu Apr 12 13:52:28 2012 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,48 +0,0 @@
-%% 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