Stochastic simulations are increasingly used to describe complex systems with uncertainties. To better characterize the uncertainty of such a simulation, we propose a random Fourier features method to fit fully Bayesian Gaussian process emulators for these simulations. The random Fourier features technique uses low-dimensional features of the correlation function of the Gaussian process to achieve a low-rank approximation of the correlation matrix. We prove the convergence of the proposed random Fourier features method. Simulation results show that the proposed method can significantly reduce computation time while maintaining high prediction accuracy. The advantages of the proposed method are also illustrated using a modern subway simulation.