In signal processing the cross-correlation (xcorr in MATLAB) is a convolution operation with one of the two sequences reversed. Since time reversal corresponds to complex conjugation in the frequency domain, you can use the DFT to compute the cross-correlation as follows:
R_xy = ifft(fft(x,N) * conj(fft(y,N)))
where N = size(x) + size(y) - 1 (preferably rounded up to a power of 2) is the length of the DFT.
Multiplication of DFTs is equivalent to circular convolution in time. Zero padding both vectors to length N keeps the circularly shifted components of y from overlapping with x, which makes the result identical to the linear convolution of x and time reversed y.
A lag of 1 is a right circular shift of y, while a lag of -1 is a left circular shift. The cross-correlation is simply the sequence of dot products for all lags. Based on standard fft ordering, these will be in an array that can be accessed as follows. Indices 0 through size(x)-1 are the positive lags. Indices N-size(y)+1 to N-1 are the negative lags in reverse order. (In Python the negative lags can be accessed conveniently with negative indices such as R_xy[-1].)
You can think of the zero-padded x and y as N-dimensional vectors. The dot product of x and y for a given lag is
|x|*|y|*cos(theta). The norms of x and y are constant for circular shifts, so dividing them out leaves just the varying cosine of the angle theta. If x and y (for a given lag) are orthogonal in N-space, the correlation is 0 (i.e. theta = 90 degrees). If they're co-linear, the value is either 1 (positively correlated) or -1 (negatively correlated, i.e. theta = 180 degrees). This leads to the cross-correlation normalized to unity:
R_xy = ifft(fft(x,N) * conj(fft(y,N))) / (norm(x) * norm(y))
This can be made unbiased by recomputing the norms for just the overlapping parts, but then you may as well do the entire computation in the time domain. Also, you'll see different versions of normalization. Instead of being normalized to unity, sometimes the cross-correlation is normalized by M (biased), where M = max(size(x), size(y)), or M-|m| (an unbiased estimate of the mth lag).
For maximum statistical significance the mean (DC bias) should be removed before computing the correlation. This is called the cross-covariance (xcov in MATLAB):
x2 = x - mean(x)
y2 = y - mean(y)
phi_xy = ifft(fft(x2,N) * conj(fft(y2,N))) / (norm(x2) * norm(y2))