This paper presents a method for determining the decomposition of a non-stationary, random vibration with a non-Gaussian distribution into several stationary portions with Gaussian distributions. The fatigue load of the initial non-stationary signal, and the summarized fatigue loads of the derived Gaussian portions are identical in both a direct comparison and in a comparison of the vibration responses caused by these signals on an arbitrary linear dynamical system. As stationary Gaussian signals can be described by power spectral densities (PSD), the analysis of the dynamical response with its corresponding fatigue load can be switched from time domain to frequency domain with an enormous performance gain regarding numerical computations. Especially, a combination of a realistic non-stationary signal with a large duration (from measurements), together with huge structural dynamics models is hard to analyze numerically in time domain because of enormous computation times. The proposed PSD-based method enables a very efficient numerical analysis in frequency domain. For the estimation of load spectra in frequency domain, well established PSD-based methods can be used, corresponding to Rainflow counting in time domain. This paper gives a detailed introduction to the procedure outlined with associated equations necessary for a numerical solution. The correct modeling of the fatigue damage load is proven by applying the fatigue damage spectrum concept.