A new numerical method for computing three-dimensional Stokes flow driven by a doubly-periodic array of regularized forces is presented. The method is based on deriving an analytical representation of a regularized Green's function in Fourier space. Then only an inverse fast Fourier transform (inverse FFT) has to be computed to determine the fluid velocity on a grid in the physical domain. The velocity at other points can be interpolated from this grid. Accuracy is verified by comparing numerical results to a solution that is independent of the method. Although the regularized forces lead to a smooth velocity field, the Green's function may contain rapid transitions that are not captured properly on a coarse grid. In that case, an Ewald splitting technique is used to compute the grid-resolved part of the flow using an inverse FFT and a sum in physical space for the localized part of the velocity. The splitting parameter can be chosen as small as a few grid cells, which makes the sum in physical space converge extremely fast. We present numerical examples that demonstrate that fact. In some cases, when the grid size is sufficiently small compared to the regularization parameter, the Ewald splitting is not needed.