Frequency-domain numerical modeling of the seismic wave equation can readily describe frequency-dependent seismic wave behaviors, yet it is computationally challenging to perform in finely discretized or large-scale geologic models. Conventional finite-difference frequency-domain (FDFD) methods for solving the Helmholtz equation usually lead to large linear systems that are difficult to solve with a direct or iterative solver. Parallel strategies and hybrid solvers can partially alleviate the computational burden by improving the performance of the linear system solver. We have developed a novel multiscale FDFD method to eventually construct a dimension-reduced linear system from the scalar Helmholtz equation based on the general framework of heterogeneous multiscale method. The methodology associated with multiscale basis functions in the multiscale finite-element method is applied to the local microscale problems of this multiscale FDFD method. Solved from frequency- and medium-dependent local Helmholtz problems, these multiscale basis functions capture fine-scale medium heterogeneities and are finally incorporated into the dimension-reduced linear system by a coupling of scalar Helmholtz problem solutions at two scales. We use several highly heterogeneous models to verify the performance in terms of the accuracy, efficiency, and memory cost of our multiscale method. The results indicate that our new method can solve the scalar Helmholtz equation in complicated models with high accuracy and quite low time and memory costs compared to the conventional FDFD methods.