Robust interpolation

1D linear interpolation despite failures of monotonicity in x or y

[Yinterp warningflag] = robustinterp(Xdata, Ydata, targetX)
Linearly interpolates in columns of X and Y values with reversals or failures of strict monotonicity in x or y (e.g. noisy and/or quantized observations of ideally monotonic relationships).
Ydata and Xdata are column vectors or matrices of the same size. TargetX is a scalar target X value, or else a vector with a target value for each data column. Robustinterp operates on each column of the data independently, to find Y values, Yinterp, corresponding to the X value(s) given in targetX.
Equalities and reversals in Xdata or Ydata are accommodated by linearly interpolating in Ydata at all points where Xdata moves to, from or through equality with targetX. The locally interpolated Y values are then simply averaged. Yinterp is not constrained to be monotonic!
Xsorted and Ysorted are the data sorted in columns so X is nondecreasing.
Each Xdata column must include at least one value less than and at least one greater than targetX: otherwise Yinterp is NaN. The optional warning flag is set to 1 where targetX coincides with the minimum (row 1) or maximum (row 2) of Xdata, otherwise 0.
Example:
x = (0:.05:1)' ;
x = x + 0.2*rand(size(x)); % disrupt monotonicity for x
y = sqrt(x);
yn = y +0.2*randn(size(y)); % ...and for y
yi = robustinterp(x,yn, 0.6)