function [result,kernel] = ifilter(data, N)

% Perform noise-suppressed integration of a time series that has a mean
% of zero over intervals greater than about N points. The result at any
% point will depend on N-1 data points before and behind that point.
% Use like you would cumsum() to integrate a data array. Example:
%  result = cumsum(data*dt);             % naive result
%  result = ifilter(data*dt, N);  % filtered result
%
% The filter kernel can optionally be returned, if you want it.

N = round(N); % it sometimes happens that a non-integer is input and causes problems.

Ndat_in = numel(data);
padded = [data, zeros(1,N)];  % prevent filter wraparound
Ndat = numel(padded);
kernel(2:N) = (1/2) * (N-1:-1:1)/N;
kernel(Ndat-(N-2):Ndat) = -(1/2) * (1:N-1)/N;

result = real( ifft( fft(kernel) .* fft(padded) ) );
result = result(1:Ndat_in);
   % That last part strips off the padding at the end.

% In octave, the last two lines can be combined, but Matlab won't parse this:
% result = real( ifft( fft(kernel) .* fft(padded) ) )(1:Ndat_in);