Finite difference methods for solving problems of time-harmonic acoustics are developed and analyzed. Multi-dimensional inhomogeneous problems with variable, possibly discontinuous, coefficients are considered, accounting for the effects of employing non-uniform grids. A weighted-average representation is less sensitive to transition in wave resolution (due to variable wave numbers or non-uniform grids) than the standard pointwise representation. Further enhancement in method performance is obtained by basing the stencils on generalizations of Padeé approximation, or generalized definitions of the derivative, reducing spurious dispersion, anisotropy, and reflection, and by improving the representation of source terms. The resulting schemes have fourth order accurate local truncation error on uniform grids and third order in the non-uniform case. Guidelines for discretization pertaining to grid orientation and resolution are presented.