A concept of Petrov–Galerkin enrichment which is appropriate for highly accurate and stable interpolation of variational solutions is introduced. In the finite element context, the setting refers to standard trial functions for the solution, while the test space will be enriched. The FEM interpolation procedure that we propose will be justified by local wavelets with vanishing moments based on Gegenbauer polynomials. For the reference Helmholtz equation, the continuous piecewise polynomial test functions are enriched using dispersion analysis on uniform meshes in 2d and 3d. From a-priori and a-posteriori numerical analysis it follows that the Petrov–Galerkin based enrichment approximates the exact interpolate solution of the Helmholtz equation with at least seventh order of accuracy.