Specialized version of the cross correlation function I am writing a specialized version of the cross correlation function as used in neuroscience. The function below is supposed to take a time series data and ask how many of its values fall in specified bins. My function xcorr works but is horrifically slow even. A test data set with 1000 points (and so 0.5*(1000*999) intervals) distributed over 400 bins takes almost ten minutes.

Bottleneck

The bottleneck is in the line counts = array([sum .... I assume it is there because each iteration of the foreach loop searches through the entire vector diffs, which is of length len(first)**2.

def xcorr(first,second,dt=0.001,window=0.2,savename=None):
    length = min(len(first),len(second))
    diffs = array([first[i]-second[j] for i in xrange(length) for j in xrange(length)])
    bins = arange(-(int(window/dt)),int(window/dt))
    counts = array[sum((diffs&gt;((i-.5)*dt)) &amp; (diffs&lt;((i+.5)*dt))) for i in bins]
    counts -= len(first)**2*dt/float(max(first[-1],second[-1])) #Normalization
    if savename:
        cPickle.dump((bins*dt,counts),open(savename,'wb'))
    return (bins,counts)
Here's a crack at my own question that uses built-in functions from NumPy 

Update

The function chain bincount(digitize(data,bins)) is equivalent to histogram, which allows for a more succinct expression.

 def xcorr(first,second,dt=0.001,window=0.2,savename=None):
      length = min(len(first),len(second))
      diffs = array([first[i]-second[j] for i in xrange(length) for j in xrange(length)])
      counts,edges = histogram(diffs,bins=arange(-window,window,dt))
      counts -= length**2*dt/float(first[length])
      if savename:
           cPickle.dump((counts,edges[:-1]),open(savename,'wb'))
      return (counts,edges[:-1])


The indexing in the return statement comes because I want to ignore the edge bins.