  Computation of high frequency solutions to wave equations is important in
many applications, and notoriously difficult in resolving wave oscillations.
Gaussian beams are asymptotically valid high frequency solutions concentrated
on a single curve through the physical domain, and superposition of Gaussian
beams provides a powerful tool to generate more general high frequency
solutions to PDEs. An alternative way to compute Gaussian beam components such
as phase, amplitude and Hessian of the phase, is to capture them in phase space
by solving Liouville type equations on uniform grids. Following \cite{LR:2009}
we present a systematic construction of asymptotic high frequency wave fields
from computations in phase space for acoustic wave equations; the superposition
of phase space based Gaussian beams over two moving domains is shown necessary.
Moreover, we prove that the $k$-th order Gaussian beam superposition converges
to the original wave field in the energy norm, at the rate of
$\ep^{\frac{k}{2}+\frac{1-n}{4}}$ in dimension $n$.
