Mercurial > hg > smallbox

comparison toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_consOMP.m @ 144:19e0af570914 release_1.5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge from branch "ivand_dev"
author Ivan <ivan.damnjanovic@eecs.qmul.ac.uk>
date Tue, 26 Jul 2011 15:14:15 +0100
parents 56d719a5fd31
children
comparison
equal deleted inserted replaced
143:8d866d96f006 144:19e0af570914
1 function y = inpaintFrame_consOMP(problemData,param)
2 % Inpainting method based on OMP with a constraint
3 % on the amplitude of the reconstructed samples an optional constraint
4 % on the maximum value of the clipped samples
5 %
6 % Usage: y = inpaintFrame_consOMP(problemData,param)
7 %
8 %
9 % Inputs:
10 % - problemData.x: observed signal to be inpainted
11 % - problemData.Imiss: Indices of clean samples
12 % - param.D - the dictionary matrix (optional if param.D_fun is set)
13 % - param.D_fun - a function handle that generates the dictionary
14 % matrix param.D if param.D is not given. See, e.g., DCT_Dictionary.m and Gabor_Dictionary.m
15 % - param.wa - Analysis window
16 % - param.Upper_Limit - if present and non-empty this fiels
17 % indicates that an upper limit constraint is active and its
18 % integer value is such that
19 %
20 % Outputs:
21 % - y: estimated frame
22 %
23 % Note that the CVX library is needed.
24 %
25 % -------------------
26 %
27 % Audio Inpainting toolbox
28 % Date: June 28, 2011
29 % By Valentin Emiya, Amir Adler, Michael Elad, Maria Jafari
30 % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
31 % ========================================================
32
33 %% Load data and parameters
34
35 x = problemData.x;
36 IObs = find(~problemData.IMiss);
37 p.N = length(x);
38 E2 = param.OMPerr^2;
39 E2M=E2*length(IObs);
40 wa = param.wa(param.N);
41
42 % build the dictionary matrix if only the dictionary generation function is given
43 if ~isfield(param,'D')
44 param.D = param.D_fun(param);
45 end
46
47
48 % clipping level detection
49 clippingLevelEst = max(abs(x(:)./wa(:)));
50
51 IMiss = true(length(x),1);
52 IMiss(IObs) = false;
53 IMissPos = find(x>=0 & IMiss);
54 IMissNeg = find(x<0 & IMiss);
55
56 DictPos=param.D(IMissPos,:);
57 DictNeg=param.D(IMissNeg,:);
58
59 % Clipping level: take the analysis window into account
60 wa_pos = wa(IMissPos);
61 wa_neg = wa(IMissNeg);
62 b_ineq_pos = wa_pos(:)*clippingLevelEst;
63 b_ineq_neg = -wa_neg(:)*clippingLevelEst;
64 if isfield(param,'Upper_Limit') && ~isempty(param.Upper_Limit)
65 b_ineq_pos_upper_limit = wa_pos(:)*param.Upper_Limit*clippingLevelEst;
66 b_ineq_neg_upper_limit = -wa_neg(:)*param.Upper_Limit*clippingLevelEst;
67 else
68 b_ineq_pos_upper_limit = Inf;
69 b_ineq_neg_upper_limit = -Inf;
70 end
71
72 %%
73 Dict=param.D(IObs,:);
74 W=1./sqrt(diag(Dict'*Dict));
75 Dict=Dict*diag(W);
76 xObs=x(IObs);
77
78 residual=xObs;
79 maxNumCoef = param.sparsityDegree;
80 indx = [];
81 currResNorm2 = E2M*2; % set a value above the threshold in order to have/force at least one loop executed
82 j = 0;
83 while currResNorm2>E2M && j < maxNumCoef,
84 j = j+1;
85 proj=Dict'*residual;
86 [dum pos] = max(abs(proj));
87 indx(j)=pos;
88 a=pinv(Dict(:,indx(1:j)))*xObs;
89 residual=xObs-Dict(:,indx(1:j))*a;
90 currResNorm2=sum(residual.^2);
91 end;
92
93
94 if isinf(b_ineq_pos_upper_limit)
95 %% CVX code
96 cvx_begin
97 cvx_quiet(true)
98 variable a(j)
99 %minimize( sum(square(xObs-Dict*a)))
100 minimize(norm(Dict(:,indx)*a-xObs))
101 subject to
102 DictPos(:,indx)*(W(indx).*a) >= b_ineq_pos
103 DictNeg(:,indx)*(W(indx).*a) <= b_ineq_neg
104 cvx_end
105 if cvx_optval>1e3
106 cvx_begin
107 cvx_quiet(true)
108 variable a(j)
109 minimize(norm(Dict(:,indx)*a-xObs))
110 cvx_end
111 end
112 else
113 %% CVX code
114 cvx_begin
115 cvx_quiet(true)
116 variable a(j)
117 %minimize( sum(square(xObs-Dict*a)))
118 minimize(norm(Dict(:,indx)*a-xObs))
119 subject to
120 DictPos(:,indx)*(W(indx).*a) >= b_ineq_pos
121 DictNeg(:,indx)*(W(indx).*a) <= b_ineq_neg
122 DictPos(:,indx)*(W(indx).*a) <= b_ineq_pos_upper_limit
123 DictNeg(:,indx)*(W(indx).*a) >= b_ineq_neg_upper_limit
124 cvx_end
125 if cvx_optval>1e3
126 cvx_begin
127 cvx_quiet(true)
128 variable a(j)
129 minimize(norm(Dict(:,indx)*a-xObs))
130 cvx_end
131 end
132 end
133
134 %% Frame Reconstruction
135 indx(length(a)+1:end) = [];
136
137 Coeff = sparse(size(param.D,2),1);
138 if (~isempty(indx))
139 Coeff(indx) = a;
140 Coeff = W.*Coeff;
141 end
142 y = param.D*Coeff;
143
144 return