function G = mfs_randrespx(Gxx, Gxy, HA, HB)

# usage: G = mfs_randrespx(Gxx, Gxy, HA, HB)
#
# Input  Gxx(:, n)  Matrix with power spectral densities
#        Gxy(:, :)  Matrix with cross spectral densities
#        HA(:, n)   Matrix with transfer functions of response A
#        HB(:, n)   Matrix with transfer functions of response B
#
# Output G(:)       Power or cross spectral density of response
#
# The function computes the power or cross spectral density of the
# response from the spectral matrix of the excitation and the transfer
# functions.
#
# The columns of all matrices correspond to frequencies. The rows
# of matrices Gxx and H correspond to load cases.
#
# Matrix Gxy contains the elements above the diagonal of the spectral
# matrix. The first n-1 rows contain the elements of the first
# row, the next n-2 rows those of the second row and so on.
#
# -------------------------------------------------------------------

# Get dimensions (compatibility is already checked)

  [nofldc, nfreq] = size(Gxx);

# Compute response spectrum

  G = zeros(1, nfreq); l = 1;

  if (isempty(HB))
     for m = 1 : nofldc
         G += Gxx(m, :) .* (HA(m, :) .* conj(HA(m, :)));
         for n = m + 1 : nofldc
             G += 2 * real(HA(m, :) .* Gxy(l++, :) .* conj(HA(n, :)));
         endfor
     endfor
  else
     for m = 1 : nofldc
         G += Gxx(m, :) .* (HA(m, :) .* conj(HB(m, :)));
         for n = m + 1 : nofldc
             G += HA(m, :) .* Gxy(l, :) .* conj(HB(n, :));
             G += HA(n, :) .* conj(Gxy(l++, :)) .* conj(HB(m, :));
         endfor
     endfor
  endif

endfunction
