function y = Phi(x,n)
% y = Phi(x)
% y = Phi(x,n)
%   where n=1,2,3,...
% "smoothing function" (B-spline with knots -n/2, -n/2+1, ..., n/2 )
% for n=1: Phi=1 for |x|<.5, 0 otherwise
% for n>=1: Phi = Phi0*...*Phi0 (convolution of n fcts Phi with n=0)
%
% Phi(x,n) has integral 1, Phi(x,n) is zero outside of [-n/2,n/2]
% sum of Phi(x-j,n) over all integers j is 1
% Phi(x,n) has Fourier transform sinc(pi*xi)^n  

% Phi_d(x,n) = Phi(x/d,n)/d;
% Phi_d(x,n) has integral 1, Phi_d(x,n) is zero outside of [-d*n/2,d*n/2]
% Phi_d(x,n) has Fourier transform sinc(pi*d*xi)^n

% returns column vector

if nargin==1
  n=1;
end
x = x + n/2;
X = repmat(x(:),1,n) - (0:n-1);
y = (X>=0 & X<1);
for k=2:n
  X = X(:,1:end-1);
  y = (X.*y(:,1:end-1) + (k-X).*y(:,2:end))/(k-1);
end

