We introduce a new multiscale Gaussian beam method for the
 numerical solution of the wave equation with smooth variable
 coefficients. The first computational question addressed in this
 paper is how to generate a Gaussian beam representation from general
 initial conditions for the wave equation. We propose fast multiscale
 Gaussian wavepacket transforms and introduce a highly efficient
 algorithm for generating the multiscale beam representation for a
 general initial condition. Starting from this multiscale
 decomposition of initial data, we propose the multiscale Gaussian
 beam method for solving the wave equation. The second question is
 how to perform long time propagation. Based on this new
 initialization algorithm, we utilize a simple reinitialization
 procedure that regenerates the beam representation when the beams
 become too wide. Numerical results in one, two, and three dimensions
 illustrate the properties of the proposed algorithm. The methodology
 can be readily generalized to treat other wave propagation problems.