We introduce a novel technique for seismic wave extrapolation in acoustic media with variable density. The method extends two-way Fourier finite differences (FFD) from the case of constant density to variable density. We derive the new FFD operator for time marching from a pseudo-analytical solution of the acoustic-wave equation, which combines accurate spectral evaluation of spatial derivatives with a temporal iteration procedure that is exact for homogeneous media. We construct the FFD operator for coupled first-order wave propagation equations using staggered spatial and temporal grids. 2-D synthetic examples demonstrate that the FFD operator is highly accurate compared with high-order finite differences (FDs) and illustrates the capability of the operator in complicated media. The method is useful for seismic modeling applications that require high accuracy.