diff toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_janssenInterpolation.m @ 138:56d719a5fd31 ivand_dev

Audio Inpaintin Toolbox

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_janssenInterpolation.m	Thu Jul 21 14:27:47 2011 +0100
@@ -0,0 +1,82 @@
+function y = inpaintFrame_janssenInterpolation(problemData,param)
+% Frame-level inpainting method based on the linear prediction by
+% Janssen.
+%
+% Usage: xEst = inpaintFrame_janssenInterpolation(problemData,param)
+%
+%
+% Inputs:
+%          - problemData.x - observed signal to be inpainted
+%          - problemData.Imiss - Indices of clean samples
+%          - param.p - Order of the autoregressive model used in
+%          for linear prediction
+%
+% Outputs:
+%          - y: estimated frame
+%
+% References:
+% 
+% -------------------
+%
+% Audio Inpainting toolbox
+% Date: June 28, 2011
+% By Valentin Emiya, Amir Adler, Maria Jafari
+% This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
+
+s = problemData.x;
+N = length(s);
+Im = find(problemData.IMiss);
+IObs = find(~problemData.IMiss);
+M = length(Im);
+Im = sort(Im); Im = Im(:); % Im: indexes of missing samples
+s(Im) = 0;
+
+if nargin<2 || ~isfield(param,'p')
+   p = min(3*M+2,round(N/3));
+else
+   p = param.p;
+end
+if nargin<2 || ~isfield(param,'GR')
+   param.GR = false;
+end
+if nargin<2 || ~isfield(param,'NIt')
+   NIt = 100;
+else
+   NIt = param.NIt;
+end
+
+
+IAA = abs(Im*ones(1,N)-ones(M,1)*(1:N));
+IAA1 = IAA<=p;
+AA = zeros(size(IAA));
+
+if param.GR
+   figure;
+   hold on
+end
+for k=1:NIt
+   
+   % Re-estimation of LPC
+   aEst = lpc(s,p).';
+   
+   % Re-estimation of the missing samples
+   b = aEst.'*hankel(aEst.',[aEst(end),zeros(1,p)]);
+   AA(IAA1) = b(IAA(IAA1)+1);
+%    xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs); % use Chol to invert matrix
+   [R flagErr] = chol(AA(:,Im));
+   if flagErr
+      %       xEst = - AA(:,Im)\(AA(:,IObs)*s(IObs));
+      xEst = -inv(AA(:,Im))*AA(:,IObs)*s(IObs);
+   else
+      xEst = -R\(R'\(AA(:,IObs)*s(IObs)));
+   end
+   s(Im) = xEst;
+   if param.GR
+      e = filter(aEst,1,s);
+      plot(k,10*log10(mean(e(p+1:end).^2)),'o')
+   end
+end
+
+y = s;
+
+return