Confused with convolutions in scipy

Question

I tried to write my own circular convolution function in python using the fact that for two signals $f$ and $g$ we have

$$ \widehat{(f * g)} = f \cdot g $$

So I tried this

from scipy import array, zeros, signal
from scipy.fftpack import fft, ifft, convolve
def conv(f, g):
  # transform f and g to frequency domain
  F = fft(f)
  G = fft(g)
  # multiply entry-wise
  C = F * G
  # transfer C to time domain
  c = ifft(C)

  return c

However, for f = array([1,2,3,4]) and g = array([5,4,3,2]) I get

conv(f,g) = [ 34.+0.j  32.+0.j  34.+0.j  40.+0.j]
convolve.convolve(f,g) = [  1.48219694e-322   2.07507571e-322   3.06320700e-322   3.26083326e-322]
signal.convolve(f,g,'same') = [14 26 40 29]

The first function is mine, the second one comes from the fft pack and the third one is from the scipy signal package. All return different results. How does this happen and what is the "correct" one?

The circular convolution I want is

$$(f*g)[n] = \sum_{p=0}^{N-1} f[p]g[n-p] $$ where $f[p] = f_{p \mod N}$. Can somebody help me?

---

Edit: I added a direct calculation for $(f*g)$

def direct(f,g):
    r = zeros(len(f))
    for k in range(len(f)):
        for p in range(len(f)):
            r[k] = r[k] + (f[p] * g[k-p])
    return r

which returns the same result as my own implementation above. So my implementation seems correct. However, why do convolve.convolve and signal.convolve return other results?

Answer 1

I do not know what convolve.convolve does but the output of signal.convolve is the linear convolution (as opposed to circular convolution) of the two sequences. And since you supplied 'same', the output has the same size as the first input vector, and the samples are taken from the center of the complete output.

The circular convolution is indeed

34 32 34 40

and the linear convolution is

5 14 26 40 29 18 8

from which you obtained samples 2 to 5 (by using the mode 'same').

Answer 2

You're doing circular convolution, which wraps both signals around in a circle before sliding them past each other. You're convolving [..., 1, 2, 3, 4, 1, 2, 3, 4, ...] with [..., 5, 4, 3, 2, 5, 4, 3, 2, ...] in other words.

What you probably want is linear convolution, which you can get by padding f and g with zeros:

f = np.concatenate((f, np.zeros(4)))
g = np.concatenate((g, np.zeros(4)))

so that when the end wraps around and overlaps the beginning, it's multiplied by zero and has no effect on the output. Then throw away the extra zeros at the end. (This is what scipy.signal.convolve does anyway)