ivan@136
|
1 %% Audio Declipping Example
|
ivan@136
|
2 %
|
ivan@136
|
3 % CHANGE!!! This example is based on the experiment suggested by Professor Pierre
|
ivan@136
|
4 % Vandergheynst on the SMALL meeting in Villars.
|
ivan@136
|
5 % The idea behind is to use patches from source image as a dictionary in
|
ivan@136
|
6 % which we represent target image using matching pursuit algorithm.
|
ivan@136
|
7 % Calling Pierre_Problem function to get src image to be used as dictionary
|
ivan@136
|
8 % and target image to be represented using MP with 3 patches from source image
|
ivan@136
|
9
|
ivan@136
|
10 %
|
ivan@136
|
11 % Centre for Digital Music, Queen Mary, University of London.
|
ivan@136
|
12 % This file copyright 2011 Ivan Damnjanovic.
|
ivan@136
|
13 %
|
ivan@136
|
14 % This program is free software; you can redistribute it and/or
|
ivan@136
|
15 % modify it under the terms of the GNU General Public License as
|
ivan@136
|
16 % published by the Free Software Foundation; either version 2 of the
|
ivan@136
|
17 % License, or (at your option) any later version. See the file
|
ivan@136
|
18 % COPYING included with this distribution for more information.
|
ivan@136
|
19 %
|
ivan@136
|
20 %%
|
ivan@136
|
21
|
ivan@136
|
22 clear all;
|
ivan@136
|
23
|
ivan@136
|
24 % Defining the Problem structure
|
ivan@136
|
25
|
ivan@136
|
26 SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2);
|
ivan@136
|
27
|
ivan@136
|
28 % % Show original image and image that is used as a dictionary
|
ivan@136
|
29 % figure('Name', 'Original and Dictionary Image');
|
ivan@136
|
30 %
|
ivan@136
|
31 % subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval);
|
ivan@136
|
32 % title('Original Image');colormap(gray);axis off; axis image;
|
ivan@136
|
33 % subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval);
|
ivan@136
|
34 % title('Dictionary image:');colormap(gray);axis off; axis image;
|
ivan@136
|
35 time=0;
|
ivan@137
|
36 error2=0.001^2;
|
ivan@136
|
37 coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
|
ivan@136
|
38
|
ivan@136
|
39 for i=1:SMALL.Problem.n
|
ivan@136
|
40
|
ivan@136
|
41 idx = find(SMALL.Problem.M(:,i));
|
ivan@137
|
42 if size(idx,1)==SMALL.Problem.m
|
ivan@137
|
43 continue
|
ivan@137
|
44 end
|
ivan@137
|
45 Dict = SMALL.Problem.B(idx,:);
|
ivan@137
|
46 wDict = 1./sqrt(diag(Dict'*Dict));
|
ivan@137
|
47
|
ivan@137
|
48 SMALL.Problem.A = Dict*diag(wDict);
|
ivan@136
|
49
|
ivan@136
|
50 SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
|
ivan@136
|
51
|
ivan@137
|
52
|
ivan@136
|
53
|
ivan@136
|
54 % Defining the parameters sparse representation
|
ivan@136
|
55 SMALL.solver=SMALL_init_solver;
|
ivan@136
|
56 SMALL.solver.toolbox='ompbox';
|
ivan@137
|
57 SMALL.solver.name='omp2Gabor';
|
ivan@136
|
58
|
ivan@136
|
59 SMALL.solver.param=struct(...
|
ivan@137
|
60 'epsilon', error2*size(idx,1),...
|
ivan@137
|
61 'maxatoms', 128);
|
ivan@136
|
62
|
ivan@136
|
63 % Find solution
|
ivan@136
|
64
|
ivan@136
|
65 SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver);
|
ivan@136
|
66
|
ivan@136
|
67
|
ivan@137
|
68 coeffFrames(:,i) = wDict .* SMALL.solver.solution;
|
ivan@136
|
69 time = time + SMALL.solver.time;
|
ivan@136
|
70
|
ivan@136
|
71
|
ivan@136
|
72
|
ivan@136
|
73 end
|
ivan@136
|
74
|
ivan@136
|
75 %% Set reconstruction function
|
ivan@136
|
76
|
ivan@136
|
77 SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
|
ivan@136
|
78 reconstructed=SMALL.Problem.reconstruct(coeffFrames);
|
ivan@136
|
79
|
ivan@136
|
80
|
ivan@137
|
81
|
ivan@137
|
82 %% Plot results
|
ivan@137
|
83
|
ivan@137
|
84 xClipped = SMALL.Problem.clipped;
|
ivan@137
|
85 xClean = SMALL.Problem.original;
|
ivan@137
|
86 xEst1 = reconstructed.audioAllSamples;
|
ivan@137
|
87 xEst2 = reconstructed.audioOnlyClipped;
|
ivan@137
|
88 fs = SMALL.Problem.fs;
|
ivan@137
|
89
|
ivan@137
|
90 figure
|
ivan@137
|
91 hold on
|
ivan@137
|
92 plot(xClipped,'r')
|
ivan@137
|
93 plot(xClean)
|
ivan@137
|
94 plot(xEst2,'--g')
|
ivan@137
|
95 plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
|
ivan@137
|
96 legend('Clipped','True solution','Estimate')
|
ivan@137
|
97
|
ivan@137
|
98 % Normalized and save sounds
|
ivan@137
|
99 normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
|
ivan@137
|
100 L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
|
ivan@137
|
101 xEst1 = xEst1(1:L)/normX;
|
ivan@137
|
102 xEst2 = xEst2(1:L)/normX;
|
ivan@137
|
103 xClipped = xClipped(1:L)/normX;
|
ivan@137
|
104 xClean = xClean(1:L)/normX;
|
ivan@137
|
105 wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']);
|
ivan@137
|
106 wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']);
|
ivan@137
|
107 wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);
|
ivan@137
|
108 wavwrite(xClean,fs,[expParam.destDir 'xClean.wav']);
|