function f = BsplineFilter(n, s)
% compute the discrete B spline of order n, optional size s

if(nargin > 1)
   x = (-(s/2)):(s/2-1);
else
   x = (-ceil((n+1)/2)+1):(+ceil((n+1)/2)-1);
end

f = zeros(size(x));
for j= 0:n+1
   coef = (-1)^j*fact(n+1)/fact(n+1-j)/fact(j)/fact(n);
   v = max(zeros(size(x)), x+(n+1)/2-j);
   if(n == 0)
      v = [v > 0];
   else
      v = v.^n;
   end
   
   f = f + coef*v;
end

