function F = four(f, N, b, t)

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

%   The method is based on a method introduced in the reference below.
%   The total number of terms in the expansion is 2N, and b is
%   a scalar parameter.   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:  H. Weber, "Numerical computation of the Fourier Transform 
%               using Laguerre functions and the FFT", Numer. Math.  
%               Vol. 36, pp. 197--209 (1981)
 
t  = t(:);
n  = [-N:N-1]';                              % Compute the nodes ....
x  = b*tan(pi*(n+1/2)/(2*N));
ff = eval(f);                                %  ... and sample the function.

a  = fft(fftshift(ff.*(b-i*x)));             %  These two lines compute the
a  = exp(-i*n*pi/(2*N)).*fftshift(a)/(2*N);  %  expansion coefficients.

n  = [0:N-1]'; s = (-1).^n;
ap = s.*a(N+1:2*N);
an = s.*flipud(a(1:N));

pp = find(t >= 0); nn = find(t < 0);         % One has to distinguish 
tp = t(pp);        tn = t(nn);               % between positive and negative
                                             % values of t.
if length(tp) > 0
Lp = laguer(ap,2*b*tp);
Fp =  2*pi*exp(-b*tp).*Lp;
F(pp) = Fp;
end;
if length(tn) > 0 
Ln = laguer(ap,-2*b*tn);
Fn =  2*pi*exp( b*tn).*Ln; 
F(nn) = conj(Fn);
end;




















