interp_r

 INTERP_R:  A no-loop 1D interpolation routine for time series data
  xa = Original time base (e.g. decimal Julian day)
  ya = Matrix of original time series data (each time series should in a column)
  x = Evenly sampled time base on which you want to interpolate the data.
      You can specify x=n for n evenly spaced points if you prefer.

0001 function [y,xn]=interp_r(xa,ya,x)
0002 % INTERP_R:  A no-loop 1D interpolation routine for time series data
0003 %  xa = Original time base (e.g. decimal Julian day)
0004 %  ya = Matrix of original time series data (each time series should in a column)
0005 %  x = Evenly sampled time base on which you want to interpolate the data.
0006 %      You can specify x=n for n evenly spaced points if you prefer.
0007 
0008 % Rich Signell (rsignell@usgs.gov)
0009 if length(x)==1
0010     n=x;
0011 %                      generate x array
0012     dx=(xa(length(xa))-xa(1))/(n-1);
0013     x=xa(1):dx:xa(length(xa));
0014     xn=x;
0015 else
0016     n=length(x);
0017     dx=(x(n)-x(1))/(n-1);
0018 end
0019 x=x(:);
0020 xa=xa(:);
0021 
0022 %%%
0023 [mm,nn]=size(ya);
0024 if(mm==1),ya=ya.';,end
0025 [mm,nn]=size(ya);
0026 y=ones(length(x),nn);
0027 %%%
0028 ind=ceil((xa-x(1))/dx)+1;
0029 a=find(ind<1);
0030 if length(a)>0 
0031          ind(a)=[];
0032 end
0033 ibelow=zeros(size(x'));
0034 ibelow(ind)=1+length(a):length(xa);
0035 ind=find(ibelow~=0);
0036 %if(~isempty(ind)),
0037  ibelow(ind)=[ibelow(ind(1)) diff(ibelow(ind))];
0038  ibelow=cumsum(ibelow(1:n));
0039 %end
0040 a=find(ibelow<1);
0041 if length(a)>0
0042          ibelow(a)=ones(size(a));
0043 end
0044 %..........................................   now upper side
0045 ind=floor((xa-x(1))/dx)+1;
0046 a=find(ind<1);
0047 if length(a)>0 
0048          ind(a)=[];
0049 end
0050 iabove=zeros(size(x'));
0051 %ind=ind(length(ind):-1:1);
0052 iabove(ind)=1+length(a):length(xa);
0053 ind=find(iabove~=0);
0054 iabove(ind)=[iabove(ind(1)) diff(iabove(ind))];
0055 a=length(x):-1:1;
0056 iabove(a)=cumsum(iabove(a));
0057 iabove=iabove(1)+1-iabove(1:length(x));
0058 a=find(iabove>length(xa));
0059 if length(a)>0
0060          iabove(a)=ones(size(a))*length(xa);
0061 end
0062 a=find(iabove-ibelow~=0);
0063 for ii=1:nn,
0064   y(:,ii)=ya(iabove,ii);
0065   if length(a)>0
0066     fac=(x(a)-xa(ibelow(a)))./(xa(iabove(a))-xa(ibelow(a)));
0067     y(a,ii)=ya(ibelow(a),ii).*(1-fac)+ya(iabove(a),ii).*fac;
0068   end
0069   ioutside=find(x<min(xa) | x>max(xa));
0070   if(length(ioutside)>0),
0071     y(ioutside,ii)=NaN*ones(size(ioutside)); 
0072   end
0073 end