To accurately assess the random fatigue damage of bimodal and trimodal Gaussian process in frequency domain, a cross-mode coupling analysis method is proposed in this paper. Different from the conventional probabilistic method, the proposed method was developed based on spectral discretization method. The response spectrum was discretized into a large number of infinitesimal frequency bands, and the fatigue damage with respective to each band could be computed individually, while a coupling coefficient ξ was introduced to account for the cross-mode coupling between any two frequency components. By using the combination rule developed in this paper, the fatigue damages caused by the infinitesimal frequency bands and their couplings were assembled to obtain the total fatigue damage. Through comprehensive case studies, by taking the results of rainflow counting (RFC) method as reference and comparing with several main frequency-domain methods, it was proved that the proposed cross-mode coupling method is more accurate and robust.