yi=cutsamples(x,y,xsi,xei)

function yi=cutsamples(x,y,xsi,xei)
% Simulates the effect of cutting a signal into samples of finite length. (An interpolation method)
% 
%  usage 1: yi=cutsamples(x,y,xsi,xei) 
%  
%       Samples the signal described by the line through x,y into pieces
%       that run from xsi to xei. yi is the mean over the intervals.
%       (xsi & xei must have equal size.)
%
%
%  usage 2: yi=cutsamples(x,y,xi,dx)
%       
%       Cut's the signal into dx-sized pieces centered around xi
%       equivalent to yi=cutsamples(x,y,xi-.5*dx,xi+.5*dx)
%       
%
% 
% Example:
%   x=sort(rand(30,1));
%   y=sin(x*5)+randn(30,1)*.3;
%   xsi=(-.2:.3:1.3)';
%   xei=xsi+.3;
%   yi=cutsamples(x,y,xei,xsi);
%   plot(x,y,'b-x',reshape([xsi xei]',[],1),reshape([yi yi]',[],1),'g',mean([xsi xei],2),yi,'go')
%   legend('measured data','samples cut')
% 
% Aslak Grinsted - Feb 2003 (updated mar2006)


[x,si,sj]=unique(x(:));
if length(si)==length(sj)
    y=y(si);
else
    warning('X values are not unique. Using mean y values.')
    ytemp=y;
    y=zeros(size(x));
    for ii=1:length(x)
       y(ii)=mean(ytemp(find(sj==ii)));
    end
end

%remove NaNs
idx=isnan(y); 
nanx=x(idx);
x(idx)=[];
y(idx)=[];

xsize=size(x);
xsi=xsi(:);
xei=xei(:);

minx=min(x);
maxx=max(x);

if length(xei)==1
    dx=abs(xei);
    xsi=xsi-.5*xei;
    xei=xsi+xei;
end


xx=unique([x;xsi;xei]); %unique is also sorting  (at least in the current version)

if (xx(1)<minx)||(xx(end)>maxx) % protect against
    dx=abs(xei(1)-xsi(1));
    x=[xx(1)-dx;x;xx(end)+dx];
    y=y([1 1:end end]);
end

yy=interp1(x,y,xx);

yint=cumtrapz(xx,yy);
ysei=interp1(xx,yint,[xsi xei]);
yi=(ysei(:,1)-ysei(:,2))./(xsi-xei);

%re-insert nans
for ii=1:length(nanx)
    yi((nanx(ii)<=xsi)&(nanx(ii)>=xei))=nan;
    %TODO: not perfect. (but pretty close)
end
yi((xei>maxx)|(xsi<minx))=nan;

if xsize(2)>1
    yi=reshape(yi,xsize); %only do it when it is necessary
end