Nonstationary random vibration problems of nonlinear fractional systems are challenging and drawing increasing attention. This study presents an efficient numerical method for the random vibration analysis of large-scale nonlinear fractional systems subjected to fully nonstationary excitations. Firstly, the fast Fourier transform (FFT)-based frequency domain method originally proposed for the nonstationary random vibration analysis of linear integer-order systems is extended to linear fractional systems, in which the explicit expression of dynamic responses of fractional systems is derived based on the load superposition principle and Newmark method. Secondly, an equivalent linearization analysis framework is constructed to solve the nonstationary responses of hysteretic systems with fractional derivatives, in which the nonstationary response analyses of equivalent linear systems are accomplished by the FFT-based frequency domain method. The presented method can account for nonstationary excitations with arbitrary forms of evolutionary power spectrum, even of the non-separable one. Two numerical examples are used to confirm its superior performance in terms of accuracy and computational efficiency.