annotate examples/AudioInpainting/Audio_Declipping_Example.m @ 155:b14209313ba4 ivand_dev

Integration of Majorization Minimisation Dictionary Learning
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Mon, 22 Aug 2011 11:46:35 +0100
parents 0de08f68256b
children 8b3c71bb44eb
rev   line source
ivan@136 1 %% Audio Declipping Example
ivan@136 2 %
ivan@140 3 % Audio declipping is a problem proposed in Audio Inpaining Toolbox and
ivan@140 4 % in [2]. This is an example of solving the problem with fast omp using
ivan@140 5 % Gabor dictionary and ompGabor implemented in SMALLbox [1].
ivan@140 6 %
ivan@140 7 % [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an
ivan@140 8 % evaluation framework for sparse representations and dictionary
ivan@140 9 % learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010,
ivan@140 10 % Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425
ivan@140 11 % [2] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, and M. D.
ivan@140 12 % Plumbley, “Audio Inpainting,” submitted to IEEE Trans. Audio, Speech,
ivan@140 13 % and Lang. Proc., 2011, http://hal.inria.fr/inria-00577079/en/.
ivan@136 14
ivan@136 15 %
ivan@136 16 % Centre for Digital Music, Queen Mary, University of London.
ivan@136 17 % This file copyright 2011 Ivan Damnjanovic.
ivan@136 18 %
ivan@136 19 % This program is free software; you can redistribute it and/or
ivan@136 20 % modify it under the terms of the GNU General Public License as
ivan@136 21 % published by the Free Software Foundation; either version 2 of the
ivan@136 22 % License, or (at your option) any later version. See the file
ivan@136 23 % COPYING included with this distribution for more information.
ivan@136 24 %
ivan@136 25 %%
ivan@136 26
ivan@136 27 clear all;
ivan@140 28 % Defining the solvers to test in Audio declipping scenario
ivan@140 29
ivan@140 30 % First solver omp2 - DCT+DST dictionary with no additional constraints
ivan@140 31
ivan@140 32 SMALL.solver(1) = SMALL_init_solver('ompbox', 'omp2', '', 0);
ivan@140 33 SMALL.solver(1).add_constraints = 0;
ivan@140 34
ivan@140 35 % Second solver omp2 - DCT+DST dictionary with additional constraints
ivan@140 36
ivan@140 37 SMALL.solver(2) = SMALL_init_solver('ompbox', 'omp2', '', 0);
ivan@140 38 SMALL.solver(2).add_constraints = 1;
ivan@140 39
ivan@140 40 % Third solver omp2 - Gabor dictionary with no additional constraints
ivan@140 41
ivan@140 42 SMALL.solver(3) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0);
ivan@140 43 SMALL.solver(3).add_constraints = 0;
ivan@140 44
ivan@140 45 % Fourth solver omp2- Gabor dictionary with no additional constraints
ivan@140 46
ivan@140 47 SMALL.solver(4) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0);
ivan@140 48 SMALL.solver(4).add_constraints = 1;
ivan@136 49
ivan@136 50 % Defining the Problem structure
ivan@136 51
ivan@155 52 SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2);
ivan@136 53
ivan@140 54 for idxSolver = 1:4
ivan@136 55
ivan@140 56 fprintf('\nStarting Audio Declipping of %s... \n', SMALL.Problem.name);
ivan@140 57 fprintf('\nClipping level %s... \n', SMALL.Problem.clippingLevel);
ivan@140 58
ivan@140 59 start=cputime;
ivan@140 60 tStart=tic;
ivan@140 61
ivan@140 62 error2=0.001^2;
ivan@140 63 coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
ivan@140 64
ivan@140 65
ivan@140 66
ivan@140 67 for i = 1:SMALL.Problem.n
ivan@140 68
ivan@140 69 idx = find(SMALL.Problem.M(:,i));
ivan@140 70 if size(idx,1)==SMALL.Problem.m
ivan@140 71 continue
ivan@140 72 end
ivan@140 73 Dict = SMALL.Problem.B(idx,:);
ivan@140 74 wDict = 1./sqrt(diag(Dict'*Dict));
ivan@140 75
ivan@140 76 SMALL.Problem.A = Dict*diag(wDict);
ivan@140 77
ivan@140 78 SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
ivan@140 79
ivan@140 80
ivan@140 81
ivan@140 82 % set solver parameters
ivan@140 83
ivan@140 84 SMALL.solver(idxSolver).param=struct(...
ivan@140 85 'epsilon', error2*size(idx,1),...
ivan@140 86 'maxatoms', 128, ...
ivan@140 87 'profile', 'off');
ivan@140 88
ivan@140 89 % Find solution
ivan@140 90
ivan@140 91 SMALL.solver(idxSolver)=SMALL_solve(SMALL.Problem, SMALL.solver(idxSolver));
ivan@140 92
ivan@140 93 % Refine solution with additional constraints for declipping scenario
ivan@140 94
ivan@140 95 if (SMALL.solver(idxSolver).add_constraints)
ivan@140 96 SMALL.solver(idxSolver)=CVX_add_const_Audio_declipping(...
ivan@140 97 SMALL.Problem, SMALL.solver(idxSolver), i);
ivan@140 98 end
ivan@140 99
ivan@140 100 %%
ivan@140 101 coeffFrames(:,i) = wDict .* SMALL.solver(idxSolver).solution;
ivan@140 102
ivan@140 103
ivan@137 104 end
ivan@137 105
ivan@140 106 %% Set reconstruction function
ivan@136 107
ivan@140 108 SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
ivan@140 109 reconstructed=SMALL.Problem.reconstruct(coeffFrames);
ivan@140 110 SMALL.Problem = rmfield(SMALL.Problem, 'reconstruct');
ivan@140 111 tElapsed=toc(tStart);
ivan@136 112
ivan@140 113 SMALL.solver(idxSolver).time = cputime - start;
ivan@140 114 fprintf('Solver %s finished task in %2f seconds (cpu time). \n', ...
ivan@140 115 SMALL.solver(idxSolver).name, SMALL.solver(idxSolver).time);
ivan@140 116 fprintf('Solver %s finished task in %2f seconds (tic-toc time). \n', ...
ivan@140 117 SMALL.solver(idxSolver).name, tElapsed);
ivan@137 118
ivan@136 119
ivan@136 120
ivan@140 121 %% Plot results
ivan@136 122
ivan@140 123 xClipped = SMALL.Problem.clipped;
ivan@140 124 xClean = SMALL.Problem.original;
ivan@140 125 xEst1 = reconstructed.audioAllSamples;
ivan@140 126 xEst2 = reconstructed.audioOnlyClipped;
ivan@140 127 fs = SMALL.Problem.fs;
ivan@136 128
ivan@140 129 figure
ivan@140 130 hold on
ivan@140 131 plot(xClipped,'r')
ivan@140 132 plot(xClean)
ivan@140 133 plot(xEst2,'--g')
ivan@140 134 plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
ivan@140 135 legend('Clipped','True solution','Estimate')
ivan@136 136 end
ivan@136 137
ivan@140 138 % % Normalized and save sounds
ivan@140 139 % normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
ivan@140 140 % L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
ivan@140 141 % xEst1 = xEst1(1:L)/normX;
ivan@140 142 % xEst2 = xEst2(1:L)/normX;
ivan@140 143 % xClipped = xClipped(1:L)/normX;
ivan@140 144 % xClean = xClean(1:L)/normX;
ivan@139 145 % wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']);
ivan@139 146 % wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']);
ivan@139 147 % wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);
ivan@139 148 % wavwrite(xClean,fs,[expParam.destDir 'xClean.wav']);