The diagonal linewidth in two-dimensional infrared spectra is often narrower than the distribution of transition frequencies. The width along the antidiagonal is broader than predicted by the lifetime broadening. These effects arise from time-dependent fluctuations of the transition frequencies. They can be accounted for with a semiclassical approach. For systems with many coupled vibrational modes, this approach, however, becomes computationally too demanding to be practically applicable. A time-averaging approximation was suggested for linear infrared absorption spectra. In this paper, we demonstrate that the averaging can be optimized to fit a broader scale of frequency fluctuations by using a Gaussian weight function instead of the originally proposed box function. We further generalize the time-averaging method to allow the simulation of two-dimensional infrared spectra and demonstrate the method on a simple system. The approximation delivers a large speed-up of the calculation without losing significant accuracy.