ivan@138: function declippingExperiment(expParam)
ivan@138: % Declip several sounds with different clipping levels, using several
ivan@138: % solvers.
ivan@138: %
ivan@138: % Usage: declippingExperiment(expParam)
ivan@138: %
ivan@138: %
ivan@138: % Inputs:
ivan@138: %          - expParam is an optional structure where the user can define
ivan@138: %          the experiment parameters.
ivan@138: %          - expParam.clippingScale: clipping values to test, as a vector
ivan@138: %           of real numbers in ]0,1[.
ivan@138: %          - expParam.soundDir: path to sound directory. All the .wav files
ivan@138: %          in this directory will be tested.
ivan@138: %          - expParam.destDir: path to store the results.
ivan@138: %          - expParam.solvers: list of solvers with their parameters
ivan@138: %
ivan@138: %
ivan@138: % -------------------
ivan@138: %
ivan@138: % Audio Inpainting toolbox
ivan@138: % Date: June 28, 2011
ivan@138: % By Valentin Emiya, Amir Adler, Maria Jafari
ivan@138: % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
ivan@138: 
ivan@138: if ~isdeployed
ivan@138:     addpath('../../Problems/');
ivan@138:     addpath('../../Solvers/');
ivan@138:     addpath('../../Utils/');
ivan@138:     addpath('../../Utils/dictionaries/');
ivan@138:     addpath('../../Utils/evaluation/');
ivan@138: %     addpath('../../Utils/TCPIP_SocketCom/');
ivan@138: %     javaaddpath('../../Utils/TCPIP_SocketCom');
ivan@138:     dbstop if error
ivan@138:     close all
ivan@138: end
ivan@138: 
ivan@138: if nargin<1
ivan@138:     expParam = [];
ivan@138: end
ivan@138: if ~isfield(expParam,'clippingScale'),
ivan@138:     expParam.clippingScale = 0.4:0.2:0.8;
ivan@138: end
ivan@138: if ~isfield(expParam,'soundDir'),
ivan@138:     expParam.soundDir = '../../Data/testSpeech8kHz_from16kHz/';
ivan@138:     expParam.soundDir = '../../Data/shortTest/';
ivan@138:     warning('AITB:soundDir','soundDir has only one sound to have faster computations. Recommended soundDir: ../../Data/testSpeech8kHz_from16kHz/');
ivan@138: end
ivan@138: if ~isfield(expParam,'destDir'),
ivan@138:     expParam.destDir = '../../tmp/declip/';
ivan@138: end
ivan@138: 
ivan@138: %% Set parameters
ivan@138: 
ivan@138: if ~isfield(expParam,'solvers'),
ivan@138:     % Choose the solver methods you would like to test: OMP, L1, Janssen
ivan@138:     warning('AITB:N','Frame length=256 is used to have faster computations. Recommended frame length is 512 at 8kHz.');
ivan@138:     warning('AITB:overlap','Overlap factor=2 is used to have faster computations. Recommended value: 4.');
ivan@138:     nSolver = 0;
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'OMP-C';
ivan@138:     expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138:     expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138:     expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138:     expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_OMP; % solver function
ivan@138:     expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138:     expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138:     expParam.solvers(nSolver).param.D_fun = @DCT_Dictionary; % Dictionary (function handle)
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138:     expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138:     expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138:     expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138:     expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'consOMP-C';
ivan@138:     expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138:     expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138:     expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138:     expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_consOMP; % solver function
ivan@138:     expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138:     expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138:     expParam.solvers(nSolver).param.D_fun = @DCT_Dictionary; % Dictionary (function handle)
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138:     expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138:     expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138:     expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138:     expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'OMP-G';
ivan@138:     expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138:     expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138:     expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138:     expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_OMP_Gabor; % solver function
ivan@138:     expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138:     expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138:     expParam.solvers(nSolver).param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138:     expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138:     expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138:     expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138:     expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'consOMP-G';
ivan@138:     expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138:     expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138:     expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138:     expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_consOMP_Gabor; % solver function
ivan@138:     expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138:     expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138:     expParam.solvers(nSolver).param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2
ivan@138:     expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138:     expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138:     expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138:     expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'Janssen';
ivan@138:     expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138:     expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_janssenInterpolation; % solver function
ivan@138:     expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138:     expParam.solvers(nSolver).param.N = 256;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138:     expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138:     expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138:     expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138: end
ivan@138: 
ivan@138: SNRClip = zeros(0,0,0);
ivan@138: fprintf('Folder %s\n',expParam.soundDir);
ivan@138: if ~exist(expParam.destDir,'dir')
ivan@138:     mkdir(expParam.destDir)
ivan@138: end
ivan@138: soundFiles = dir([expParam.soundDir '*.wav']);
ivan@138: 
ivan@138: for kf = 1:length(soundFiles)
ivan@138:     soundfile = [expParam.soundDir soundFiles(kf).name];
ivan@138:     fprintf(' File %s\n',soundfile);
ivan@138:     %% Read test signal
ivan@138:     [x fs] = wavread(soundfile);
ivan@138:     
ivan@138:     for kClip = 1:length(expParam.clippingScale)
ivan@138:         clippingLevel = expParam.clippingScale(kClip);
ivan@138:         fprintf('  Clip level %g\n',clippingLevel);
ivan@138:         
ivan@138:         %% Generate the problem
ivan@138:         [problemData, solutionData] = generateDeclippingProblem(x,clippingLevel);
ivan@138:         
ivan@138:         for nSolver = 1:length(expParam.solvers)
ivan@138:             %% Declip with solver
ivan@138:             solverParam = expParam.solvers(nSolver).param;
ivan@138:             [xEst1 xEst2] = expParam.solvers(nSolver).function(problemData,solverParam);
ivan@138:             
ivan@138:             %% compute performance
ivan@138:             L = length(xEst1);
ivan@138:             N = solverParam.N;
ivan@138:             [SNRAll, SNRmiss] = ...
ivan@138:                 SNRInpaintingPerformance(...
ivan@138:                 solutionData.xClean(N:L-N),...
ivan@138:                 problemData.x(N:L-N),...
ivan@138:                 xEst2(N:L-N),...
ivan@138:                 problemData.IMiss(N:L-N));
ivan@138:             SNRClip(kf,kClip,nSolver) = SNRmiss(2);
ivan@138:             
ivan@138:             % normalize and save both the reference and the estimates!
ivan@138:             normX = 1.1*max(abs([xEst1(:);xEst2(:);solutionData.xClean(:)]));
ivan@138:             
ivan@138:             L = min([length(xEst2),length(xEst1),length(solutionData.xClean),length(problemData.x)]);
ivan@138:             xEst1 = xEst1(1:L)/normX;
ivan@138:             xEst2 = xEst2(1:L)/normX;
ivan@138:             xClipped = problemData.x(1:L)/normX;
ivan@138:             xClean = solutionData.xClean(1:L)/normX;
ivan@138:             wavwrite(xEst1,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Est1',clippingLevel));
ivan@138:             wavwrite(xEst2,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Est2',clippingLevel));
ivan@138:             wavwrite(xClipped,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Clipped',clippingLevel));
ivan@138:             wavwrite(xClean,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Ref',clippingLevel));
ivan@138:             
ivan@138:             fprintf('\n');
ivan@138:             clear a xEst1 xEst2 xClipped xClean IClipped
ivan@138:             save([expParam.destDir 'clippingExp.mat']);
ivan@138:         end
ivan@138:     end
ivan@138: end
ivan@138: 
ivan@138: %% Plot results
ivan@138: averageSNR = squeeze(mean(SNRClip,1));
ivan@138: disp(averageSNR)
ivan@138: figure,
ivan@138: plot(averageSNR)
ivan@138: legend(arrayfun(@(x)x.name,expParam.solvers,'UniformOutput',false));
ivan@138: xlabel('Clipping level')
ivan@138: ylabel('SNR')
ivan@138: return