  Due to its highly oscillating solution, the Helmholtz equation is numerically
challenging to solve. To obtain a reasonable solution, a mesh size that is much
smaller than the reciprocal of the wavenumber is typically required (known as
the pollution effect). High order schemes are desirable, because they are
better in mitigating the pollution effect. In this paper, we present a high
order compact finite difference method for 2D Helmholtz equations with singular
sources, which can also handle any possible combinations of boundary conditions
(Dirichlet, Neumann, and impedance) on a rectangular domain. Our method
achieves a sixth order consistency for a constant wavenumber, and a fifth order
consistency for a piecewise constant wavenumber. To reduce the pollution
effect, we propose a new pollution minimization strategy that is based on the
average truncation error of plane waves. Our numerical experiments demonstrate
the superiority of our proposed finite difference scheme with reduced pollution
effect to several state-of-the-art finite difference schemes, particularly in
the critical pre-asymptotic region where $\textsf{k} h$ is near $1$ with
$\textsf{k}$ being the wavenumber and $h$ the mesh size.
