diff examples/AudioInpainting/Audio_Declipping_Example.m @ 140:31d2864dfdd4 ivand_dev

Audio Impainting additional constraints with cvx added

--- a/examples/AudioInpainting/Audio_Declipping_Example.m	Thu Jul 21 16:37:14 2011 +0100
+++ b/examples/AudioInpainting/Audio_Declipping_Example.m	Mon Jul 25 17:27:05 2011 +0100
@@ -1,11 +1,16 @@
 %%  Audio Declipping Example
 %   
-%   CHANGE!!! This example is based on the experiment suggested by Professor Pierre
-%   Vandergheynst on the SMALL meeting in Villars.
-%   The idea behind is to use patches from source image as a dictionary in
-%   which we represent target image using matching pursuit algorithm.
-%   Calling Pierre_Problem function to get src image to be used as dictionary
-%   and target image to be represented using MP with 3 patches from source image
+%   Audio declipping is a problem proposed in Audio Inpaining Toolbox and
+%   in [2]. This is an example of solving the problem with fast omp using
+%   Gabor dictionary and ompGabor implemented in SMALLbox [1].
+%
+%   [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an 
+%   evaluation framework for sparse representations and dictionary 
+%   learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, 
+%   Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425
+%   [2] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, and M. D.
+%   Plumbley, “Audio Inpainting,” submitted to IEEE Trans. Audio, Speech,
+%   and Lang. Proc., 2011, http://hal.inria.fr/inria-00577079/en/.
 
 %
 %   Centre for Digital Music, Queen Mary, University of London.
@@ -20,88 +25,123 @@
 %%
 
 clear all;
+%   Defining the solvers to test in Audio declipping scenario
+
+% First solver omp2 - DCT+DST dictionary  with no additional constraints
+
+SMALL.solver(1) = SMALL_init_solver('ompbox', 'omp2', '', 0);
+SMALL.solver(1).add_constraints = 0;
+
+% Second solver omp2 - DCT+DST dictionary  with additional constraints
+
+SMALL.solver(2) = SMALL_init_solver('ompbox', 'omp2', '', 0);
+SMALL.solver(2).add_constraints = 1;
+
+% Third solver omp2 - Gabor dictionary with no additional constraints
+
+SMALL.solver(3) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0);
+SMALL.solver(3).add_constraints = 0;
+
+% Fourth solver omp2- Gabor dictionary with no additional constraints
+
+SMALL.solver(4) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0);
+SMALL.solver(4).add_constraints = 1;
 
 %   Defining the Problem structure
 
 SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2);
 
-% %   Show original image and image that is used as a dictionary
-% figure('Name', 'Original and Dictionary Image');
-% 
-% subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval);
-% title('Original Image');colormap(gray);axis off; axis image;
-% subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval);
-% title('Dictionary image:');colormap(gray);axis off; axis image;
-time=0;
-error2=0.001^2;
-coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
-
-for i=1:SMALL.Problem.n
+for idxSolver = 1:4
     
-    idx = find(SMALL.Problem.M(:,i));
-    if size(idx,1)==SMALL.Problem.m
-        continue
+    fprintf('\nStarting Audio Declipping of %s... \n', SMALL.Problem.name);
+    fprintf('\nClipping level %s... \n', SMALL.Problem.clippingLevel);
+    
+    start=cputime;
+    tStart=tic;
+    
+    error2=0.001^2;
+    coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
+    
+    
+    
+    for i = 1:SMALL.Problem.n
+        
+        idx = find(SMALL.Problem.M(:,i));
+        if size(idx,1)==SMALL.Problem.m
+            continue
+        end
+        Dict = SMALL.Problem.B(idx,:);
+        wDict = 1./sqrt(diag(Dict'*Dict));
+        
+        SMALL.Problem.A  = Dict*diag(wDict);
+        
+        SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
+        
+        
+        
+        % set solver parameters
+        
+        SMALL.solver(idxSolver).param=struct(...
+            'epsilon', error2*size(idx,1),...
+            'maxatoms', 128, ...
+            'profile', 'off');
+        
+        % Find solution
+        
+        SMALL.solver(idxSolver)=SMALL_solve(SMALL.Problem, SMALL.solver(idxSolver));
+        
+        % Refine solution with additional constraints for declipping scenario
+        
+        if (SMALL.solver(idxSolver).add_constraints)
+            SMALL.solver(idxSolver)=CVX_add_const_Audio_declipping(...
+                SMALL.Problem, SMALL.solver(idxSolver), i);
+        end
+        
+        %%
+        coeffFrames(:,i) = wDict .* SMALL.solver(idxSolver).solution;
+        
+        
     end
-    Dict = SMALL.Problem.B(idx,:);
-    wDict = 1./sqrt(diag(Dict'*Dict));
     
-    SMALL.Problem.A  = Dict*diag(wDict);
+    %%   Set reconstruction function
     
-    SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
+    SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
+    reconstructed=SMALL.Problem.reconstruct(coeffFrames);
+    SMALL.Problem = rmfield(SMALL.Problem, 'reconstruct');
+    tElapsed=toc(tStart);
     
+    SMALL.solver(idxSolver).time = cputime - start;
+    fprintf('Solver %s finished task in %2f seconds (cpu time). \n', ...
+        SMALL.solver(idxSolver).name, SMALL.solver(idxSolver).time);
+    fprintf('Solver %s finished task in %2f seconds (tic-toc time). \n', ...
+        SMALL.solver(idxSolver).name, tElapsed);
     
     
-    %   Defining the parameters sparse representation
-    SMALL.solver=SMALL_init_solver;
-    SMALL.solver.toolbox='ompbox';
-    SMALL.solver.name='omp2Gabor';
     
-    SMALL.solver.param=struct(...
-        'epsilon', error2*size(idx,1),...
-        'maxatoms', 128); 
+    %% Plot results
     
-    % Find solution
+    xClipped = SMALL.Problem.clipped;
+    xClean   = SMALL.Problem.original;
+    xEst1    = reconstructed.audioAllSamples;
+    xEst2    = reconstructed.audioOnlyClipped;
+    fs       = SMALL.Problem.fs;
     
-    SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver);
-    
-    
-    coeffFrames(:,i) = wDict .* SMALL.solver.solution;
-    time = time + SMALL.solver.time;
-    
-    
-      
+    figure
+    hold on
+    plot(xClipped,'r')
+    plot(xClean)
+    plot(xEst2,'--g')
+    plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
+    legend('Clipped','True solution','Estimate')
 end
 
-%%   Set reconstruction function
-
-SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
-reconstructed=SMALL.Problem.reconstruct(coeffFrames);
-
-
-
-%% Plot results
-
-xClipped = SMALL.Problem.clipped;
-xClean   = SMALL.Problem.original;
-xEst1    = reconstructed.audioAllSamples;
-xEst2    = reconstructed.audioOnlyClipped;
-fs       = SMALL.Problem.fs;
-
-figure
-hold on
-plot(xClipped,'r')
-plot(xClean)
-plot(xEst2,'--g')
-plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
-legend('Clipped','True solution','Estimate')
-
-% Normalized and save sounds
-normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
-L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
-xEst1 = xEst1(1:L)/normX;
-xEst2 = xEst2(1:L)/normX;
-xClipped = xClipped(1:L)/normX;
-xClean = xClean(1:L)/normX;
+% % Normalized and save sounds
+% normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
+% L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
+% xEst1 = xEst1(1:L)/normX;
+% xEst2 = xEst2(1:L)/normX;
+% xClipped = xClipped(1:L)/normX;
+% xClean = xClean(1:L)/normX;
 % wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']);
 % wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']);
 % wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);