function h = hilb1(F, N, b, y)

%   The function h = hilb1(F, N, b, y) computes the Hilbert transform
%   of a function F(x) defined on the real line, at specified
%   values of y (y could be a scalar, vector, or matrix.)

%   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 function F(x) has to be supplied as a string function of x, e.g.,  
%   F(x) = '1./(1+x.^2)'

%   Reference:  J.A.C. Weideman, "Computing the Hilbert Transform on
%               the Real Line", Math. Comp. Vol. 64, pp. 745--762 (1995)
 

n  = [-N:N-1]';                        %  Compute the collocation points ...
x  = b*tan(pi*(n+1/2)/(2*N));
FF = eval(F);                          %  ... and sample the function.

a  = fft(fftshift(FF.*(b-i*x)));       %  These three lines compute the
a  = exp(-i*n*pi/(2*N)).*fftshift(a);  %  expansion coefficients.
a  = flipud(i*(sign(n+1/2).*a))/(2*N);

z  = (b+i*y)./(b-i*y);                 %  The evaluation of the transform
h  = polyval(a,z)./(z.^N.*(b-i*y));    %  reduces to polynomial evaluation
                                       %  in the variable z.




















