title: Cross-correlation of non-periodic function with NumPy
body: <p>I have two data sets that I'm trying to cross-correlate. They look similar to the <code>arctan</code> function, so I've been using it as a model to work out how to do my signal processing.</p>
  
  <pre><code>x = linspace(-15, 15, 2**13)
  f1 = arctan(x)
  f2 = arctan(x + 2)
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/1Uwcn.png" alt="enter image description here"></p>
  
  <p>The question I need to answer is, how much do I need to shift the green signal to get it to (mostly) overlap with the blue one? I thought it would be as simple as finding the maximum in the cross-correlation function of <code>f1</code> and <code>f2</code>, and I broadly followed the advice here: <a href="https://stackoverflow.com/questions/5130808/how-to-correlate-two-time-series-with-gaps-and-different-time-bases">How to correlate two time series with gaps and different time bases?</a>. This is what I've been trying</p>
  
  <pre><code>c = correlate(f1, f2, 'full')
  s = arange(1-2**13, 2**13)
  dx = 30/2**13
  shift = s[c.argmax()]*dx
  </code></pre>
  
  <p>I would expect <code>shift</code> to equal more or less exactly 2, but in fact it's only <code>0.234</code>. This doesn't make any sense to me; I've found the x-coordinate of the maximum in cross-correlation, which should be found where there is maximal overlap of the two signals.</p>
  
  <p>Any ideas on how to calculate this quantity for this kind of function?</p>
  
  <p>EDIT: I should add, for my real data, all of the values will be between zero and one</p>
  
  <p>EDIT EDIT: The following functions are actually more like my real data:</p>
  
  <pre><code>x = linspace(-15, 15, 400)
  f1 = (arctan(-x) + pi/2) / pi
  f2 = (arctan(-x + 2) + pi/2) / pi
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/PZndi.png" alt="enter image description here"></p>
  
  <p>So using the formula given here: <a href="http://paulbourke.net/miscellaneous/correlate/" rel="nofollow noreferrer">http://paulbourke.net/miscellaneous/correlate/</a> I can write a cross-correlation function that pads the data to add ones to the left and zeros to the right:</p>
  
  <pre><code>def xcorr(x, y);
      mx = x.mean()
      my = y.mean()
      sx = x.std()
      sy = y.std()
      r = zeros(2*len(x))
  
      for d in range(-len(x), len(x)):
          csum = 0
          for i in range(0, len(x):
              yindx = i - d
              if i - d &lt; 0:
                  yval = 1
              elif i - d &gt;= len(x):
                  yval = 0
              else:
                  yval = y[yindx]
              csum += (x[i] - mx) * (yval - my)
          r[d + len(x)] = csum / (sx * sy)
      return r
  </code></pre>
  
  <p>With this function, I can now do</p>
  
  <pre><code>c = xcorr(f1, f2)
  s = arange(-400, 400)
  dx = 30/400
  shift = s[c.argmax()] * dx
  </code></pre>
  
  <p>Which comes out to 2.025, which is as close as you can get to 2 with this precision. So it looks like Jamie was correct, the issue is in how numpy <code>correlate</code> does the padding of signals.</p>
  
  <p>So, obviously my <code>xcorr</code> function is really slow as it stands. The question now is, is there a way I can make NumPy do a similar thing, or should I just proceed to write my algorithm in C using <code>ctypes</code>?</p>