function [xI, yI] = smooth(x, y, std_y, xI)
% SMOOTH does smoothed data interpolation
% If the data vector is small enough then this is done by a call to sspline
% which fits a smoothed cubic spline. However, for large vector sspline becomes
% unacceptably slow. In this case smooth does a series of smoothing spline fits
% using smaller sections of the data. These 'windows' are chosen to overlap to
% minimize the obvious 'joins' on the final smoothed curve'. The inputs and
% outputs are like those in sspline but the xI may actually equal the end points
% of x. Also the xI and x will be sorted.
%
% Usage: 
%         [xI, yI] = smooth(x, y, std_y, xI)
%
% where the inputs are: 
%          
%             x and y:  the given data pairs, in which x is not 
%      		        necessary equally spaced, and y may contain errors; 
%
%             std_y:    an estimate of the standard deviation of y;
%
%             xI:       the inserted x points where you want y's value to be
%                       interpolated; xI can be a vector too, as is often the
%                       case. Unlike SSPLINE, the first element of xI, xI(1) may
%                       be equal to the first point of given x, x(1), and the
%                       last point of xI, say, xI(N), may be equal to the last
%                       point of given x, say, x(n).
%
%       the outputs are:
%     
%              xI:      a sorted version of the above
%              yI:      the interpolated values at xI

% $Id: smooth.m,v 1.1 1998/10/30 06:43:04 mansbrid Exp $
% Copyright J. V. Mansbridge, CSIRO, Thu Oct 29 17:06:35 EST 1998

merge_type = 'linear';
[x, ind] = sort(x(:));
y = y(ind);
y = y(:);
xI = sort(xI(:));

if xI(1) < x(1)
  error('xI must always be >= the smallest value of x')
elseif xI(1) == x(1)
  x = [2*x(1) - x(2); x];
  y = [2*y(1) - y(2); y];
end
if xI(end) > x(end)
  error('xI must always be <= the largest value of x')
elseif xI(end) == x(end)
  x = [x; 2*x(end) - x(end-1)];
  y = [y; 2*y(end) - y(end-1)];
end


len_x = length(x);
window = 100;
overlap = 20;
wmo = window - overlap;

if len_x <= 2*window
  [xI, yI, a, b, c, d] = sspline(x, y, std_y, xI);
  return
end

num_slabs = ceil((len_x - overlap)/wmo);
st_ind = 1 + (0:num_slabs-1)*wmo;
fin_ind = st_ind + window - 1;
fin_ind(end) = len_x;
y_vals = repmat(NaN, length(xI), num_slabs);
for ii = 1:num_slabs
  ran = st_ind(ii):fin_ind(ii);
  ff = find((xI > x(st_ind(ii))) & (xI < x(fin_ind(ii))));
  xx = xI(ff);
  [xret, yret, a, b, c, d] =  sspline(x(ran), y(ran), std_y, xx);
  y_vals(ff, ii) = yret;
end

% merge the splines in the overlap region.

switch merge_type
  case 'linear'
    yI = y_vals(:, 1);
    ff_old = find(~isnan(y_vals(:, 1)));
    for ii = 2:num_slabs
      ff = find(~isnan(y_vals(:, ii)));
      yI(ff) = y_vals(ff, ii);
      ff_i = intersect(ff, ff_old);
      len_xx = length(ff_i);
      xx = (1:len_xx)'/(len_xx + 1);
      yI(ff_i) = xx.*y_vals(ff_i, ii) + (1 - xx).*y_vals(ff_i, ii-1);
      ff_old = ff;
    end
  case 'mean'
    yI = nanmean(y_vals, 2);
end
