An approximation of the solution to a hyperbolic equation with a damping term is introduced. It is built as the composition of Fourier integral operators (FIO). We prove the convergence of this approximation in the sense of Sobolev norms as well as for the wavefront set of the solution. We apply the introduced method to numerically image seismic data.