  A new algorithm for the approximation and simulation of twofold iterated
stochastic integrals together with the corresponding L\'{e}vy areas driven by a
multidimensional Brownian motion is proposed. The algorithm is based on a
truncated Fourier series approach. However, the approximation of the remainder
terms differs from the approach considered by Wiktrosson (2001). As the main
advantage, the presented algorithm makes use of a diagonal covariance matrix
for the approximation of one part of the remainder term and has a higher
accuracy due to an exact approximation of the other part of the remainder. This
results in a significant reduction of the computational cost compared to, e.g.,
the algorithm introduced by Wiktorsson. Convergence in $L^p(\Omega)$-norm with
$p \geq 2$ for the approximations calculated with the new algorithm as well as
for approximations calculated by the basic truncated Fourier series algorithm
is proved and the efficiency of the new algorithm is analyzed.
