Abstract: To address the problems of large number of random variables and time-consuming computation in the simulation of non-stationary non-Gaussian stochastic processes， a fast computation method of non-stationary non-Gaussian stochastic processes is proposed based on sample interpolation by combining the stochastic harmonic function. With the known of the target evolutionary power spectrum and target density function of non-Gaussian stochastic processes， the correlation function equations of non-Gaussian stochastic processes and underlying Gaussian stochastic processes are established through Mehler’s formula， and a fast calculation method for the evolutionary power spectrum of underlying Gaussian stochastic processes is proposed through interpolation method. Subsequently， a fast simulation method for non-stationary non-Gaussian stochastic processes is proposed by combining stochastic harmonic functions， The effectiveness of this method is verified by simulating single-point uniformly modulated non-Gaussian stochastic process and multi-point non-uniformly modulated non-Gaussian stochastic processes. The results show that， when calculating the evolutionary power spectrum of the underlying Gaussian random process under the condition of ensuring accuracy， the calculation time of interpolation solution is lower than that of Mehler’s formula solution， and as the number of excitations increases， the efficiency of interpolation solution in calculating the evolutionary power spectrum of the underlying Gaussian random process is more obvious. The proposed fast computational method of non-stationary non-Gaussian stochastic processes can effectively simulate the non-Gaussian stochastic processes with the target evolutionary power spectrum and the target density function.