Frequency-domain fatigue damage prediction based on spectral moments provides a framework in which anticipated life calculated over an entire structure subject to vibratory random loading (typically in the high cycle fatigue regime) can be rapidly obtained. However, the basis of the methods of spectral fatigue assume stationary, Gaussian, zero-mean, narrow-band (single dominant frequency) input, without the presence of overloads (stresses that exceed the initial yield stress), a significant set of restrictions. Given the importance of overloads in determining fatigue lifeI, we propose a novel “bilinear” formulation of spectral fatigue equations, that separates damage due to small and large strain amplitudes, is developed that matches or significantly outperforms existing stress-based HCF approaches (including for multiaxial elastoplastic loading) while avoiding non-conservative predictions suffered by an existing strain-based implicit formulation when the power spectral density include excursions into plastic loading due to the presence of overloads. Comparisons with synthetic and experimental data sets demonstrate the efficacy of the approach in a variety of different loading conditions.