annotate toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_janssenInterpolation.m @ 147:65fc57f3903c ivand_dev

Merge from the default branch
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Tue, 26 Jul 2011 15:55:14 +0100
parents 56d719a5fd31
children
rev   line source
ivan@138 1 function y = inpaintFrame_janssenInterpolation(problemData,param)
ivan@138 2 % Frame-level inpainting method based on the linear prediction by
ivan@138 3 % Janssen.
ivan@138 4 %
ivan@138 5 % Usage: xEst = inpaintFrame_janssenInterpolation(problemData,param)
ivan@138 6 %
ivan@138 7 %
ivan@138 8 % Inputs:
ivan@138 9 % - problemData.x - observed signal to be inpainted
ivan@138 10 % - problemData.Imiss - Indices of clean samples
ivan@138 11 % - param.p - Order of the autoregressive model used in
ivan@138 12 % for linear prediction
ivan@138 13 %
ivan@138 14 % Outputs:
ivan@138 15 % - y: estimated frame
ivan@138 16 %
ivan@138 17 % References:
ivan@138 18 %
ivan@138 19 % -------------------
ivan@138 20 %
ivan@138 21 % Audio Inpainting toolbox
ivan@138 22 % Date: June 28, 2011
ivan@138 23 % By Valentin Emiya, Amir Adler, Maria Jafari
ivan@138 24 % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
ivan@138 25
ivan@138 26 s = problemData.x;
ivan@138 27 N = length(s);
ivan@138 28 Im = find(problemData.IMiss);
ivan@138 29 IObs = find(~problemData.IMiss);
ivan@138 30 M = length(Im);
ivan@138 31 Im = sort(Im); Im = Im(:); % Im: indexes of missing samples
ivan@138 32 s(Im) = 0;
ivan@138 33
ivan@138 34 if nargin<2 || ~isfield(param,'p')
ivan@138 35 p = min(3*M+2,round(N/3));
ivan@138 36 else
ivan@138 37 p = param.p;
ivan@138 38 end
ivan@138 39 if nargin<2 || ~isfield(param,'GR')
ivan@138 40 param.GR = false;
ivan@138 41 end
ivan@138 42 if nargin<2 || ~isfield(param,'NIt')
ivan@138 43 NIt = 100;
ivan@138 44 else
ivan@138 45 NIt = param.NIt;
ivan@138 46 end
ivan@138 47
ivan@138 48
ivan@138 49 IAA = abs(Im*ones(1,N)-ones(M,1)*(1:N));
ivan@138 50 IAA1 = IAA<=p;
ivan@138 51 AA = zeros(size(IAA));
ivan@138 52
ivan@138 53 if param.GR
ivan@138 54 figure;
ivan@138 55 hold on
ivan@138 56 end
ivan@138 57 for k=1:NIt
ivan@138 58
ivan@138 59 % Re-estimation of LPC
ivan@138 60 aEst = lpc(s,p).';
ivan@138 61
ivan@138 62 % Re-estimation of the missing samples
ivan@138 63 b = aEst.'*hankel(aEst.',[aEst(end),zeros(1,p)]);
ivan@138 64 AA(IAA1) = b(IAA(IAA1)+1);
ivan@138 65 % xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs); % use Chol to invert matrix
ivan@138 66 [R flagErr] = chol(AA(:,Im));
ivan@138 67 if flagErr
ivan@138 68 % xEst = - AA(:,Im)\(AA(:,IObs)*s(IObs));
ivan@138 69 xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs);
ivan@138 70 else
ivan@138 71 xEst = -R\(R'\(AA(:,IObs)*s(IObs)));
ivan@138 72 end
ivan@138 73 s(Im) = xEst;
ivan@138 74 if param.GR
ivan@138 75 e = filter(aEst,1,s);
ivan@138 76 plot(k,10*log10(mean(e(p+1:end).^2)),'o')
ivan@138 77 end
ivan@138 78 end
ivan@138 79
ivan@138 80 y = s;
ivan@138 81
ivan@138 82 return