4 Correction to the some dodgy english. edited May 11 '12 at 15:45 Henry Gomersall 50622 silver badges66 bronze badges 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 the negation of$$-1$$ times the rotations of the positive frequencies (for corresponding frequency pairs). 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 the negation of the rotations of the positive frequencies. 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). 3 small correction to the period stuff. edited May 11 '12 at 15:19 Henry Gomersall 50622 silver badges66 bronze badges 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-1$$$$N$$ (the shift is implemented in the discrete Fourier domain). $$\tau$$ is in range $$[-N+1, N-1]$$. 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-1$$ (the shift is implemented in the discrete Fourier domain). $$\tau$$ is in range $$[-N+1, N-1]$$. 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]$$. 2 Changed a variable name in the code to make it more intuitive. edited May 11 '12 at 11:56 Henry Gomersall 50622 silver badges66 bronze badges 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) min_argmax_arg = optimize.fminbound(lambda tau: -correlate_point(tau), start_arg, end_arg) return min_argmax_arg  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) min_arg = optimize.fminbound(lambda tau: -correlate_point(tau), start_arg, end_arg) return min_arg  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  1 answered May 11 '12 at 11:51 Henry Gomersall 50622 silver badges66 bronze badges