comparison examples/SMALL_test_mocod.m @ 177:714fa7b8c1ad danieleb

added ramirez dl (to be completed) and MOCOD dictionary update

+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); 
+eta  	  = logspace(-2,2,10); 
+
+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
+nInitDicts  = length(initDicts);		%number of initial dictionaries
+nZetas = length(zeta);
+nEtas  = length(eta);
+
+SMALL.DL(nTrials,nInitDicts,nZetas,nEtas) = SMALL_init_DL(toolbox); %create dictionary learning structures
+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,...
+					'Tdata',nActiveAtoms,...
+					'dictsize',dictSize,...
+					'iternum',iterNum,...
+					'memusage','high',...
+					'solver',solver,...
+					'initdict',initDicts(iInitDicts),...
+					'zeta',zeta(iZetas),...
+					'eta',eta(iEtas));
+				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 for the various methods
+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
+
+save('MOCOD.mat')
+
+%% Plot results
+minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1)));	%lowe bound on coherence
+initDictsNames = {'Data','Gabor'};
+lineStyles     = {'k.-','r*-','b+-'};
+for iInitDict=1:nInitDicts
+	figure, hold on, grid on
+	title([initDictsNames{iInitDict} ' Initialisation']);
+	coherenceLevels = squeeze(mean(mu(:,iInitDict,:,:),1));
+	meanSNRs		= squeeze(mean(sr(:,iInitDict,:,:),1));
+	%stdSNRs				= squeeze(std(sr(:,iInitDict,iZetas,iEtas),0,1));
+	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')
+	
+	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')
+	
+	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