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. However, the usual Gaussian beam ansatz, once mapped into phase space is no longer an asymptotic solution of the wave equation. In this work we present a recovery theory of high frequency wave fields from phase space based measurements. The construction of these components use essentially the idea of Gaussian beams, that is to allow complex Hessian so that the wave is localized, and globally bounded amplitude and the Hessian can thus be obtained in phase space. We show that a superposition of these ansatz over a moving domain remains an asymptotic solution. Convergence and convergence rates are obtained. Moreover, we prove that the $k^{th}$ order phase space based Gaussian beam superposition converges to the original wave field in $L^2$ at the rate of $\ep^{\frac{k}{2}- \frac{n}{4}}$ in dimension $n$. Though some calculations are carried out only for linear Schr\"{o}dinger equations, our results and main arguments apply to more general linear wave equations. This work is in collaboration with James Ralston (UCLA).