ivan@138: function MissingSampleTopologyExperiment(expParam)
ivan@138: % For a total number of missing samples C in a frame, create several
ivan@138: % configuration of B holes with length A, where A*B=C (i.e. the total 
ivan@138: % number of missing samples is constant). Test several values of C, several
ivan@138: % solvers. For each C, test all possible combination of (A,B) such that
ivan@138: % A*B=C.
ivan@138: % Note that for each combination (A,B), a number of frames are tested at
ivan@138: % random and SNR results are then averaged.
ivan@138: %
ivan@138: % Usage: MissingSampleTopologyExperiment(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.soundDir: path to sound directory. All the .wav files
ivan@138: %          in this directory will be tested at random.
ivan@138: %          - expParam.destDir: path to store the results.
ivan@138: %          - expParam.N: frame length
ivan@138: %          - expParam.NFramesPerHoleSize: number of frames to use for each
ivan@138: %          testing configuration (A,B). Results are then averaged.
ivan@138: %          - expParam.totalMissSamplesList: list of all tested values C for
ivan@138: %          the total number of missing samples in a frame
ivan@138: %          - expParam.solvers: list of solvers with their parameters
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: %% Set parameters
ivan@138: if nargin<1
ivan@138:     expParam = [];
ivan@138: end
ivan@138: % Path to audio files
ivan@138: if ~isfield(expParam,'soundDir'),
ivan@138:     expParam.soundDir = '../../Data/testSpeech8kHz_from16kHz/';
ivan@138: end
ivan@138: if ~isfield(expParam,'destDir'),
ivan@138:     expParam.destDir = '../../tmp/missSampTopoExp/';
ivan@138: end
ivan@138: if ~exist(expParam.destDir,'dir')
ivan@138:     mkdir(expParam.destDir)
ivan@138: end
ivan@138: 
ivan@138: 
ivan@138: % frame length
ivan@138: if ~isfield(expParam,'N'),
ivan@138:     expParam.N = 512;
ivan@138:     expParam.N = 256;
ivan@138:     warning('AITB:N','Frame length=256 is used to have faster computations. Recommended frame length is 512 at 8kHz.');
ivan@138: end
ivan@138: 
ivan@138: % Number of random frames to test
ivan@138: if ~isfield(expParam,'NFramesPerHoleSize'),
ivan@138:     expParam.NFramesPerHoleSize = 20;
ivan@138:     warning('AITB:NFrames','expParam.NFramesPerHoleSize = 20 is used to have faster computations. Recommended value: several hundreds.');
ivan@138: end
ivan@138: 
ivan@138: % Number of missing samples: which numbers to test?
ivan@138: if ~isfield(expParam,'totalMissSamplesList')
ivan@138:     expParam.totalMissSamplesList = [12,36,60,120,180,240];
ivan@138:     expParam.totalMissSamplesList = [12,36];
ivan@138:     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: end
ivan@138: 
ivan@138: % Choose the solver methods you would like to test: OMP, L1, Janssen
ivan@138: if ~isfield(expParam,'solvers'),
ivan@138:     nSolver = 0;
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'OMP-G';
ivan@138:     expParam.solvers(nSolver).inpaintFrame = @inpaintFrame_OMP_Gabor; % solver function
ivan@138:     expParam.solvers(nSolver).param.N = expParam.N; % frame length
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.redundancyFactor = 2; % Dictionary redundancy
ivan@138:     expParam.solvers(nSolver).param.wa = @wRect; % Analysis window
ivan@138:     
ivan@138:     nSolver = nSolver+1;
ivan@138:     expParam.solvers(nSolver).name = 'Janssen';
ivan@138:     expParam.solvers(nSolver).inpaintFrame = @inpaintFrame_janssenInterpolation; % solver function
ivan@138:     expParam.solvers(nSolver).param.N = expParam.N; % frame length
ivan@138: end
ivan@138: 
ivan@138: 
ivan@138: 
ivan@138: 
ivan@138: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ivan@138: soundDir = expParam.soundDir;
ivan@138: wavFiles = dir([soundDir '*.wav']);
ivan@138: wavFiles = arrayfun(@(x)[soundDir x.name],wavFiles,'UniformOutput',false);
ivan@138: 
ivan@138: %% Draw a list of random frames
ivan@138: % Choose an audio file at random
ivan@138: frameParam.kFrameFile = randi(length(wavFiles),expParam.NFramesPerHoleSize);
ivan@138: 
ivan@138: % For each audio file, find maximum mean energy among all frames
ivan@138: [dum fs] = wavread([soundDir wavFiles{1}],'size');
ivan@138: Ne = round(512/16000*fs);
ivan@138: E2m = zeros(length(wavFiles),1);
ivan@138: for kf = 1:length(wavFiles)
ivan@138:     x=wavread(wavFiles{kf});
ivan@138:     xm = filter(ones(Ne,1)/Ne,1,abs(x.^2));
ivan@138:     E2m(kf) = 10*log10(max(xm));
ivan@138: end
ivan@138: 
ivan@138: % Choose the location of a frame at random, with a minimum energy
ivan@138: maxDiffE2m = 10;
ivan@138: frameParam.kFrameBegin = NaN(expParam.NFramesPerHoleSize,1);
ivan@138: for kf = 1:expParam.NFramesPerHoleSize
ivan@138:     siz = wavread(wavFiles{frameParam.kFrameFile(kf)},'size');
ivan@138:     while true
ivan@138:         frameParam.kFrameBegin(kf) = randi(siz(1)-expParam.N+1);
ivan@138:         x = wavread(wavFiles{frameParam.kFrameFile(kf)},[0,expParam.N-1]+frameParam.kFrameBegin(kf));
ivan@138:         E2m0 = 10*log10(mean(abs(x.^2)));
ivan@138:         if E2m(frameParam.kFrameFile(kf))-E2m0 <= maxDiffE2m
ivan@138:             break
ivan@138:         end
ivan@138:     end
ivan@138: end
ivan@138: 
ivan@138: %% Test each number of missing samples
ivan@138: PerfRes = cell(length(expParam.totalMissSamplesList),length(expParam.solvers));
ivan@138: factorsToTest = cell(length(expParam.totalMissSamplesList),length(expParam.solvers));
ivan@138: outputFile = [expParam.destDir 'missSampTopoExp.mat'];
ivan@138: for kSolver = 1:length(expParam.solvers)
ivan@138:     fprintf('\n ------ Solver: %s ------\n\n',...
ivan@138:         expParam.solvers(kSolver).name);
ivan@138:     for kMiss = 1:length(expParam.totalMissSamplesList)
ivan@138:         NMissSamples = expParam.totalMissSamplesList(kMiss);
ivan@138:         factorsToTest{kMiss} = allFactors(NMissSamples);
ivan@138:         PerfRes{kMiss,kSolver} = zeros([length(factorsToTest{kMiss}),expParam.NFramesPerHoleSize]);
ivan@138:         for kFactor = 1:length(factorsToTest{kMiss})
ivan@138:             holeSize = factorsToTest{kMiss}(kFactor);
ivan@138:             NHoles = NMissSamples/holeSize;
ivan@138:             fprintf('%d %d-length holes (%d missing samples = %.1f%%)\n',...
ivan@138:                 NHoles,holeSize,NMissSamples,NMissSamples/expParam.N*100)
ivan@138:             problemParameters.holeSize = holeSize;
ivan@138:             problemParameters.NHoles = NHoles;
ivan@138:             for kFrame = 1:expParam.NFramesPerHoleSize
ivan@138:                 %% load audio frame
ivan@138:                 xFrame = wavread(...
ivan@138:                     wavFiles{frameParam.kFrameFile(kFrame)},...
ivan@138:                     frameParam.kFrameBegin(kFrame)+[0,expParam.N-1]);
ivan@138:                 
ivan@138:                 %% generate problem
ivan@138:                 [problemData, solutionData] = ...
ivan@138:                     generateMissingGroupsProblem(xFrame,problemParameters);
ivan@138:                 
ivan@138:                 %% solve problem
ivan@138:                 xEst = ...
ivan@138:                     expParam.solvers(kSolver).inpaintFrame(...
ivan@138:                     problemData,...
ivan@138:                     expParam.solvers(kSolver).param);
ivan@138:                 
ivan@138:                 %% compute and store performance
ivan@138:                 [SNRAll, SNRmiss] = ...
ivan@138:                     SNRInpaintingPerformance(...
ivan@138:                     solutionData.xClean,...
ivan@138:                     problemData.x,...
ivan@138:                     xEst,...
ivan@138:                     problemData.IMiss);
ivan@138:                 PerfRes{kMiss,kSolver}(kFactor,kFrame) = SNRmiss(2);
ivan@138:                 
ivan@138:             end
ivan@138:         end
ivan@138:         save(outputFile,'PerfRes','expParam');
ivan@138:     end
ivan@138: end
ivan@138: 
ivan@138: figure
ivan@138: Nrows = floor(sqrt(length(expParam.solvers)));
ivan@138: Ncols = ceil(sqrt(length(expParam.solvers))/Nrows);
ivan@138: cmap = lines;
ivan@138: for kSolver = 1:length(expParam.solvers)
ivan@138:     subplot(Nrows,Ncols,kSolver)
ivan@138:     hold on,grid on
ivan@138:     for kMiss = 1:length(expParam.totalMissSamplesList)
ivan@138:         plot(factorsToTest{kMiss},mean(PerfRes{kMiss,kSolver},2),...
ivan@138:             'color',cmap(kMiss,:));
ivan@138:     end
ivan@138:     title(expParam.solvers(kSolver).name)
ivan@138: end
ivan@138: return
ivan@138: 
ivan@138: function m = allFactors(n)
ivan@138: % Find the list of all factors (not only prime factors)
ivan@138: 
ivan@138: primeFactors = factor(n);
ivan@138: 
ivan@138: degrees = zeros(size(primeFactors));
ivan@138: 
ivan@138: for k=1:length(degrees)
ivan@138:     degrees(k) = sum(primeFactors==primeFactors(k));
ivan@138: end
ivan@138: 
ivan@138: [primeFactors, I] = unique(primeFactors);
ivan@138: degrees = degrees(I);
ivan@138: 
ivan@138: D = (0:degrees(1)).';
ivan@138: for k=2:length(degrees)
ivan@138:     Dk = ones(size(D,1),1)*(0:degrees(k));
ivan@138:     D = [repmat(D,degrees(k)+1,1),Dk(:)];
ivan@138: end
ivan@138: 
ivan@138: m = unique(sort(prod((ones(size(D,1),1)*primeFactors).^D,2)));
ivan@138: 
ivan@138: return