daniele@177: clc, clear, close all daniele@177: daniele@177: %% Parameteres daniele@177: nTrials = 10; %number of trials of the experiment daniele@177: daniele@177: % Dictionary learning parameters daniele@177: toolbox = 'TwoStepDL'; %dictionary learning toolbox daniele@177: dicUpdate = 'mocod'; %dictionary learning updates daniele@177: zeta = logspace(-2,2,10); daniele@177: eta = logspace(-2,2,10); daniele@177: daniele@177: iterNum = 20; %number of iterations daniele@177: epsilon = 1e-6; %tolerance level daniele@177: dictSize = 512; %number of atoms in the dictionary daniele@177: 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@177: initDicts = {[],gaborDict}; daniele@177: daniele@177: %% Generate audio approximation problem daniele@177: signal = buffer(signal,blockSize,blockSize*overlap,@rectwin); %buffer frames of audio into columns of the matrix S daniele@177: SMALL.Problem.b = signal.S; daniele@177: SMALL.Problem.b1 = SMALL.Problem.b; % copy signals from training set b to test set b1 (needed for later functions) daniele@177: daniele@177: % omp2 sparse representation solver daniele@177: ompParam = struct('X',SMALL.Problem.b,'epsilon',epsilon,'maxatoms',nActiveAtoms); %parameters daniele@177: solver = SMALL_init_solver('ompbox','omp2',ompParam,false); %solver structure daniele@177: daniele@177: %% Test daniele@177: nInitDicts = length(initDicts); %number of initial dictionaries daniele@177: nZetas = length(zeta); daniele@177: nEtas = length(eta); daniele@177: daniele@177: SMALL.DL(nTrials,nInitDicts,nZetas,nEtas) = SMALL_init_DL(toolbox); %create dictionary learning structures 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@177: struct('data',SMALL.Problem.b,... daniele@177: 'Tdata',nActiveAtoms,... daniele@177: 'dictsize',dictSize,... daniele@177: 'iternum',iterNum,... daniele@177: 'memusage','high',... daniele@177: 'solver',solver,... daniele@177: 'initdict',initDicts(iInitDicts),... daniele@177: 'zeta',zeta(iZetas),... daniele@177: 'eta',eta(iEtas)); 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@177: %% Evaluate coherence and snr of representation for the various methods daniele@177: sr = zeros(size(SMALL.DL)); %signal to noise ratio daniele@177: mu = zeros(iTrial,iInitDicts,iZetas,iEtas); %coherence daniele@177: 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@177: snr(SMALL.Problem.b,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: save('MOCOD.mat') daniele@177: daniele@177: %% Plot results daniele@177: minMu = sqrt((dictSize-blockSize)/(blockSize*(dictSize-1))); %lowe bound on coherence daniele@177: initDictsNames = {'Data','Gabor'}; daniele@177: lineStyles = {'k.-','r*-','b+-'}; daniele@177: for iInitDict=1:nInitDicts daniele@177: figure, hold on, grid on daniele@177: title([initDictsNames{iInitDict} ' Initialisation']); daniele@177: coherenceLevels = squeeze(mean(mu(:,iInitDict,:,:),1)); daniele@177: meanSNRs = squeeze(mean(sr(:,iInitDict,:,:),1)); daniele@177: %stdSNRs = squeeze(std(sr(:,iInitDict,iZetas,iEtas),0,1)); 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@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@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