# Is up-sampling prior to cross-correlation useless?

Consider a simple case where two signals from two different sensors are cross-correlated, and the time-delay-of-arrival computed from the absissa of the peak of their cross-correlation function.

Now let us further assume that due to the dimensionality constraints of both antennas and the constraints on maximum possible sampling rate, the maximum attainable delay possible is $D$, corresponding to 10 samples.

The problem:

Because of those constraints, your computed delay may vary from any integer value between 0 and 10 samples, that is: $0 \le D \le 10$. This is problematic because what I really want is fractional-delay discrimination of the delay between the two signals impinging on my antennas, and changing the dimensions or the sampling rate are not an option.

Some thoughts:

• Naturally, the first thing I think of for this case is upsampling the signals before performing a cross-correlation. However I think this is 'cheating' somehow, because I am not really adding any new information into the system.

• I do not understand how upsampling is not 'cheating' in a sense. Yes, we are reconstructing our signal based on its currently observed frequency information, but how does this give one knowledge of where a signal truly started between, say, $D=7$ and $D=8$? Where was this information contained in the original signal that determined that the true fractional-delay start of the signal was actually at $D=7.751$?

The question(s):

• Is this truly 'cheating'?

• If not, then where is this new 'information' coming from?
• If yes, then what other options are available for estimating fractional-delay times?
• I am aware of upsampling the result of the cross-correlation, in an attempt to garner sub-sample answers to the delay, but is this too not also a form of 'cheating'? Why is it different from upsampling prior to the cross-correlation?

If indeed it is the case that the upsampling is not 'cheating', then why would we ever need to increase our sampling rate? (Isnt having a higher sampling rate always better in a sense than interpolating a low sampled signal?)

It would seem then that we could just sample at a very low rate and interpolate as much as we want. Would this then not make increasing the sample rate 'useless' in light of simply interpolating a signal to our heart's desire? I realize that interpolation takes computational time and simply starting with a higher sample rate would not, but is that then the only reason?

Thanks.

• I doubt there's any difference in accuracy, since the amount of information is the same either way, but it's certainly cheaper to interpolate after the cross-correlation only in the region of interest than to upsample everything first and then do all those extra multiplications. – endolith May 8 '12 at 18:12
• @endolith Good point(s). I am now clearer as to why/how this works, and yes, upsampling the result would be the way to go in this case. – Spacey May 9 '12 at 15:43

It's not cheating, and it's also not adding any new information. What you are doing is the same thing that any upsampling LPF is doing- adding zeros and then reconstructing the waveform with the already known frequency information. Thus, there is no new information, but there is still finer time resolution.

Upsampling the result is similar- no new information but finer time resolution. You can do something very similar through quadratic interpolation.

All of these methods- upsampling and polynomial interpolation- get their information on where the fractional peak is from both the peak itself and its neighbors. A quick pictorial example.

The blue line in the picture above is my simulated cross-correlation data (though it could be any result, not just a cross-correlation). It is what I call a "balanced" peak because the neighbors are symmetric. As you might expect, the resulting quadratic interpolation (red line) indicates that the true peak is at zero.

The image below, on the other hand, shows an unbalanced peak. Please note that nothing has changed in the result except for the values of the two nearest neighbors. This causes the interpolator, though, to shift its estimate of the fractional peak.

A nifty side benefit of these methods (polynomial interpolation and upsampling) is that it also gives you an estimate of the true peak value, though we are usually more interested in the location.

If indeed it is the case that the upsampling is not 'cheating', then why would we ever need to increase our sampling rate?

To satisfy the Nyquist criterion.

Isn't having a higher sampling rate always better in a sense than interpolating a low sampled signal?

No. From a theoretical standpoint, as long as the Nyquist criterion is satisfied it doesn't matter what the sample rate is. From a practical standpoint you generally go with as low a sample rate as you can get away with to reduce the storage requirements and computational load, which in turn reduces the resources that are needed and power consumption.

• @Jim, I am aware of Nyquist Criterion. :-) What I meant was in the context of increasing $F_s$, would it not be better to have a higher sample rate to begin with, so as not to have to up-sample, and from what I am hearing, this is true. It seems that once Nyquist is satisfied, all the information you are ever going to get is 'already there', and up-sampling to some level will give you whatever time-resolution you want. In that sense, does this make it open ended? The higher you can afford to upsample the better? – Spacey May 8 '12 at 19:54
• @Mohammad Yes and no. The noise- either from the results themselves or quantization noise- will eventually make increasing the time resolution meaningless. Until that point, though, yes, upsampling more should improve the accuracy of the estimate. – Jim Clay May 8 '12 at 20:03
• notes and example code for quadratic/parabolic interpolation: gist.github.com/255291#file_parabolic.md and some alternative interpolation methods: dspguru.com/dsp/howtos/how-to-interpolate-fft-peak – endolith May 8 '12 at 20:23
• @JimClay Got a chance to sleep on it. It is clear now - the information is always there - its just encoded in the relationship among the samples so to speak. And this is what the poly-fitting in fact uses in its interpolation. And since the signal is band-limited, (ie, can only change so quickly within a time span), there are only so many ways in which it can exist between samples. – Spacey May 9 '12 at 15:36

Any bandlimited signal can be interpolated. The additional information "between the samples" is contained in the adjacent samples plus the fact that the signal was bandlimited before sampling (which tends to spread information among adjacent samples). If two signals are bandlimited, than so will be the cross-correlation, so the cross-correlation can be interpolated as well. Upsampling is just another form of interpolation, a very accurate form of interpolation for bandlimited signals; but you could use Sinc interpolation as well (both of which can be more accurate than quadratic or parabolic interpolation).

Interpolation may show a peak between samples. Thus perhaps not useless.

If you have a signal containing a wider spectrum, then it can contain more information. Sampling it at a higher rate will thus provide more information, but only up to just below half the new band limit frequency, and only if the signal contained actual useful spectral frequency content above the old band limit, and if you can now obtain this additional spectrum by using a new more wide-band band-limiting process or filter, instead of the old more lossy one. Sampling data at a much higher frequency of a signal that was already bandlimited to a much lower frequency below Fs/2 will only buy you interpolation, not any more information content.

If the sampling is quantized, then sampling at a higher rate might buy you a fraction of an LSB more information, due to dithering or noise shaping of the quantization error. But this depends on the S/N ratio and accuracy of the sampler and the exact quantization process used in sampling.

If two signals are not properly bandlimited prior to sampling and cross-correlation, than not only might both upsampling or interpolation buy you a garbage result, but so might be the original non-interpolated cross-correlation.

• Thanks hotpaw2. So it is then correct to say that it doesnt really matter if you upsample both signals and then correlate, or correlate and then upsample the result? Because of the band-limiting, the two methods should give you the same results? – Spacey May 9 '12 at 15:29
• @Mohammad: I would think it contains the same information either way, but since interpolation is not perfect, the results will be slightly different depending on implementation. – endolith May 9 '12 at 17:07

I think that the best answer I can give you is: you have all the means to find out by yourself. Build an example "backwards". Using Matlab, start with two signals sampled with very small sampling periods (so that they are almost continuous-time signals). Compute the cross-correlation, and find the peak (if that's what you want), which you will be able to do with high precision. Then, downsample both signals and repeat the process. Compare the second peak's location and height with the first one. I'm sure the second will be worse. The improvement from the second to the first one is what you gain, if you upsample before cross-correlating.

To upsample the right way, both signals need to be band limited, and you need to know those bandwidths. The "new" information that you mention in your question comes from the adjacent samples, and the fact that the signals are band limited.

• Thanks Telaclavo. One thing that is really not clear to me is, the terminology of being 'band-limited'. I know what is means, but I do not understand why it is being mentioned here. ANY system, except possibly noise is 'band-limited' so why is it being mentioned over and over in this sense? – Spacey May 10 '12 at 3:25

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

• I'll need to think about whether your question is correct (and I don't have time right now). The way I think of it is that you're shifting the group delay of your signal, which you can do by any amount you like. This is exactly equivalent to circular convolution with a sinc function in the time domain, with the sinc offset by $\tau$ (but still sampled at the same locations as the original signal). It's worth noting that with an integer $\tau$, all the zero crossings of the sinc line up with sample times except for at $\tau$ (where it's $1$), which is why the signal is trivially shifted. – Henry Gomersall May 11 '12 at 15:16
• Oh, your question disappeared! Still, I'll leave the reply there. – Henry Gomersall May 11 '12 at 15:17
• Thanks Henry - (Sorry for deleting my question I was trying to make it clearer! :-) ). (I am a little slow on the python uptake but working in it). Anyway, yeah, I think I understand your method - and I think the heart of it lies in that an interpolation of the phase-response of the cross-correlation result is taken, and the corresponding $tau$ value deciphered from there. Maybe I need more time to digest, but interesting nonetheless. Where/Why have you used it instead of time-domain interpolation? The context might help. Thanks! – Spacey May 11 '12 at 15:20
• It was naive pursuit of speed, though I'm happy to be told that there is a faster time domain interpolation algorithm. My rationale is that to perform the non-integer time shift in the time domain, you'll need $N^2$ multiplications to perform the convolution, versus $N\left(\log N + 1\right)$ using the Fourier method (or thereabouts). The correlation is then taken in the time domain in both cases. If you try to resample, that strikes me as a very slow thing to be doing (and I'd do that in the Fourier domain as well!). – Henry Gomersall May 11 '12 at 15:32
• Also, I find thinking in the Fourier domain actually much simpler. But maybe that's not normal! – Henry Gomersall May 11 '12 at 15:38

There is an intuitive proof that upsampling before cross-correlation is equivalent to doing it afterwards:

Cross-correlation is convolution with the other signal time-reversed. Time reversal does not affect the bandwidth. Convolution is multiplication in the frequency domain, which also does not increase the bandwidth. If the original signals are properly band-limited, to half the sampling frequency, then so will the cross-correlation result be. No aliasing is introduced to ruin the result. Interpolation afterwards saves work.