function inpaint = SinModInt(inpaint)

%% Function for interpolation of missing samples
% input parameters:
% inpaint.signal.samples - input signal
% inpaint.signal.fs - sampling frequency of input signal
% inpaint.miss - parameters of missing region
% inpaint.solver.param.KB - order of AR model of the signal left from gap 
% inpaint.solver.param.KBKA - order of AR model of the signal right from gap 
% inpaint.solver.param.KBKNB - order of AR model of the residual signal left from gap 
% inpaint.solver.param.KBKNA - order of AR model of the residual signal right from gap 
% inpaint.solver.param.KBthrF - threshold for frequency vector distance
% inpaint.solver.param.KBthrA - threshold for amplitude vector distance
% inpaint.solver.param.KBparM - maximum length of amplitude vector for amplitude estimate
% inpaint.solver.param.KBlB - length of extrapolated vectors from the left side
% inpaint.solver.param.KBlA - length of extrapolated vectors from the right side
% inpaint.solver.param.KBsolver: Long Interpolation based on Sinusoidal Modeling
% output parameter:
% inpaint.sol.err - signal with missing region
% inpaint.sol.dr1 - interpolated signal without residual noise
% inpaint.sol.samples - interpolated signal
%%
% clear all
% close all
addpath('./sinemodel')
addpath('./funkce')

if nargin < 1
    error('Nedostatek vstupních parametrů')
end

if (inpaint.miss(1,1)+inpaint.miss(1,2))>length(inpaint.signal.samples)
    error('Velikost mezery přesahuje délku signálu')
end

param = inpaint.solver.param;

if ~isstruct(param)
    error('Parametry nemají strukturu pole')
end

if ~isfield(param,'thrF')
    param.thrF = 2;  %default
end

if ~isfield(param,'thrA')
    param.thrA = 0.5;  %default
end

if ~isfield(param,'parM')
    param.parM = 30;  %default
end



%vektor vzorků
d = inpaint.signal.samples;
%vzorkovací kmitočet
sr = inpaint.signal.fs;

inpaint.sol.orig=d;


%==========================================================================
%==========================================================================

%Mezera
Me = inpaint.miss(1,2);  %délka
l = inpaint.miss(1,1); %počátek mezery

d(l:l+Me-1)=0; %mezera

%výstup signálu s mezerou
inpaint.sol.err = d;


%sinusoidální model, pro potřeby výpočtu délky mezery
[F,M,P]=sinmodel(d,sr);


%==========================================================================
%==========================================================================
%Levá strana

dL = d(1:l-1);

[FB,MB,PB]=sinmodel(dL,sr);

%unwrapp fáze
PB = mod(PB,2*pi);

%přidá 0 na okraje matic kvůli ztrátám okrajů při syntetizaci
FB0 = [0*FB(:,1),FB,0*FB(:,end)];
MB0 = [0*MB(:,1),MB,0*MB(:,end)];
PB0 = [0*PB(:,1),PB,0*PB(:,end)];

%výpočet residua
dLsynth = synthphtrax(FB0,MB0,PB0,sr,256,128);
dLres = (dL - dLsynth(1:length(dL))')';

%==========================================================================
%Pravá strana

dR = d(l+Me:length(d));
inpaint.sol.right = dR;


[FA,MA,PA]=sinmodel(dR,sr);

%unwrapp fáze
PA = mod(PA,2*pi);

%přidá 0 na okraje matic kvůli ztrátám okrajů při syntetizaci
FA0 = [0*FA(:,1),FA,0*FA(:,end)];
MA0 = [0*MA(:,1),MA,0*MA(:,end)];
PA0 = [0*PA(:,1),PA,0*PA(:,end)];

%výpočet residua
dRsynth = synthphtrax(FA0,MA0,PA0,sr,256,128);
dRres = (dR - dRsynth(1:length(dR))')';

%=========================================================================
%INTERPOLACE

%Délka mezery
FPlen=size(F,2)-size(FA,2)-size(FB,2);
%==========================================================================

%LEVÁ STRANA

%alokace paměti
FBI=zeros(size(FB,1),FPlen);
MBI=zeros(size(MB,1),FPlen);

for i=1:size(FB,1);
    
    if isnan(FB(i,end)) %přeskočí řádek pokud na kraji nejsou hodnoty
        continue;
    end          
    %Zjištění indexu, na kterém začínají hodnoty složky
    iB=0; %index matice na kterém začínají hodnoty (frekvence, amplitudy)
    for columnB=size(FB,2):-1:1;
        if isnan(FB(i,columnB))
            iB=columnB;
            break;
        end
    end
    
%     if ~isfield(param,'KB')
%     param.KB = round(min(40,size(FB,2)-iB)/2);  %default
%     end

    %pokud je na kraji pouze jedna hodnota, nakopíruje ji,
    %jinak provede predikci
    if length(FB(i,iB+1:end))~=1

    KB=round(min(40,size(FB,2)-iB)/2); %řád AR modelu
   
    %ODHAD FREKVENCE
    
    aFB=BurgMeth(FB(i,iB+1:end),KB); %AR koeficienty
    %predikce hodnot frekvence
    FBI(i,:) = interleft(FB(i,iB+1:end),aFB,KB,FPlen);  
    
    %ODHAD AMPLITUDY
    
    aMB=BurgMeth(MB(i,iB+1:end),KB);  %AR koeficienty

    MBI(i,:) = interleft(MB(i,iB+1:end),aMB,KB,FPlen);

    else %nakopíruje hodnotu
       FBI(i,:) = FB(i,end);
       MBI(i,:) = MB(i,end);      
    end   
end

%ODHAD RESIDUA
 if ~isfield(param,'KNB') %nastavení parametru řád AR modelu
    param.KNB = round(length(dLres)/2);  %default
 end
 
KNB = param.KNB; %Řád AR levé strany
aNB = BurgMeth(dLres,KNB);
NBI = interleft(dLres,aNB,KNB,Me);


%==========================================================================
%PRAVÁ STRANA
%alokace paměti
FAI=zeros(size(FA,1),FPlen);
MAI=zeros(size(MA,1),FPlen);

for i=1:size(FA,1);    
    if isnan(FA(i,1)) %přeskočí řádek pokud na kraji nejsou hodnoty
            continue;
    end
    %Zjištění indexu, na kterém začínají hodnoty složky
    iA=0; %index matice na kterém končí hodnoty (frekvence, amplitudy)
    for columnA=1:size(FA,2);
        if (isnan(FA(i,columnA)))
            break;
        else
            iA=columnA;            
        end
    end
    
%     if ~isfield(param,'KA')
%     param.KA = round(min(40,iA)/2);  %default
%     end
    
    %pokud je na kraji pouze jedna hodnota, nakopíruje ji,
    %jinak provede predikci
    if length(FA(i,1:iA))~=1

    KA=round(min(40,iA)/2); %řád AR modelu

    %ODHAD FREKVENCE
    
    aFA=BurgMeth(fliplr(FA(i,1:iA)),KA);    %AR koeficienty
    %predikce hodnot frekvence
    FAI(i,:) = fliplr(interleft(fliplr(FA(i,1:iA)),aFA,KA,FPlen));
  
    %ODHAD AMPLITUDY

    aMA=BurgMeth(fliplr(MA(i,1:iA)),KA);%AR koeficienty
    %predikce hodnot amplitudy
    MAI(i,:) = fliplr(interleft(fliplr(MA(i,1:iA)),aMA,KA,FPlen));
    
    else %nakopíruje hodnotu    
       FAI(i,:) = FA(i,1);
       MAI(i,:) = MA(i,1);       
    end
end

%ODHAD RESIDUA
if ~isfield(param,'KNA')%nastavení parametru řád AR modelu
    param.KNA = round(length(dRres)/2);  %default
end

KNA=param.KNA; %Řád AR levé strany

aNA=BurgMeth(fliplr(dRres),KNA);
NAI = fliplr(interleft(fliplr(dRres),aNA,KNA,Me));




%Najdi shodu
%==========================================================================
%==========================================================================
%INTERPOLACE


thrF = param.thrF; %threshold F 
thrA = param.thrA; %threshold A           
parM = param.parM; %parametr M pro délu vektoru hodnot amplitudy 
%váhovací okno
w=okno(FPlen);
% pocitadlo radku
counterRow = 1;
% vector indexu "i" a "j" (použitých řádků)
indexI=[];
indexJ=[];

FI=[];
MI=[];
PI=[];

% spojení odhadovaných hodnot z matice FBI a FAI
sizeI=size(FBI,1);
sizeJ=size(FAI,1);
for i=1:sizeI
    
    %přeskočí tádek pokud je již zpracován nebo obsahuje nuly
    if  (isempty(find((indexI - i) == 0)) == 0) || (FBI(i,1)==0)
        continue;
    end
    
    for j=1:sizeJ
        
        %přeskočí tádek pokud je již zpracován nebo obsahuje nuly
        if  (isempty(find((indexJ - j) == 0)) == 0) || (FAI(j,1)==0)
        continue;
        end
        
        %Podmínky pro spojení protějších vektorů
        if (matchvect(FBI(i,:),FAI(j,:))) < thrF && (matchvect(MBI(i,:),...
            MAI(j,:))) < thrA && (matchvect(FBI(i,:),FAI(j,:))) ~=0
        
            %váhový součet vektorů frekvence
            FP=(w.*FBI(i,:))+((1.-w).*FAI(j,:));
            
            iB=0; %index matice na kterém začínají hodnoty (frekvence, amplitudy)
            for columnB=size(MB,2):-1:1;
                 if isnan(MB(i,columnB))
                     iB=columnB;
                 break;
                 end
            end
            
            iA=0; %konec vektoru hodnot frekvence
            for columnA=1:size(MA,2);
                if (isnan(MA(j,columnA)))
                break;
                else
                    iA=columnA;            
                end
            end

            %%váhový součet vektorů amplitudy
            deltaB = vypDeli(MBI(i,:),MA(j,1:iA),parM,FPlen);
            deltaA = vypDelj(MAI(j,:),MB(i,iB+1:end),parM,FPlen);
            MP=(w.*(MBI(i,:)+deltaB))+((1.-w).*(MAI(j,:)+deltaA));

            %výpočet vektoru fáze
            PBI = interphase(PB(i,iB+1:end),PA(j,1:iA),FB(i,iB+1:end),FBI(i,:),FA(j,1),FPlen,sr);

            %zápis do matice odhadovaných parametrů
            FI(counterRow,:) = [0*FP(:,1),FP,0*FP(:,end)];
            MI(counterRow,:) = [0*MP(:,1),MP,0*MP(:,end)];                                
            PI(counterRow,:) = [0*PBI(:,1),PBI,0*PBI(:,end)];
            
            counterRow = counterRow + 1;
            
            %Zapíše indexy řádků, které už jsou zpracovány
            if isempty(indexI)
                indexI(1)=i;
            else
                indexI(end+1)=i;
            end;
            if isempty(indexJ)
                indexJ(1)=j;
            else
                indexJ(end+1)=j;
            end
            
        end %if        
    end %for j
end %for i
%==========================================================================
%EXTRAPOLACE

%LEVÁ STRANA
if ~isfield(param,'lB')
    param.lB = size(MBI,2)-1;  %default
 end
lB = param.lB; %délka extrapolovaného vektoru

for i=1:size(FB,1)
    
    %přeskočí tádek pokud je již zpracován
    if  (isempty(find((indexI - i) == 0)) == 0)
        continue;
    end
    
    %vybere řádek na jehož kraji jsou hodnoty
    if isnan(FB(i,end))
    else 
        
         %extrapolace amplitudy, indexS je konec vektoru
         [MP,indexS]=extraampB(MBI(i,:),lB);          
         MP = vecnanB(MP,indexS); %doplní NaN 
         
         %Pokud je celý vektor NaN, pokračuje na další řádek
         if isnan(MP(1))
             continue;
         end
         
         FP=vecnanB(FBI(i,:),indexS); %doplní NaN 
         %extrapolace fáze
         PBI = extraphase(PB(i,:),FB(i,:),FBI(i,:),FPlen,sr);
         PBI = vecnanB(PBI,indexS); %doplní NaN 
         
         %zápis do matice odhadovaných parametrů      
         FI(counterRow,:) = [0*FP(:,1),FP,0*FP(:,end)];
         MI(counterRow,:) = [0*MP(:,1),MP,0*MP(:,end)];                                
         PI(counterRow,:) = [0*PBI(:,1),PBI,0*PBI(:,end)];
         counterRow = counterRow + 1;
         
         %zapíše zpracovaný řádek
          if isempty(indexI)
                indexI(1)=i;
          else
                indexI(end+1)=i;
          end;
    end    
end

%PRAVÁ STRANA

if ~isfield(param,'lA')
    param.lA = round((size(MAI,2)-1)/2);  %default
end
lA=param.lA;
 

for j=1:size(FA,1)
    if  (isempty(find((indexJ - j) == 0)) == 0)
        continue;
    end
    if isnan(FA(j,1))
    else
                    
         [MP, indexS]=extraampA(MAI(j,:),lA);
         MP=vecnanA(MP,indexS);
         
         if isnan(MP(end))
             continue;
         end
         
         FP=vecnanA(FAI(j,:),indexS);        
         PAI = fliplr(extraphase(fliplr(PA(j,:)),fliplr(FA(j,:)),fliplr(FAI(j,:)),FPlen,sr));
         PAI = vecnanA(PAI,indexS);

         FI(counterRow,:) = [0*FP(:,1),FP,0*FP(:,end)];
         MI(counterRow,:) = [0*MP(:,1),MP,0*MP(:,end)];                                
         PI(counterRow,:) = [0*PAI(:,1),PAI,0*PAI(:,end)];
         counterRow = counterRow + 1;
         
         % jiz zakazany index
         if isempty(indexJ)
                indexJ(1)=j;
         else
                indexJ(end+1)=j;
         end
    end    
end
%==========================================================================
%váhový součet residua
w=okno(Me);
NP = (w.*NBI)+((1.-w).*NAI);
%==========================================================================

%syntéza nepoškozeného signálu
[F,M,P]=sinmodel(inpaint.sol.orig,sr);

F0 = [0*F(:,1),F,0*F(:,end)];
M0 = [0*M(:,1),M,0*M(:,end)];
P0 = [0*P(:,1),P,0*P(:,end)];

inpaint.sol.drorig = synthphtrax(F0,M0,P0,sr,256,128,length(d)/sr);


%výstup rekonstruovaného signálu
if ~isempty(FI)
dr1miss = synthphtrax(FI,MI,PI,sr,256,128);
d(l:l+Me-1)=dr1miss(128:Me+127);
end
inpaint.sol.dr1=d;
d(l:l+Me-1)=d(l:l+Me-1)+NP';

inpaint.sol.samples = d;






end














