# HG changeset patch # User luisf # Date 1334235203 -3600 # Node ID f0441ef3ae8fe88fab7ac387fd2e6db1ff3cd3b4 # Parent fd0b5d36f6ad544a7da2fec393276634e773f14f Updating this branch in order to reflect the default branch. diff -r fd0b5d36f6ad -r f0441ef3ae8f DL/dl_ramirez.m --- 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) 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f examples/SMALL_test_coherence.m --- 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') diff -r fd0b5d36f6ad -r f0441ef3ae8f examples/SMALL_test_coherence2.m --- 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)); diff -r fd0b5d36f6ad -r f0441ef3ae8f examples/SMALL_test_mocod.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@audio/audio.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@audio/buffer.m --- 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); diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@audio/plot.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@audio/unbuffer.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@dictionary/cumcoherence.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@dictionary/dictionary.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@dictionary/mtimes.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@dictionary/plotcumcoherencebounds.m --- 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/@dictionary/spark.m --- 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=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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/dictionaryMatrices/iterativeprojections.m --- 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 iItereps - 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 diff -r fd0b5d36f6ad -r f0441ef3ae8f util/classes/dictionaryMatrices/shrinkgram.m --- 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 iIterbuffer 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