To add a bit to the previous answers, you can get the equivalent of an upsampled band-limited cross-correlation by making your correlation variable a non-integer.

The following (python) code computes $\tau$, where
\tau = \arg \max_{\tau}\sum_{n=0}^{N-1}f\left(n\right)g\left(n+\tau\right)

That is, it finds the maximum of the cross correlation.

The input variables `a` and `b` describe $f\left(n\right)$ and $g\left(n\right)$ for $n = \{0, 1, ... , N-1\}$ and are both assumed to be band limited and periodic with period $N$ (the shift is implemented in the discrete Fourier domain). $\tau$ is in range $[-N+1, N-1]$.

The intention is to show how the cross-correlation can be performed for non-integer $\tau$, which is defined by the closure `correlate_point`. This uses the `omega` array, which describes the rotation of the complex phasor at each discrete frequency corresponding to a time-shift $\tau=1$. $\tau$ then scales this for each shift. It should be apparent that to maintain a real time signal, the rotations of the negative frequencies are just $-1$ times the rotations of the positive frequencies (for corresponding frequency pairs).

The one subtlety is in how you treat the $\frac{N}{2}$ sample (the nyquist frequency), as this is shared between the positive and negative bands. The solution used here is to interpolate between the positive rotation phasor and the negative rotation phasor (which are reflections on the real axis), which is to project either unit rotation phasor onto the real axis, which is a cos function (the `pi` is because that is value of `omega` corresponding to the nyquist frequency). Clearly this value needs to be real to maintain a real time domain signal.

You can use this to compute the cross-correlation for any arbitrarily precise value of $\tau$. Just call the closure (which can be returned as a callable) with whatever value of $\tau$ you fancy.

    import numpy
    from numpy import fft
    from scipy import optimize

    def arg_max_corr(a, b):

        if len(a.shape) > 1:
            raise ValueError('Needs a 1-dimensional array.')

        length = len(a)
        if not length % 2 == 0:
            raise ValueError('Needs an even length array.')

        if not a.shape == b.shape:
            raise ValueError('The 2 arrays need to be the same shape')

        # Start by finding the coarse discretised arg_max
        coarse_max = numpy.argmax(numpy.correlate(a, b, mode='full')) - length+1

        omega = numpy.zeros(length)
        omega[0:length/2] = (2*numpy.pi*numpy.arange(length/2))/length
        omega[length/2+1:] = (2*numpy.pi*
                (numpy.arange(length/2+1, length)-length))/length

        fft_a = fft.fft(a)

        def correlate_point(tau):
            rotate_vec = numpy.exp(1j*tau*omega)
            rotate_vec[length/2] = numpy.cos(numpy.pi*tau)

            return numpy.sum((fft.ifft(fft_a*rotate_vec)).real*b)

        start_arg, end_arg = (float(coarse_max)-1, float(coarse_max)+1)

        max_arg = optimize.fminbound(lambda tau: -correlate_point(tau), 
                start_arg, end_arg)

        return max_arg