ivan@136: %%  Audio Declipping Example
ivan@136: %   
ivan@136: %   CHANGE!!! This example is based on the experiment suggested by Professor Pierre
ivan@136: %   Vandergheynst on the SMALL meeting in Villars.
ivan@136: %   The idea behind is to use patches from source image as a dictionary in
ivan@136: %   which we represent target image using matching pursuit algorithm.
ivan@136: %   Calling Pierre_Problem function to get src image to be used as dictionary
ivan@136: %   and target image to be represented using MP with 3 patches from source image
ivan@136: 
ivan@136: %
ivan@136: %   Centre for Digital Music, Queen Mary, University of London.
ivan@136: %   This file copyright 2011 Ivan Damnjanovic.
ivan@136: %
ivan@136: %   This program is free software; you can redistribute it and/or
ivan@136: %   modify it under the terms of the GNU General Public License as
ivan@136: %   published by the Free Software Foundation; either version 2 of the
ivan@136: %   License, or (at your option) any later version.  See the file
ivan@136: %   COPYING included with this distribution for more information.
ivan@136: %
ivan@136: %%
ivan@136: 
ivan@136: clear all;
ivan@136: 
ivan@136: %   Defining the Problem structure
ivan@136: 
ivan@136: SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2);
ivan@136: 
ivan@136: % %   Show original image and image that is used as a dictionary
ivan@136: % figure('Name', 'Original and Dictionary Image');
ivan@136: % 
ivan@136: % subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval);
ivan@136: % title('Original Image');colormap(gray);axis off; axis image;
ivan@136: % subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval);
ivan@136: % title('Dictionary image:');colormap(gray);axis off; axis image;
ivan@136: time=0;
ivan@137: error2=0.001^2;
ivan@136: coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
ivan@136: 
ivan@136: for i=1:SMALL.Problem.n
ivan@136:     
ivan@136:     idx = find(SMALL.Problem.M(:,i));
ivan@137:     if size(idx,1)==SMALL.Problem.m
ivan@137:         continue
ivan@137:     end
ivan@137:     Dict = SMALL.Problem.B(idx,:);
ivan@137:     wDict = 1./sqrt(diag(Dict'*Dict));
ivan@137:     
ivan@137:     SMALL.Problem.A  = Dict*diag(wDict);
ivan@136:     
ivan@136:     SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
ivan@136:     
ivan@137:     
ivan@136:     
ivan@136:     %   Defining the parameters sparse representation
ivan@136:     SMALL.solver=SMALL_init_solver;
ivan@136:     SMALL.solver.toolbox='ompbox';
ivan@137:     SMALL.solver.name='omp2Gabor';
ivan@136:     
ivan@136:     SMALL.solver.param=struct(...
ivan@137:         'epsilon', error2*size(idx,1),...
ivan@137:         'maxatoms', 128); 
ivan@136:     
ivan@136:     % Find solution
ivan@136:     
ivan@136:     SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver);
ivan@136:     
ivan@136:     
ivan@137:     coeffFrames(:,i) = wDict .* SMALL.solver.solution;
ivan@136:     time = time + SMALL.solver.time;
ivan@136:     
ivan@136:     
ivan@136:       
ivan@136: end
ivan@136: 
ivan@136: %%   Set reconstruction function
ivan@136: 
ivan@136: SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
ivan@136: reconstructed=SMALL.Problem.reconstruct(coeffFrames);
ivan@136: 
ivan@136: 
ivan@137: 
ivan@137: %% Plot results
ivan@137: 
ivan@137: xClipped = SMALL.Problem.clipped;
ivan@137: xClean   = SMALL.Problem.original;
ivan@137: xEst1    = reconstructed.audioAllSamples;
ivan@137: xEst2    = reconstructed.audioOnlyClipped;
ivan@137: fs       = SMALL.Problem.fs;
ivan@137: 
ivan@137: figure
ivan@137: hold on
ivan@137: plot(xClipped,'r')
ivan@137: plot(xClean)
ivan@137: plot(xEst2,'--g')
ivan@137: plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
ivan@137: legend('Clipped','True solution','Estimate')
ivan@137: 
ivan@137: % Normalized and save sounds
ivan@137: normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
ivan@137: L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
ivan@137: xEst1 = xEst1(1:L)/normX;
ivan@137: xEst2 = xEst2(1:L)/normX;
ivan@137: xClipped = xClipped(1:L)/normX;
ivan@137: xClean = xClean(1:L)/normX;
ivan@137: wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']);
ivan@137: wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']);
ivan@137: wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);
ivan@137: wavwrite(xClean,fs,[expParam.destDir 'xClean.wav']);