ivan@138: function y = inpaintFrame_OMP(problemData,param)
ivan@138: % Inpainting method based on OMP 
ivan@138: %
ivan@138: % Usage: y = inpaintFrame_OMP(problemData,param)
ivan@138: %
ivan@138: %
ivan@138: % Inputs:
ivan@138: %          - problemData.x: observed signal to be inpainted
ivan@138: %          - problemData.Imiss: Indices of clean samples
ivan@138: %          - param.D - the dictionary matrix (optional if param.D_fun is set)
ivan@138: %          - param.D_fun - a function handle that generates the dictionary 
ivan@138: %          matrix param.D if param.D is not given. See, e.g., DCT_Dictionary.m and Gabor_Dictionary.m
ivan@138: %          - param.wa - Analysis window
ivan@138: %
ivan@138: % Outputs:
ivan@138: %          - y: estimated frame
ivan@138: %
ivan@138: %
ivan@138: % -------------------
ivan@138: %
ivan@138: % Audio Inpainting toolbox
ivan@138: % Date: June 28, 2011
ivan@138: % By Valentin Emiya, Amir Adler, Michael Elad, 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: % To do next: use a faster implementation of OMP
ivan@138: 
ivan@138: %% Load data and parameters
ivan@138: 
ivan@138: x = problemData.x;
ivan@138: IObs = find(~problemData.IMiss);
ivan@138: p.N = length(x);
ivan@138:    E2 = param.OMPerr^2;
ivan@138:    E2M=E2*length(IObs);
ivan@138: 
ivan@138:    wa = param.wa(param.N);
ivan@138: 
ivan@138: 
ivan@138: %% Build and normalized dictionary
ivan@138: % build the dictionary matrix if only the dictionary generation function is given
ivan@138: if ~isfield(param,'D')
ivan@138:    param.D = param.D_fun(param);
ivan@138: end
ivan@138: Dict=param.D(IObs,:);
ivan@138: W=1./sqrt(diag(Dict'*Dict));
ivan@138: Dict=Dict*diag(W);
ivan@138: xObs=x(IObs);
ivan@138: 
ivan@138: %% OMP iterations
ivan@138: residual=xObs;
ivan@138: maxNumCoef = param.sparsityDegree;
ivan@138: indx = [];
ivan@138: currResNorm2 = E2M*2; % set a value above the threshold in order to have/force at least one loop executed
ivan@138: j = 0;
ivan@138: while currResNorm2>E2M && j < maxNumCoef,
ivan@138:    j = j+1;
ivan@138:    proj=Dict'*residual;
ivan@138:    [dum pos] = max(abs(proj));
ivan@138:    indx(j)=pos;
ivan@138:    
ivan@138:    a=pinv(Dict(:,indx(1:j)))*xObs;
ivan@138:    
ivan@138:    residual=xObs-Dict(:,indx(1:j))*a;
ivan@138:    currResNorm2=sum(residual.^2);
ivan@138: end
ivan@138: 
ivan@138: %% Frame Reconstruction
ivan@138: indx(length(a)+1:end) = [];
ivan@138: 
ivan@138: Coeff = sparse(size(param.D,2),1);
ivan@138: if (~isempty(indx))
ivan@138:    Coeff(indx) = a;
ivan@138:    Coeff = W.*Coeff;
ivan@138: end
ivan@138: y = param.D*Coeff;
ivan@138: 
ivan@138: return