A harmonic wavelets based approximate analytical technique for determining the response evolutionary power spectrum of linear and non-linear (time-variant) oscillators endowed with fractional derivative elements is developed. Specifically, time- and frequency-dependent harmonic wavelets based frequency response functions are defined based on the localization properties of harmonic wavelets. This leads to a closed form harmonic wavelets based excitation-response relationship which can be viewed as a natural generalization of the celebrated Wiener–Khinchin spectral relationship of the linear stationary random vibration theory to account for fully non-stationary in time and frequency stochastic processes. Further, relying on the orthogonality properties of harmonic wavelets an extension via statistical linearization of the excitation-response relationship for the case of non-linear systems is developed. This involves the novel concept of determining optimal equivalent linear elements which are both time- and frequency-dependent. Several linear and non-linear oscillators with fractional derivative elements are studied as numerical examples. Comparisons with pertinent Monte Carlo simulations demonstrate the reliability of the technique.