changeset 136:1334d2302dd9 ivand_dev

Added Audio declipping problem (problem, reconstruct and example function)
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Thu, 14 Jul 2011 16:26:07 +0100 (2011-07-14)
parents 10343fb66448
children 9207d56c5547
files Problems/AudioDeclipping_reconstruct.m Problems/generateAudioDeclippingProblem.m examples/AudioInpainting/Audio_Declipping_Example.m util/Pierre_reconstruct.m
diffstat 4 files changed, 249 insertions(+), 33 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/AudioDeclipping_reconstruct.m	Thu Jul 14 16:26:07 2011 +0100
@@ -0,0 +1,50 @@
+function reconstructed=AudioDeclipping_reconstruct(y, Problem, SparseDict)
+%%  Audio declipping Problem reconstruction function
+%   
+%   This reconstruction function is using sparse representation y 
+%   in dictionary Problem.A to reconstruct the patches of the denoised
+%   image.
+
+%
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2009 Ivan Damnjanovic.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   published by the Free Software Foundation; either version 2 of the
+%   License, or (at your option) any later version.  See the file
+%   COPYING included with this distribution for more information.
+%%
+
+windowSize = Problem.windowSize;
+overlap = Problem.overlap;
+ws = Problem.ws(windowSize);
+wa = Problem.wa(windowSize);
+A = Problem.B;
+
+orig     = Problem.original;
+clipped  = Problem.clipped;
+clipMask = Problem.clipMask;
+
+% reconstruct audio frames
+
+xFrames = diag(ws)*(A*y);
+wNormFrames = (ws.*wa)'*ones(1,size(xFrames,2));
+
+%   overlap and add
+
+rec   = col2imstep(xFrames, size(clipped), [windowSize 1], [windowSize*overlap 1]);
+wNorm = col2imstep(wNormFrames, size(clipped), [windowSize 1], [windowSize*overlap 1]); 
+wNorm(find(wNorm==0)) = 1; 
+recN  = rec./wNorm;
+
+% change only clipped samples
+
+recSignal = orig.*double(~clipMask) + recN.*double(clipMask); 
+
+%% output structure image+psnr %%
+reconstructed.audioAllSamples  = recN;
+reconstructed.audioOnlyClipped = recSignal;
+[reconstructed.snrAll , reconstructed.snrMiss] = SNRInpaintingPerformance(orig, clipped, recSignal, clipMask, 1);
+
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/generateAudioDeclippingProblem.m	Thu Jul 14 16:26:07 2011 +0100
@@ -0,0 +1,120 @@
+function data = generateAudioDeclippingProblem(soundfile, clippingLevel, windowSize, overlap, wa, ws, wd, Dict_fun, redundancyFactor)
+%%  Generate Audio Declipping Problem
+%   
+%   CHANGE!!!generateAMT_Learning_Problem is a part of the SMALLbox and generates
+%   a problem that can be used for comparison of Dictionary Learning/Sparse
+%   Representation techniques in automatic music transcription scenario.
+%   The function prompts a user for an audio file (mid, wav, mat) reads it
+%   and generates a spectrogram given fft size (default nfft=4096), analysis
+%   window size (windowSize=2822), and analysis window overlap (overlap =
+%   0.5).
+%   
+%   The output of the function is stucture with following fields:
+%       b - matrix with magnitudes of the spectrogram
+%       f - vector of frequencies at wihch spectrogram is computed
+%       windowSize - analysis window size
+%       overlap - analysis window overlap
+%       fs - sampling frequency
+%       m - number of frequenciy points in spectrogram
+%       n - number of time points in the spectrogram
+%       p - number of dictionary elements to be learned (eg 88 for piano)
+%       notesOriginal - notes of the original audio to be used for
+%                       comparison (if midi of the original exists)
+%       name - name of the audio file to transcribe
+
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2011 Ivan Damnjanovic.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   published by the Free Software Foundation; either version 2 of the
+%   License, or (at your option) any later version.  See the file
+%   COPYING included with this distribution for more information.
+%  
+%%
+FS=filesep;
+TMPpath=pwd;
+
+if ~ exist( 'soundfile', 'var' ) || isempty(soundfile)
+    %ask for file name 
+    [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+    cd([pathstr1,FS,'data',FS,'audio']);
+    [filename,pathname] = uigetfile({'*.mat; *.mid; *.wav'},'Select a file to transcribe');
+    [pathstr, name, ext, versn] = fileparts(filename);
+    data.name=name;
+
+    if strcmp(ext,'.mid')
+        midi=readmidi(filename);
+%         data.notesOriginal=midiInfo(midi);
+        y=midi2audio(midi);
+        wavwrite(y, 44100, 16, 'temp.wav');
+        [x.signal, x.fs, x.nbits]=wavread('temp.wav');
+        delete('temp.wav');
+    elseif strcmp(ext,'.wav')
+%         cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']);
+%         filename1=[name, '.mid'];
+%         if exist(filename1, 'file')
+%             midi=readmidi(filename1);
+%             data.notesOriginal=midiInfo(midi);
+%         end
+        cd([pathstr1,FS, 'data', FS, 'audio', FS, 'wav']);
+        [x.signal, x.fs, x.nbits]=wavread(filename);
+    else
+%         cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']);
+%         filename1=[name, '.mid'];
+%         if exist(filename1, 'file')
+%             midi=readmidi(filename1);
+%             data.notesOriginal=midiInfo(midi);
+%         end
+        cd([pathstr1,FS, 'data', FS, 'audio', FS, 'mat']);
+        x=load([pathname,filename]);
+    end
+else
+    [x.signal, x.fs, x.nbits]=wavread(soundfile);
+    [pathstr, name, ext, versn] = fileparts(soundfile);
+    data.name=name;
+end
+
+if ~ exist( 'clippingLevel', 'var' ) || isempty(clippingLevel), clippingLevel = 0.6; end
+if ~ exist( 'windowSize', 'var' ) || isempty(windowSize), windowSize = 256; end
+if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 0.5; end
+if ~ exist( 'wa', 'var' ) || isempty(wa), wa = @wRect; end % Analysis window
+if ~ exist( 'ws', 'var' ) || isempty(ws), ws = @wSine; end % Synthesis window
+if ~ exist( 'wd', 'var' ) || isempty(wd), wd = @wRect; end % Weighting window for dictionary atoms
+
+%% preparing signal
+
+[problemData, solutionData] = generateDeclippingProblem(x.signal,clippingLevel);
+
+x_clip = im2colstep(problemData.x,[windowSize 1],[overlap*windowSize 1]);
+x_clip= diag(wa(windowSize)) * x_clip;
+blkMask=im2colstep(double(~problemData.IMiss),[256 1],[128 1]);
+
+%% Building dictionary
+if ~exist( 'redundancyFactor', 'var' ) || isempty(redundancyFactor), redundancyFactor = 2; end % Weighting window for dictionary atoms
+if exist('Dict_fun', 'var')&&~isempty(Dict_fun)
+    param=struct('N', windowSize, 'redundancyFactor', redundancyFactor, 'wd', wd);
+	data.B = Dict_fun(param);
+end
+
+data.b = x_clip;
+data.M = blkMask;
+data.original = solutionData.xClean;
+data.clipped = problemData.x;
+data.clipMask = problemData.IMiss;
+data.clippingLevel = clippingLevel;
+data.windowSize = windowSize;
+data.overlap = overlap;
+data.ws = ws;
+data.wa = wa;
+data.wd = wd;
+
+data.fs = x.fs;
+data.nbits = x.nbits;
+
+[data.m, data.n] = size(x_clip);
+data.p = windowSize*redundancyFactor; %number of dictionary elements 
+
+cd(TMPpath);
+
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/AudioInpainting/Audio_Declipping_Example.m	Thu Jul 14 16:26:07 2011 +0100
@@ -0,0 +1,79 @@
+%%  Audio Declipping Example
+%   
+%   CHANGE!!! This example is based on the experiment suggested by Professor Pierre
+%   Vandergheynst on the SMALL meeting in Villars.
+%   The idea behind is to use patches from source image as a dictionary in
+%   which we represent target image using matching pursuit algorithm.
+%   Calling Pierre_Problem function to get src image to be used as dictionary
+%   and target image to be represented using MP with 3 patches from source image
+
+%
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2011 Ivan Damnjanovic.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   published by the Free Software Foundation; either version 2 of the
+%   License, or (at your option) any later version.  See the file
+%   COPYING included with this distribution for more information.
+%
+%%
+
+clear all;
+
+%   Defining the Problem structure
+
+SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2);
+
+% %   Show original image and image that is used as a dictionary
+% figure('Name', 'Original and Dictionary Image');
+% 
+% subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval);
+% title('Original Image');colormap(gray);axis off; axis image;
+% subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval);
+% title('Dictionary image:');colormap(gray);axis off; axis image;
+time=0;
+coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
+
+for i=1:SMALL.Problem.n
+    
+    idx = find(SMALL.Problem.M(:,i));
+    SMALL.Problem.A  = SMALL.Problem.B(idx,:);
+    
+    SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
+    
+   
+    
+    %   Defining the parameters sparse representation
+    SMALL.solver=SMALL_init_solver;
+    SMALL.solver.toolbox='ompbox';
+    SMALL.solver.name='omp2';
+    
+    SMALL.solver.param=struct(...
+        'epsilon', 0.001,...
+        'maxatoms', 64); 
+    
+    % Find solution
+    
+    SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver);
+    
+    
+    coeffFrames(:,i) = SMALL.solver.solution;
+    time = time + SMALL.solver.time;
+    
+    
+      
+end
+
+%%   Set reconstruction function
+
+SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
+reconstructed=SMALL.Problem.reconstruct(coeffFrames);
+
+%%  plot time and psnr given dictionary size %%
+figure('Name', 'time and psnr');
+
+subplot(1,2,1); plot(dictsize(1,:), time(1,:), 'ro-');
+title('Time vs number of source image patches used');
+subplot(1,2,2); plot(dictsize(1,:), psnr(1,:), 'b*-');
+title('PSNR vs number of source image patches used');
\ No newline at end of file
--- a/util/Pierre_reconstruct.m	Mon Jul 11 13:42:00 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-function reconstructed=Pierre_reconstruct(y, Problem)
-%%  Pierre Villars Example - reconstruction function
-%
-%   using sparse representation y in dictionary Problem.A reconstruct the
-%   patches from the target image
-%   This example is based on the experiment suggested by Professor Pierre
-%   Vandergheynst on the SMALL meeting in Villars.
-
-%
-%   Centre for Digital Music, Queen Mary, University of London.
-%   This file copyright 2009 Ivan Damnjanovic.
-%
-%   This program is free software; you can redistribute it and/or
-%   modify it under the terms of the GNU General Public License as
-%   published by the Free Software Foundation; either version 2 of the
-%   License, or (at your option) any later version.  See the file
-%   COPYING included with this distribution for more information.
-%   
-%%
-imout=Problem.A*y;
-
-%   combine the patches into reconstructed image
-
-im=col2imstep(imout,size(Problem.imageTrg),Problem.blocksize,Problem.blocksize);
-
-%   bound the pixel values to [0,255] range 
-im(im<0)=0;
-im(im>255)=255;
-
-%% output structure image+psnr %%
-reconstructed.image=im;
-reconstructed.psnr = 20*log10(Problem.maxval * sqrt(numel(Problem.imageTrg(:))) / norm(Problem.imageTrg(:)-im(:)));
-end
\ No newline at end of file