ivan@138: function y = inpaintFrame_janssenInterpolation(problemData,param)
ivan@138: % Frame-level inpainting method based on the linear prediction by
ivan@138: % Janssen.
ivan@138: %
ivan@138: % Usage: xEst = inpaintFrame_janssenInterpolation(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.p - Order of the autoregressive model used in
ivan@138: %          for linear prediction
ivan@138: %
ivan@138: % Outputs:
ivan@138: %          - y: estimated frame
ivan@138: %
ivan@138: % References:
ivan@138: % 
ivan@138: % -------------------
ivan@138: %
ivan@138: % Audio Inpainting toolbox
ivan@138: % Date: June 28, 2011
ivan@138: % By Valentin Emiya, Amir Adler, 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: s = problemData.x;
ivan@138: N = length(s);
ivan@138: Im = find(problemData.IMiss);
ivan@138: IObs = find(~problemData.IMiss);
ivan@138: M = length(Im);
ivan@138: Im = sort(Im); Im = Im(:); % Im: indexes of missing samples
ivan@138: s(Im) = 0;
ivan@138: 
ivan@138: if nargin<2 || ~isfield(param,'p')
ivan@138:    p = min(3*M+2,round(N/3));
ivan@138: else
ivan@138:    p = param.p;
ivan@138: end
ivan@138: if nargin<2 || ~isfield(param,'GR')
ivan@138:    param.GR = false;
ivan@138: end
ivan@138: if nargin<2 || ~isfield(param,'NIt')
ivan@138:    NIt = 100;
ivan@138: else
ivan@138:    NIt = param.NIt;
ivan@138: end
ivan@138: 
ivan@138: 
ivan@138: IAA = abs(Im*ones(1,N)-ones(M,1)*(1:N));
ivan@138: IAA1 = IAA<=p;
ivan@138: AA = zeros(size(IAA));
ivan@138: 
ivan@138: if param.GR
ivan@138:    figure;
ivan@138:    hold on
ivan@138: end
ivan@138: for k=1:NIt
ivan@138:    
ivan@138:    % Re-estimation of LPC
ivan@138:    aEst = lpc(s,p).';
ivan@138:    
ivan@138:    % Re-estimation of the missing samples
ivan@138:    b = aEst.'*hankel(aEst.',[aEst(end),zeros(1,p)]);
ivan@138:    AA(IAA1) = b(IAA(IAA1)+1);
ivan@138: %    xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs); % use Chol to invert matrix
ivan@138:    [R flagErr] = chol(AA(:,Im));
ivan@138:    if flagErr
ivan@138:       %       xEst = - AA(:,Im)\(AA(:,IObs)*s(IObs));
ivan@138:       xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs);
ivan@138:    else
ivan@138:       xEst = -R\(R'\(AA(:,IObs)*s(IObs)));
ivan@138:    end
ivan@138:    s(Im) = xEst;
ivan@138:    if param.GR
ivan@138:       e = filter(aEst,1,s);
ivan@138:       plot(k,10*log10(mean(e(p+1:end).^2)),'o')
ivan@138:    end
ivan@138: end
ivan@138: 
ivan@138: y = s;
ivan@138: 
ivan@138: return