title: Convolution computations in Numpy/Scipy
body: <p>Profiling some computational work I'm doing showed me that one bottleneck in my program was a function that basically did this (<code>np</code> is <code>numpy</code>, <code>sp</code> is <code>scipy</code>):</p>
  
  <pre><code>def mix1(signal1, signal2):
      spec1 = np.fft.fft(signal1, axis=1)
      spec2 = np.fft.fft(signal2, axis=1)
      return np.fft.ifft(spec1*spec2, axis=1)
  </code></pre>
  
  <p>Both signals have shape <code>(C, N)</code> where <code>C</code> is the number of sets of data (usually less than 20) and <code>N</code> is the number of samples in each set (around 5000). The computation for each set (row) is completely independent of any other set.</p>
  
  <p>I figured that this was just a simple convolution, so I tried to replace it with:</p>
  
  <pre><code>def mix2(signal1, signal2):
      outputs = np.empty_like(signal1)
  
      for idx, row in enumerate(outputs):
          outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same')
  
      return outputs
  </code></pre>
  
  <p>...just to see if I got the same results. But I didn't, and my questions are:</p>
  
  <ol>
  <li>Why not?</li>
  <li>Is there a better way to compute the equivalent of <code>mix1()</code>?</li>
  </ol>
  
  <p>(I realise that <code>mix2</code> probably wouldn't have been faster as-is, but it might have been a good starting point for parallelisation.)</p>
  
  <p>Here's the full script I used to quickly check this:</p>
  
  <pre><code>import numpy as np
  import scipy as sp
  import scipy.signal
  
  N = 4680
  C = 6
  
  def mix1(signal1, signal2):
      spec1 = np.fft.fft(signal1, axis=1)
      spec2 = np.fft.fft(signal2, axis=1)
      return np.fft.ifft(spec1*spec2, axis=1)
  
  def mix2(signal1, signal2):
      outputs = np.empty_like(signal1)
  
      for idx, row in enumerate(outputs):
          outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same')
  
      return outputs
  
  def test(num, chans):
      sig1 = np.random.randn(chans, num)
      sig2 = np.random.randn(chans, num)
      res1 = mix1(sig1, sig2)
      res2 = mix2(sig1, sig2)
  
      np.testing.assert_almost_equal(res1, res2)
  
  if __name__ == "__main__":
      np.random.seed(0x1234ABCD)
      test(N, C)
  </code></pre>