  A practical and simple stable method for calculating Fourier integrals is
proposed, effective both at low and at high frequencies. An approach based on
the fruitful idea of Levin, to use of the collocation method to approximate the
slowly oscillating part of the antiderivative of the desired integral, allows
reducing the calculation of the integral of a highly oscillating function (with
a linear phase) to solving a system of linear algebraic equations with a
three-diagonal triangular or five-diagonal band Hermitian matrix. The choice of
Gauss-Lobatto grid nodes as collocation points makes it possible to use the
properties of discrete "orthogonality" of Chebyshev differentiation matrices in
physical and spectral spaces. This is realized by increasing the efficiency of
the numerical algorithm for solving the problem. The system pre-conditioning
procedure leads to significantly less cumbersome and more economical
calculation formulas. To avoid possible numerical instability of the algorithm,
we proceed to the solution of a normal system of linear algebraic equations.
