y=kernelResample2(N,x,start,step)

function y=kernelResample2(N,x,start,step)
% y=rResampleF(N,x,start,step)
% general purpose resampler
% a matrix is resampled columnwise
if nargin ~= 4,
    error('Need four inputs');
end
if nargout ~=1,
    error('Only one output is returned.');
end
if step<=0,
    error('Step needs to be a positive number.');
end
if start<-0.5,
    error('Start needs to be greater than -0.5.');
end

load kernel;  % defines h which is a table and kernelLength
hN = numel(h);
d = max([1 step]);
kL = kernelLength*d;
toIndex = (numel(h)-1)/(kL);
toIndex1 = ceil(toIndex);
toIndex2 = toIndex1 - toIndex;

% calc starting point
pos = kL/2; pos=floor(pos)-pos+start+1;
start_ind = floor(pos)+1;
start_frac = (pos - start_ind)*toIndex;
start_frac1 = floor(start_frac);
start_frac2 = start_frac - start_frac1;

% calc length of kernel in normalized units

indN = ceil(kL);
ind=0:indN-1;

% split up step
step_ind = floor(step);
step_frac = (step - step_ind)*toIndex;
step_frac1 = floor(step_frac);
step_frac2 = step_frac-step_frac1;
% generate table index
table_index = ind*toIndex+1;

[m,n]=size(x);
y=zeros(N,n);
clear m;

for k=1:N
    ind1 = table_index-start_frac1-start_frac2;
    if ind1(indN)>hN,
        ind1(indN)=hN;
    end
    ind2 = ceil(ind1);
%    ind1 = ind2 - ind1;
%    weights = (1-ind1).*h(ind2)+ind1.*h(ind2-1);
    weights = h(ind2);
    y(k,:) =  weights * x(ind+start_ind,:);
    start_ind = start_ind + step_ind;
    start_frac1 = start_frac1 + step_frac1;
    start_frac2 = start_frac2 + step_frac2;
    if start_frac2 >= 1,
        start_frac2 = start_frac2 - 1;
        start_frac1 = start_frac1 + 1;
    end
    if start_frac1 >= 0,
        start_frac1 = start_frac1 - toIndex1;
        start_ind = start_ind + 1;
        start_frac2 = start_frac2 + toIndex2;
        if start_frac2 >= 1,
            start_frac2 = start_frac2 - 1;
            start_frac1 = start_frac1 + 1;
            % we do not need to check start_frac1 
            % because it is very low at this point
        end
    end
end