function I = mfs_qosfun1(x, Y, w)

# usage: I = mfs_qosfun1(x, Y, w)
#
# Input  x(:)     Discretization of x-axis
#        Y(:, :)  Function values (rows correspond to functions)
#        w        Circular frequency
#
# Output I(:)     Integrals
#
# The function computes the integral
#
#     I = int f(x) * exp(-i * w * x) dx
#
# using a linear interpolation of f(x) within each interval
# (first-order Filon formula).
#
# -------------------------------------------------------------------

# Initialize

  nfun = rows(Y);
  I    = zeros(nfun, 1);
  x    = x(:);

# Weights

  nx  = length(x); N = nx - 1;
  h   = diff(x);
  iw  = i * w;
  iwh = iw * h;

  v = (1 - iwh - exp(-iwh)) ./ (w^2 * h);

  w        = zeros(nx, 1);
  w(1)     = v(1);
  w(2 : N) = conj(v(1 : N - 1)) + v(2 : N);
  w(nx)    = conj(v(N));

  w = w .* exp(-iw * x);

# Integrals

  I = Y * w;
  
end
