Mercurial > hg > smallbox

comparison toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_consOMP_Gabor.m @ 138:56d719a5fd31 ivand_dev

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