Vibro-acoustic analysis at high frequencies is challenging due to a large sensitivity to spatial variations and short acoustical wavelengths, rendering deterministic methods expensive. At these frequencies, the wave field is usually assumed diffuse. A realization of a diffuse field can be obtained considering that the mode shapes of the system are Gaussian random fields, and that its squared eigenfrequency spacings distribution conforms to the Gaussian orthogonal ensemble (GOE) eigenvalue spacings distribution. These diffuse field properties were for example used to extend statistical energy analysis (SEA) towards energetic variance prediction. However, energetic methods such as SEA lack the propagation of phase information, which means that, e.g., time-domain reconstruction is not possible. In this contribution, instead of total energies, expressions for the ensemble average and the cross-frequency covariance of frequency response functions are presented. These expressions are obtained from the generalized definition of a diffuse field, treating the eigenfrequencies of a diffuse subsystem as a collection of points randomly located on the frequency axis, making the analysis amenable to random point process theory. These expressions are numerically validated by comparison with computationally costly detailed models where random wave scatterers are modeled explicitly.