function h = hilb2(F, x, N, b)

%   The function h = hilb2(F, x, N, b) computes the Hilbert transform
%   of the function F(x) at the collocation points 
%    
%                x = b*tan(pi*(n+0.5)/2N),   n = -N,..., N-1.
%
%   The method is based on an expansion in rational eigenfunctions
%   of the Hilbert transform, described in the reference below.
%   The total number of terms in the expansion is 2N, and b is
%   a scalar parameter (called L in the reference.)   Note that 
%   the input F has to be supplied at the collocation points x.

%   Reference:  J.A.C. Weideman, "Computing the Hilbert Transform on
%               the Real Line", Math. Comp. Vol. 64, pp. 745--762 (1995)
 
F = F(:); x = x(:);     % Ensure that both F and x are column vectors
n = [-N:N-1]';                  
a = fftshift(fft(F.*(b-i*x)));  
h = ifft(fftshift(i*(sign(n+1/2).*a)))./(b-i*x);  



















