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. In this talk I will present a recovery theory of high frequency wave fields from phase space based measurements. The construction use essentially the idea of Gaussian beams, level set description in phase space as well as the geometric optics. Our main result asserts that the kth order phase space based Gaussian beam superposition converges to the original wave field in L2 at the rate of epsilon^{k/2-n/4} in dimension n. The damage done by caustics is accurately quantified. This work is in collaboration with James Ralston (UCLA).