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
|