Home > NoiseTools > fillgap.m

# fillgap

## PURPOSE

y=fillgap(x,w,maxGapSize) - fill gaps using simple interpolation

## SYNOPSIS

function [y,w]=fillgap(x,w,maxGapSize)

## DESCRIPTION

    y=fillgap(x,w,maxGapSize) - fill gaps using simple interpolation

      y: interpolated data

      x: data to interpolate
      w: weighting function
      maxGapSize: largest expected gap size

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- nt_inpaint function y=nt_inpaint(x,w) - weighted interpolation based on correlation structure
- nt_outliers [w,y]=nt_outliers(x,w,thresh) - detect outliers based on weighted correlation structure

## SOURCE CODE

    0001 function [y,w]=fillgap(x,w,maxGapSize)
    0002 %y=fillgap(x,w,maxGapSize) - fill gaps using simple interpolation
    0003 %
    0004 %  y: interpolated data
    0005 %
    0006 %  x: data to interpolate
    0007 %  w: weighting function
    0008 %  maxGapSize: largest expected gap size
    0009 if nargin<2; error('!'); end
    0010 if nargin<3||isempty(maxGapSize); maxGapSize=1; end
    0011 if size(x,2)>1; error('!'); end
    0012 if size(x) ~= size(w); error('!'); end
    0013 y=x;
    0014 if all(w); return; end
    0015 % simple case size one
    0016 iToFix=1+find(~w(2:end-1)&w(1:end-2)&w(3:end));
    0017 y(iToFix)=(y(iToFix-1)+y(iToFix+1))/2;
    0018 w(iToFix)=1;
    0019 % general case size > 1
    0020 iStart=find(w(1:end-2) & ~w(2:end-1));  % position preceding gap
    0021 iStop=find(~w(1:end-1) & w(2:end));     % last position in gap
    0022 if isempty(iStart)||isempty(iStop); return; end
    0023 if iStop(1)<iStart(1);
    0024     iStop=iStop(2:end);                 % ignore gap at beginning
    0025 end
    0026 iStart=iStart(1:numel(iStop));          % ignore gap at end
    0027 for gapSize=2:maxGapSize
    0028     idx=find(iStop-iStart==gapSize);
    0029     for k=1:gapSize
    0030         % interpolate between samples on either side of gap
    0031         y(iStart(idx)+k) = ( y(iStart(idx)) * (gapSize-k+1) + y(iStart(idx)+gapSize+1) * k ) / (gapSize+1);
    0032         w(iStart(idx)+k) = 1;
    0033     end
    0034 end