Models for spatial autocorrelation and cross-correlation depend on the distance and direction separating two locations, and are constrained so that for all possible sets of locations, the covariance matrices implied from the models remain nonnegative-definite. Based on spatial correlation, optimal linear predictors can be constructed that yield complete maps of spatial fields from incomplete and noisy spatial data. This methodology is called kriging if the data are of only one variable type, and it is called cokriging if it is of two or more variable types. Historically, to satisfy the nonnegative-definite condition, cokriging has used coregionalization models for cross-variograms, even though this class of models is not very flexible. Recent research has shown that moving-average functions may be used to generate a large class of valid, flexible variogram models, and that they can also be used to generate valid cross-variograms that are compatible with component variograms. There are several problems with the moving-average approach, including large numbers of parameters and difficulties with integration. This article shows how the fast Fourier Transform (FFT) solves these problems. The flexible moving-average function that we consider is composed of many small rectangles, which eliminates the integration problem. The FFT allows us to compute the cross-variogram on a set of discrete lags; we show how to interpolate the cross-variogram for any continuous lag, which allows us to fit flexible models using standard minimization routines. Simulation examples are given to demonstrate the methods.