annotate toolboxes/AudioInpaintingToolbox/Experiments/DeclippingExperiment/declippingExperiment.m @ 171:e8428989412f danieleb

Added dictionary decorrelation functions and test script for Letters paper.
author Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
date Thu, 06 Oct 2011 14:33:52 +0100
parents 56d719a5fd31
children
rev   line source
ivan@138 1 function declippingExperiment(expParam)
ivan@138 2 % Declip several sounds with different clipping levels, using several
ivan@138 3 % solvers.
ivan@138 4 %
ivan@138 5 % Usage: declippingExperiment(expParam)
ivan@138 6 %
ivan@138 7 %
ivan@138 8 % Inputs:
ivan@138 9 % - expParam is an optional structure where the user can define
ivan@138 10 % the experiment parameters.
ivan@138 11 % - expParam.clippingScale: clipping values to test, as a vector
ivan@138 12 % of real numbers in ]0,1[.
ivan@138 13 % - expParam.soundDir: path to sound directory. All the .wav files
ivan@138 14 % in this directory will be tested.
ivan@138 15 % - expParam.destDir: path to store the results.
ivan@138 16 % - expParam.solvers: list of solvers with their parameters
ivan@138 17 %
ivan@138 18 %
ivan@138 19 % -------------------
ivan@138 20 %
ivan@138 21 % Audio Inpainting toolbox
ivan@138 22 % Date: June 28, 2011
ivan@138 23 % By Valentin Emiya, Amir Adler, Maria Jafari
ivan@138 24 % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
ivan@138 25
ivan@138 26 if ~isdeployed
ivan@138 27 addpath('../../Problems/');
ivan@138 28 addpath('../../Solvers/');
ivan@138 29 addpath('../../Utils/');
ivan@138 30 addpath('../../Utils/dictionaries/');
ivan@138 31 addpath('../../Utils/evaluation/');
ivan@138 32 % addpath('../../Utils/TCPIP_SocketCom/');
ivan@138 33 % javaaddpath('../../Utils/TCPIP_SocketCom');
ivan@138 34 dbstop if error
ivan@138 35 close all
ivan@138 36 end
ivan@138 37
ivan@138 38 if nargin<1
ivan@138 39 expParam = [];
ivan@138 40 end
ivan@138 41 if ~isfield(expParam,'clippingScale'),
ivan@138 42 expParam.clippingScale = 0.4:0.2:0.8;
ivan@138 43 end
ivan@138 44 if ~isfield(expParam,'soundDir'),
ivan@138 45 expParam.soundDir = '../../Data/testSpeech8kHz_from16kHz/';
ivan@138 46 expParam.soundDir = '../../Data/shortTest/';
ivan@138 47 warning('AITB:soundDir','soundDir has only one sound to have faster computations. Recommended soundDir: ../../Data/testSpeech8kHz_from16kHz/');
ivan@138 48 end
ivan@138 49 if ~isfield(expParam,'destDir'),
ivan@138 50 expParam.destDir = '../../tmp/declip/';
ivan@138 51 end
ivan@138 52
ivan@138 53 %% Set parameters
ivan@138 54
ivan@138 55 if ~isfield(expParam,'solvers'),
ivan@138 56 % Choose the solver methods you would like to test: OMP, L1, Janssen
ivan@138 57 warning('AITB:N','Frame length=256 is used to have faster computations. Recommended frame length is 512 at 8kHz.');
ivan@138 58 warning('AITB:overlap','Overlap factor=2 is used to have faster computations. Recommended value: 4.');
ivan@138 59 nSolver = 0;
ivan@138 60
ivan@138 61 nSolver = nSolver+1;
ivan@138 62 expParam.solvers(nSolver).name = 'OMP-C';
ivan@138 63 expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138 64 expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138 65 expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138 66 expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_OMP; % solver function
ivan@138 67 expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138 68 expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138 69 expParam.solvers(nSolver).param.D_fun = @DCT_Dictionary; % Dictionary (function handle)
ivan@138 70 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138 71 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138 72 expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138 73 expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138 74 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 75 expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138 76 expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138 77 expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138 78
ivan@138 79 nSolver = nSolver+1;
ivan@138 80 expParam.solvers(nSolver).name = 'consOMP-C';
ivan@138 81 expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138 82 expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138 83 expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138 84 expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_consOMP; % solver function
ivan@138 85 expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138 86 expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138 87 expParam.solvers(nSolver).param.D_fun = @DCT_Dictionary; % Dictionary (function handle)
ivan@138 88 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138 89 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138 90 expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138 91 expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138 92 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 93 expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138 94 expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138 95 expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138 96
ivan@138 97 nSolver = nSolver+1;
ivan@138 98 expParam.solvers(nSolver).name = 'OMP-G';
ivan@138 99 expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138 100 expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138 101 expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138 102 expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_OMP_Gabor; % solver function
ivan@138 103 expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138 104 expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138 105 expParam.solvers(nSolver).param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
ivan@138 106 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138 107 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2;
ivan@138 108 expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138 109 expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138 110 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 111 expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138 112 expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138 113 expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138 114
ivan@138 115 nSolver = nSolver+1;
ivan@138 116 expParam.solvers(nSolver).name = 'consOMP-G';
ivan@138 117 expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138 118 expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138 119 expParam.solvers(nSolver).param.N = 256; % frame length
ivan@138 120 expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_consOMP_Gabor; % solver function
ivan@138 121 expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138 122 expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138 123 expParam.solvers(nSolver).param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
ivan@138 124 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138 125 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2
ivan@138 126 expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138 127 expParam.solvers(nSolver).param.wd = @wRect; % Weighting window for dictionary atoms
ivan@138 128 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 129 expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138 130 expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138 131 expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138 132
ivan@138 133 nSolver = nSolver+1;
ivan@138 134 expParam.solvers(nSolver).name = 'Janssen';
ivan@138 135 expParam.solvers(nSolver).function = @inpaintSignal_IndependentProcessingOfFrames;
ivan@138 136 expParam.solvers(nSolver).param.inpaintFrame = @inpaintFrame_janssenInterpolation; % solver function
ivan@138 137 expParam.solvers(nSolver).param.N = 512; % frame length
ivan@138 138 expParam.solvers(nSolver).param.N = 256;
ivan@138 139 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 4;
ivan@138 140 expParam.solvers(nSolver).param.OLA_frameOverlapFactor = 2
ivan@138 141 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 142 expParam.solvers(nSolver).param.OLA_ws = @wSine; % Synthesis window
ivan@138 143 expParam.solvers(nSolver).param.SKIP_CLEAN_FRAMES = true; % do not process frames where there is no missing samples
ivan@138 144 expParam.solvers(nSolver).param.MULTITHREAD_FRAME_PROCESSING = false; % not implemented yet
ivan@138 145 end
ivan@138 146
ivan@138 147 SNRClip = zeros(0,0,0);
ivan@138 148 fprintf('Folder %s\n',expParam.soundDir);
ivan@138 149 if ~exist(expParam.destDir,'dir')
ivan@138 150 mkdir(expParam.destDir)
ivan@138 151 end
ivan@138 152 soundFiles = dir([expParam.soundDir '*.wav']);
ivan@138 153
ivan@138 154 for kf = 1:length(soundFiles)
ivan@138 155 soundfile = [expParam.soundDir soundFiles(kf).name];
ivan@138 156 fprintf(' File %s\n',soundfile);
ivan@138 157 %% Read test signal
ivan@138 158 [x fs] = wavread(soundfile);
ivan@138 159
ivan@138 160 for kClip = 1:length(expParam.clippingScale)
ivan@138 161 clippingLevel = expParam.clippingScale(kClip);
ivan@138 162 fprintf(' Clip level %g\n',clippingLevel);
ivan@138 163
ivan@138 164 %% Generate the problem
ivan@138 165 [problemData, solutionData] = generateDeclippingProblem(x,clippingLevel);
ivan@138 166
ivan@138 167 for nSolver = 1:length(expParam.solvers)
ivan@138 168 %% Declip with solver
ivan@138 169 solverParam = expParam.solvers(nSolver).param;
ivan@138 170 [xEst1 xEst2] = expParam.solvers(nSolver).function(problemData,solverParam);
ivan@138 171
ivan@138 172 %% compute performance
ivan@138 173 L = length(xEst1);
ivan@138 174 N = solverParam.N;
ivan@138 175 [SNRAll, SNRmiss] = ...
ivan@138 176 SNRInpaintingPerformance(...
ivan@138 177 solutionData.xClean(N:L-N),...
ivan@138 178 problemData.x(N:L-N),...
ivan@138 179 xEst2(N:L-N),...
ivan@138 180 problemData.IMiss(N:L-N));
ivan@138 181 SNRClip(kf,kClip,nSolver) = SNRmiss(2);
ivan@138 182
ivan@138 183 % normalize and save both the reference and the estimates!
ivan@138 184 normX = 1.1*max(abs([xEst1(:);xEst2(:);solutionData.xClean(:)]));
ivan@138 185
ivan@138 186 L = min([length(xEst2),length(xEst1),length(solutionData.xClean),length(problemData.x)]);
ivan@138 187 xEst1 = xEst1(1:L)/normX;
ivan@138 188 xEst2 = xEst2(1:L)/normX;
ivan@138 189 xClipped = problemData.x(1:L)/normX;
ivan@138 190 xClean = solutionData.xClean(1:L)/normX;
ivan@138 191 wavwrite(xEst1,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Est1',clippingLevel));
ivan@138 192 wavwrite(xEst2,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Est2',clippingLevel));
ivan@138 193 wavwrite(xClipped,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Clipped',clippingLevel));
ivan@138 194 wavwrite(xClean,fs,sprintf('%s%s%s%g.wav',expParam.destDir,soundFiles(kf).name(1:end-4),'Ref',clippingLevel));
ivan@138 195
ivan@138 196 fprintf('\n');
ivan@138 197 clear a xEst1 xEst2 xClipped xClean IClipped
ivan@138 198 save([expParam.destDir 'clippingExp.mat']);
ivan@138 199 end
ivan@138 200 end
ivan@138 201 end
ivan@138 202
ivan@138 203 %% Plot results
ivan@138 204 averageSNR = squeeze(mean(SNRClip,1));
ivan@138 205 disp(averageSNR)
ivan@138 206 figure,
ivan@138 207 plot(averageSNR)
ivan@138 208 legend(arrayfun(@(x)x.name,expParam.solvers,'UniformOutput',false));
ivan@138 209 xlabel('Clipping level')
ivan@138 210 ylabel('SNR')
ivan@138 211 return