The simulation of multivariate ergodic stochastic processes is critical for structural dynamic analysis and reliability evaluation. Although the traditional spectral representation method (SRM) has a wide application in many areas, it is highly inefficient in simulating stochastic processes with many simulation points or long durations due to the significant computational cost associated with matrix factorizations concerning frequency. To address the encountered challenge, this paper presents an efficient approach for simulating ergodic stochastic processes with limited frequencies. Central to this approach is a fusion of the adaptive spectral sampling and the non-uniform fast Fourier transform (NUFFT) techniques. The adaptive spectral sampling of the envelope spectrum enables the determination of limited non-equispaced frequencies, which are randomly sampled according to a uniform distribution. Thus, the Cholesky decomposition is only required at limited specific frequencies, which dramatically reduces the computational cost of matrix factorizations. Since the randomly sampled frequencies are not equispaced, utilizing FFT to accelerate the summation of trigonometric functions becomes impractical. Then, the NUFFT that adapts the non-equispaced sampling points is employed instead to expedite this process with the non-uniform increment approximated through reduced interpolation. By taking the wind field simulation of a long-span suspension bridge as an example, a parametric analysis is conducted to investigate the effect of random frequencies on the simulation error of the developed approach and the convergence of spectra. Finally, the developed approach is further validated by focusing on the spectra and probabilistic density functions of the simulated wind samples, and the simulation performance is compared with that of the traditional approach. The analytical results demonstrate the efficiency and accuracy of the developed approach in simulating ergodic stochastic processes.