annotate toolboxes/AudioInpaintingToolbox/Experiments/MissingSampleTopologyExperiment/MissingSampleTopologyExperiment.m @ 150:30872eb3d160 release_1.51

Merge from branch "ivand_dev"
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Tue, 26 Jul 2011 16:01:59 +0100
parents 56d719a5fd31
children
rev   line source
ivan@138 1 function MissingSampleTopologyExperiment(expParam)
ivan@138 2 % For a total number of missing samples C in a frame, create several
ivan@138 3 % configuration of B holes with length A, where A*B=C (i.e. the total
ivan@138 4 % number of missing samples is constant). Test several values of C, several
ivan@138 5 % solvers. For each C, test all possible combination of (A,B) such that
ivan@138 6 % A*B=C.
ivan@138 7 % Note that for each combination (A,B), a number of frames are tested at
ivan@138 8 % random and SNR results are then averaged.
ivan@138 9 %
ivan@138 10 % Usage: MissingSampleTopologyExperiment(expParam)
ivan@138 11 %
ivan@138 12 %
ivan@138 13 % Inputs:
ivan@138 14 % - expParam is an optional structure where the user can define
ivan@138 15 % the experiment parameters.
ivan@138 16 % - expParam.soundDir: path to sound directory. All the .wav files
ivan@138 17 % in this directory will be tested at random.
ivan@138 18 % - expParam.destDir: path to store the results.
ivan@138 19 % - expParam.N: frame length
ivan@138 20 % - expParam.NFramesPerHoleSize: number of frames to use for each
ivan@138 21 % testing configuration (A,B). Results are then averaged.
ivan@138 22 % - expParam.totalMissSamplesList: list of all tested values C for
ivan@138 23 % the total number of missing samples in a frame
ivan@138 24 % - expParam.solvers: list of solvers with their parameters
ivan@138 25 %
ivan@138 26 % -------------------
ivan@138 27 %
ivan@138 28 % Audio Inpainting toolbox
ivan@138 29 % Date: June 28, 2011
ivan@138 30 % By Valentin Emiya, Amir Adler, Maria Jafari
ivan@138 31 % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
ivan@138 32
ivan@138 33 if ~isdeployed
ivan@138 34 addpath('../../Problems/');
ivan@138 35 addpath('../../Solvers/');
ivan@138 36 addpath('../../Utils/');
ivan@138 37 addpath('../../Utils/dictionaries/');
ivan@138 38 addpath('../../Utils/evaluation/');
ivan@138 39 % addpath('../../Utils/TCPIP_SocketCom/');
ivan@138 40 % javaaddpath('../../Utils/TCPIP_SocketCom');
ivan@138 41 dbstop if error
ivan@138 42 close all
ivan@138 43 end
ivan@138 44
ivan@138 45 %% Set parameters
ivan@138 46 if nargin<1
ivan@138 47 expParam = [];
ivan@138 48 end
ivan@138 49 % Path to audio files
ivan@138 50 if ~isfield(expParam,'soundDir'),
ivan@138 51 expParam.soundDir = '../../Data/testSpeech8kHz_from16kHz/';
ivan@138 52 end
ivan@138 53 if ~isfield(expParam,'destDir'),
ivan@138 54 expParam.destDir = '../../tmp/missSampTopoExp/';
ivan@138 55 end
ivan@138 56 if ~exist(expParam.destDir,'dir')
ivan@138 57 mkdir(expParam.destDir)
ivan@138 58 end
ivan@138 59
ivan@138 60
ivan@138 61 % frame length
ivan@138 62 if ~isfield(expParam,'N'),
ivan@138 63 expParam.N = 512;
ivan@138 64 expParam.N = 256;
ivan@138 65 warning('AITB:N','Frame length=256 is used to have faster computations. Recommended frame length is 512 at 8kHz.');
ivan@138 66 end
ivan@138 67
ivan@138 68 % Number of random frames to test
ivan@138 69 if ~isfield(expParam,'NFramesPerHoleSize'),
ivan@138 70 expParam.NFramesPerHoleSize = 20;
ivan@138 71 warning('AITB:NFrames','expParam.NFramesPerHoleSize = 20 is used to have faster computations. Recommended value: several hundreds.');
ivan@138 72 end
ivan@138 73
ivan@138 74 % Number of missing samples: which numbers to test?
ivan@138 75 if ~isfield(expParam,'totalMissSamplesList')
ivan@138 76 expParam.totalMissSamplesList = [12,36,60,120,180,240];
ivan@138 77 expParam.totalMissSamplesList = [12,36];
ivan@138 78 warning('AITB:Miss','expParam.totalMissSamplesList = [12,36] is used to have faster computations. Recommended list: expParam.totalMissSamplesList = [12,36,60,120,180,240].');
ivan@138 79 end
ivan@138 80
ivan@138 81 % Choose the solver methods you would like to test: OMP, L1, Janssen
ivan@138 82 if ~isfield(expParam,'solvers'),
ivan@138 83 nSolver = 0;
ivan@138 84 nSolver = nSolver+1;
ivan@138 85 expParam.solvers(nSolver).name = 'OMP-G';
ivan@138 86 expParam.solvers(nSolver).inpaintFrame = @inpaintFrame_OMP_Gabor; % solver function
ivan@138 87 expParam.solvers(nSolver).param.N = expParam.N; % frame length
ivan@138 88 expParam.solvers(nSolver).param.OMPerr = 0.001;
ivan@138 89 expParam.solvers(nSolver).param.sparsityDegree = expParam.solvers(nSolver).param.N/4;
ivan@138 90 expParam.solvers(nSolver).param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
ivan@138 91 expParam.solvers(nSolver).param.redundancyFactor = 2; % Dictionary redundancy
ivan@138 92 expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138 93
ivan@138 94 nSolver = nSolver+1;
ivan@138 95 expParam.solvers(nSolver).name = 'Janssen';
ivan@138 96 expParam.solvers(nSolver).inpaintFrame = @inpaintFrame_janssenInterpolation; % solver function
ivan@138 97 expParam.solvers(nSolver).param.N = expParam.N; % frame length
ivan@138 98 end
ivan@138 99
ivan@138 100
ivan@138 101
ivan@138 102
ivan@138 103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ivan@138 104 soundDir = expParam.soundDir;
ivan@138 105 wavFiles = dir([soundDir '*.wav']);
ivan@138 106 wavFiles = arrayfun(@(x)[soundDir x.name],wavFiles,'UniformOutput',false);
ivan@138 107
ivan@138 108 %% Draw a list of random frames
ivan@138 109 % Choose an audio file at random
ivan@138 110 frameParam.kFrameFile = randi(length(wavFiles),expParam.NFramesPerHoleSize);
ivan@138 111
ivan@138 112 % For each audio file, find maximum mean energy among all frames
ivan@138 113 [dum fs] = wavread([soundDir wavFiles{1}],'size');
ivan@138 114 Ne = round(512/16000*fs);
ivan@138 115 E2m = zeros(length(wavFiles),1);
ivan@138 116 for kf = 1:length(wavFiles)
ivan@138 117 x=wavread(wavFiles{kf});
ivan@138 118 xm = filter(ones(Ne,1)/Ne,1,abs(x.^2));
ivan@138 119 E2m(kf) = 10*log10(max(xm));
ivan@138 120 end
ivan@138 121
ivan@138 122 % Choose the location of a frame at random, with a minimum energy
ivan@138 123 maxDiffE2m = 10;
ivan@138 124 frameParam.kFrameBegin = NaN(expParam.NFramesPerHoleSize,1);
ivan@138 125 for kf = 1:expParam.NFramesPerHoleSize
ivan@138 126 siz = wavread(wavFiles{frameParam.kFrameFile(kf)},'size');
ivan@138 127 while true
ivan@138 128 frameParam.kFrameBegin(kf) = randi(siz(1)-expParam.N+1);
ivan@138 129 x = wavread(wavFiles{frameParam.kFrameFile(kf)},[0,expParam.N-1]+frameParam.kFrameBegin(kf));
ivan@138 130 E2m0 = 10*log10(mean(abs(x.^2)));
ivan@138 131 if E2m(frameParam.kFrameFile(kf))-E2m0 <= maxDiffE2m
ivan@138 132 break
ivan@138 133 end
ivan@138 134 end
ivan@138 135 end
ivan@138 136
ivan@138 137 %% Test each number of missing samples
ivan@138 138 PerfRes = cell(length(expParam.totalMissSamplesList),length(expParam.solvers));
ivan@138 139 factorsToTest = cell(length(expParam.totalMissSamplesList),length(expParam.solvers));
ivan@138 140 outputFile = [expParam.destDir 'missSampTopoExp.mat'];
ivan@138 141 for kSolver = 1:length(expParam.solvers)
ivan@138 142 fprintf('\n ------ Solver: %s ------\n\n',...
ivan@138 143 expParam.solvers(kSolver).name);
ivan@138 144 for kMiss = 1:length(expParam.totalMissSamplesList)
ivan@138 145 NMissSamples = expParam.totalMissSamplesList(kMiss);
ivan@138 146 factorsToTest{kMiss} = allFactors(NMissSamples);
ivan@138 147 PerfRes{kMiss,kSolver} = zeros([length(factorsToTest{kMiss}),expParam.NFramesPerHoleSize]);
ivan@138 148 for kFactor = 1:length(factorsToTest{kMiss})
ivan@138 149 holeSize = factorsToTest{kMiss}(kFactor);
ivan@138 150 NHoles = NMissSamples/holeSize;
ivan@138 151 fprintf('%d %d-length holes (%d missing samples = %.1f%%)\n',...
ivan@138 152 NHoles,holeSize,NMissSamples,NMissSamples/expParam.N*100)
ivan@138 153 problemParameters.holeSize = holeSize;
ivan@138 154 problemParameters.NHoles = NHoles;
ivan@138 155 for kFrame = 1:expParam.NFramesPerHoleSize
ivan@138 156 %% load audio frame
ivan@138 157 xFrame = wavread(...
ivan@138 158 wavFiles{frameParam.kFrameFile(kFrame)},...
ivan@138 159 frameParam.kFrameBegin(kFrame)+[0,expParam.N-1]);
ivan@138 160
ivan@138 161 %% generate problem
ivan@138 162 [problemData, solutionData] = ...
ivan@138 163 generateMissingGroupsProblem(xFrame,problemParameters);
ivan@138 164
ivan@138 165 %% solve problem
ivan@138 166 xEst = ...
ivan@138 167 expParam.solvers(kSolver).inpaintFrame(...
ivan@138 168 problemData,...
ivan@138 169 expParam.solvers(kSolver).param);
ivan@138 170
ivan@138 171 %% compute and store performance
ivan@138 172 [SNRAll, SNRmiss] = ...
ivan@138 173 SNRInpaintingPerformance(...
ivan@138 174 solutionData.xClean,...
ivan@138 175 problemData.x,...
ivan@138 176 xEst,...
ivan@138 177 problemData.IMiss);
ivan@138 178 PerfRes{kMiss,kSolver}(kFactor,kFrame) = SNRmiss(2);
ivan@138 179
ivan@138 180 end
ivan@138 181 end
ivan@138 182 save(outputFile,'PerfRes','expParam');
ivan@138 183 end
ivan@138 184 end
ivan@138 185
ivan@138 186 figure
ivan@138 187 Nrows = floor(sqrt(length(expParam.solvers)));
ivan@138 188 Ncols = ceil(sqrt(length(expParam.solvers))/Nrows);
ivan@138 189 cmap = lines;
ivan@138 190 for kSolver = 1:length(expParam.solvers)
ivan@138 191 subplot(Nrows,Ncols,kSolver)
ivan@138 192 hold on,grid on
ivan@138 193 for kMiss = 1:length(expParam.totalMissSamplesList)
ivan@138 194 plot(factorsToTest{kMiss},mean(PerfRes{kMiss,kSolver},2),...
ivan@138 195 'color',cmap(kMiss,:));
ivan@138 196 end
ivan@138 197 title(expParam.solvers(kSolver).name)
ivan@138 198 end
ivan@138 199 return
ivan@138 200
ivan@138 201 function m = allFactors(n)
ivan@138 202 % Find the list of all factors (not only prime factors)
ivan@138 203
ivan@138 204 primeFactors = factor(n);
ivan@138 205
ivan@138 206 degrees = zeros(size(primeFactors));
ivan@138 207
ivan@138 208 for k=1:length(degrees)
ivan@138 209 degrees(k) = sum(primeFactors==primeFactors(k));
ivan@138 210 end
ivan@138 211
ivan@138 212 [primeFactors, I] = unique(primeFactors);
ivan@138 213 degrees = degrees(I);
ivan@138 214
ivan@138 215 D = (0:degrees(1)).';
ivan@138 216 for k=2:length(degrees)
ivan@138 217 Dk = ones(size(D,1),1)*(0:degrees(k));
ivan@138 218 D = [repmat(D,degrees(k)+1,1),Dk(:)];
ivan@138 219 end
ivan@138 220
ivan@138 221 m = unique(sort(prod((ones(size(D,1),1)*primeFactors).^D,2)));
ivan@138 222
ivan@138 223 return