An approach to represent a stochastic process by the combination of finite stochastic harmonic functions is proposed. The conditions that should be satisfied to make sure that the power spectral density function of the stochastic harmonic function process is identical to the target power spectral density are firstly studied. Then, two kinds of stochastic harmonic functions, of which the distribution of the amplitudes and the random frequencies are different, are discussed. The probabilistic characteristics of the two kinds of stochastic harmonic functions, including the asymptotic distribution, the one-dimensional probability density function, and the rate of approaching the asymptotic distribution, etc., are studied in detail by theoretical treatment and numerical examples. Responses of a nonlinear structure subjected to strong earthquake excitation are investigated. The studies show that the proposed approach can capture the target power spectral density exactly with any number of components. The reduction of the components provides flexibility and reduces the computational cost. Finally, problems that need further investigations are discussed.