function [dataout] = ft_interpolatenan(cfg, datain)

% FT_INTERPOLATENAN interpolates time series that contains segments of nans obtained

% by replacing artifactual data with nans using, for example, FT_REJECTARTIFACT, or

% by redefining trials with FT_REDEFINETRIAL resulting in trials with gaps.

%

% Use as

% outdata = ft_interpolatenan(cfg, indata)

% where cfg is a configuration structure and the input data is obtained from FT_PREPROCESSING.

%

% The configuration should contain

% cfg.method = string, interpolation method, see INTERP1 (default = 'linear')

% cfg.prewindow = value, length of data prior to interpolation window, in seconds (default = 1)

% cfg.postwindow = value, length of data after interpolation window, in seconds (default = 1)

% cfg.feedback = string, 'no', 'text', 'textbar', 'gui' (default = 'text')

%

% This function only interpolates over time, not over space. If you want to

% interpolate using spatial information, e.g. using neighbouring channels, you should

% use FT_CHANNELREPAIR.

%

% To facilitate data-handling and distributed computing, you can use

% cfg.inputfile = ...

% cfg.outputfile = ...

% If you specify one of these (or both) the input data will be read from a *.mat

% file on disk and/or the output data will be written to a *.mat file. These mat

% files should contain only a single variable, corresponding with the

% input/output structure.

%

% See also FT_REJECTARTIFACT, FT_REDEFINETRIAL, FT_CHANNELREPAIR

% Copyright (C) 2003-2020, Jan-Mathijs Schoffelen & Robert Oostenveld

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

%

% FieldTrip is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

% GNU General Public License for more details.

%

% You should have received a copy of the GNU General Public License

% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the initial part deals with parsing the input options and data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% these are used by the ft_preamble/ft_postamble function and scripts

ft_revision = '$Id$';

ft_nargin = nargin;

ft_nargout = nargout;

% do the general setup of the function

ft_defaults

ft_preamble init

ft_preamble debug

ft_preamble loadvar datain

ft_preamble provenance datain

% the ft_abort variable is set to true or false in ft_preamble_init

if ft_abort

return

end

% check if the input data is valid for this function

datain = ft_checkdata(datain, 'datatype', {'raw+comp', 'raw'}, 'feedback', 'yes', 'hassampleinfo', 'yes');

% get the options

cfg.method = ft_getopt(cfg, 'method', 'linear'); % default is linear, can be 'nearest', 'linear', 'spline', 'pchip', 'cubic', 'v5cubic', 'makima'

cfg.prewindow = ft_getopt(cfg, 'prewindow', 1); % default is 1 second

cfg.postwindow = ft_getopt(cfg, 'postwindow', 1); % default is 1 seconds

cfg.feedback = ft_getopt(cfg, 'feedback', 'etf');

% check if the input is valid

cfg = ft_checkopt(cfg, 'prewindow', 'numericscalar');

cfg = ft_checkopt(cfg, 'postwindow', 'numericscalar');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the actual computation is done in the middle part

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prewindow = round(cfg.prewindow * datain.fsample); % Express window in samples

postwindow = round(cfg.postwindow * datain.fsample); % Express window in samples

% Start with a copy of the input data

dataout = datain;

% Let users know that the interpolation will start and initialize the progress indicator

ntrl = numel(datain.trial);

nchan = numel(datain.label);

fprintf('Initializing %s interpolation of %d trials\n', cfg.method, ntrl);

ft_progress('init', cfg.feedback, 'Processing trial...');

for i=1:ntrl

ft_progress(i/ntrl, 'Processing trial %d from %d', i, ntrl);

tim = datain.time{i};

for j=1:nchan

dat = datain.trial{i}(j,:);

replace = isnan(dat); % Find samples that have been replaced by nans

if ~any(replace)

continue

end

onset = find(diff([0 replace])>0);

offset = find(diff([replace 0])<0);

for k=1:numel(onset)

begsample = onset(k)-prewindow;

endsample = offset(k)+postwindow;

if begsample<1

ft_warning('not enough samples for prewindow')

begsample = 1;

end

if endsample>numel(dat)

ft_warning('not enough samples for postwindow')

endsample = numel(dat);

end

x = tim(begsample:endsample);

y = dat(begsample:endsample);

xx = x; % this is where we want to know the interpolated values

x = x(~replace(begsample:endsample)); % remove the part that needs to be interpolated

y = y(~replace(begsample:endsample)); % remove the part that needs to be interpolated

yy = interp1(x, y, xx, cfg.method); % this may contain nans

% The default extrapolation behavior of INTERP1 with four input arguments is to

% extrapolate for 'spline', 'pchip' and 'makima', and to use nan for other

% methods.

if begsample==1

% there may be nans at the beginning, replace the data with mean of the values that are not nan

f = find(~isnan(yy), 1, 'first');

yy(1:f-1) = nanmean(yy);

elseif endsample==numel(dat)

% there may be nans at the end, replace the data with mean of the values that are not nan

f = find(~isnan(yy), 1, 'last');

yy(f+1:end) = nanmean(yy);

end

% insert the interpolated data

dataout.trial{i}(j,begsample:endsample) = yy;

end % for all nan-segments

end % for all channels

end % for all trials

ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cleanup

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ft_postamble debug

ft_postamble previous datain

ft_postamble provenance dataout

ft_postamble history dataout

ft_postamble savevar dataout