Stochastic response of non-linear systems to joint time-frequency non-stationary excitation has been always the focus of random vibration theory. A procedure on determination of response evolutionary power spectrum (EPS) density of nonlinear systems is presented, via the concept of statistical linearization on the time-frequency domain. Based on the newly proposed periodic generalized harmonic wavelet and its connection coefficients, the EPS relationship between the excitation and the response of linear system is obtained, by the wavelet-Galerkin solution of the linear stochastic differential equation of the motion. The response EPS dependent equivalent linear parameters of non-linear system are obtained, by employing the statistical linearization on distinct time-frequency domains. The mutually dependency between the response EPS and the equivalent linear parameters constitutes a procedure to obtain the response EPS approximately in an iterative manner. Pertinent Monte Carlo simulations demonstrate the reliability and the efficiency of the proposed approach.