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