Time harmonic waves, modeled by the Helmholtz equation, are used in several engineering processes including, for instance, radar and seismic imaging. In the context of seismic imaging, waves propagate through the earth, which can be represented (in the simplest case) as an heterogeneous acoustic medium. Depending on the application, the simulation of high-frequency waves can be required, especially for high-resolution imaging. Numerical approximation of high-frequency waves is a challenging problem, even in homogeneous media. Indeed, numerical approximations suffer from the so-called pollution effect: if the number of discretization points per wave length is kept constant, the numerical solution diverges from the best approximation the scheme is capable of when the frequency is increasing. As a result, drastic conditions are imposed on the mesh at high frequency: the number of points per wave length must be increased when the frequency increases. In the case of homogeneous media, it has been shown (and observed numerically) that high order methods are able to reduce the pollution effect, making them cheaper than low order methods to solve for high-frequency. However, the application of high order methods to highly heterogeneous media is not trivial. It turns out that high order methods are build on coarser meshes, so that they do not capture fine scale variations of the propagation medium if the parameters are taken to be constant in each cell. The aim of this work is the study of a multiscale medium approximation strategy for high order methods. We propose a theoretical analysis of the pollution effect in the context of multiscale medium approximation when using polynomial shape functions and we present numerical experiments (including geophysical benchmarks).