The full-wave computation of transient acoustic fields with sizes in the order of 100x100x100 wavelengths by 100 periods requires a numerical method that is extremely efficient in terms of storage and computation. Iterative integral equation methods offer a good performance on these points, provided that the recurring spatiotemporal convolutions are computed with a coarse sampling and relatively few computational operations. This paper describes a method for the numerical evaluation of very large-scale, four-dimensional convolutions that employs a fast Fourier transformation and that uses a sampling rate close to or at the limit of two points per wavelength and per period. To achieve this, the functions involved are systematically filtered, windowed, and zero-padded with respect to all relevant coordinates prior to sampling. The method is developed in the context of the Neumann iterative solution of the acoustic contrast source problem for an inhomogeneous medium. The implementation of the method on a parallel computer is discussed. The obtained numerical results have a relative root mean square error of a few percent when sampling at two points per wavelength and per period. Further, the results prove that the method enables the computation of transient fields in the order of the indicated size.