why does my convolution routine differ from numpy & scipy's?

Question

I wanted to manually code a 1D convolution because I was playing around with kernels for time series classification, and I decided to make the famous Wikipedia convolution image, as seen here.

Here's my script. I'm using the standard formula for convolution for a digital signal.

import numpy as np 
import matplotlib.pyplot as plt
import scipy.ndimage

plt.style.use('ggplot')

def convolve1d(signal, ir):
    """
    we use the 'same' / 'constant' method for zero padding. 
    """
    n = len(signal)
    m = len(ir)
    output = np.zeros(n)

    for i in range(n):
        for j in range(m):
            if i - j < 0: continue
            output[i] += signal[i - j] * ir[j]

    return output

def make_square_and_saw_waves(height, start, end, n):
    single_square_wave = []
    single_saw_wave = []
    for i in range(n):
        if start <= i < end:
            single_square_wave.append(height)
            single_saw_wave.append(height * (end-i) / (end-start))
        else:
            single_square_wave.append(0)
            single_saw_wave.append(0)

    return single_square_wave, single_saw_wave

# create signal and IR
start = 40
end = 60
single_square_wave, single_saw_wave = make_square_and_saw_waves(
    height=10, start=start, end=end, n=100)

# convolve, compare different methods
np_conv = np.convolve(
    single_square_wave, single_saw_wave, mode='same')

convolution1d = convolve1d(
    single_square_wave, single_saw_wave)

sconv = scipy.ndimage.convolve1d(
    single_square_wave, single_saw_wave, mode='constant')

# plot them, scaling by the height
plt.clf()
fig, axs = plt.subplots(5, 1, figsize=(12, 6), sharey=True, sharex=True)

axs[0].plot(single_square_wave / np.max(single_square_wave), c='r')
axs[0].set_title('Single Square')
axs[0].set_ylim(-.1, 1.1)

axs[1].plot(single_saw_wave / np.max(single_saw_wave), c='b')
axs[1].set_title('Single Saw')
axs[2].set_ylim(-.1, 1.1)

axs[2].plot(convolution1d / np.max(convolution1d), c='g')
axs[2].set_title('Our Convolution')
axs[2].set_ylim(-.1, 1.1)

axs[3].plot(np_conv / np.max(np_conv), c='g')
axs[3].set_title('Numpy Convolution')
axs[3].set_ylim(-.1, 1.1)

axs[4].plot(sconv / np.max(sconv), c='purple')
axs[4].set_title('Scipy Convolution')
axs[4].set_ylim(-.1, 1.1)

plt.show()

And here's the plot I get:

As you can see, for some reason, my convolution is shifted. The numbers in the curve (y values) are identical, but shifted about half the size of the filter itself.

Anyone know what's going on here?

Answer 1

Like in the formula you are linking to, convolution summs up the indices from minus to plus infinity. For finite sequences you somehow have to deal with the boundary effects that will inevitable ocure. Numpy and scipy offer different ways to do so:

numpy convolve:

mode : {‘full’, ‘valid’, ‘same’}, optional

scipy convolve:

mode : {‘reflect’,’constant’,’nearest’,’mirror’, ‘wrap’}, optional

The next point is where you place your origin. In the implementation you provide, you start your signal at t=0 and discard summands for negative t. Scipy offers a parameter origin to take account for that.

origin : array_like, optional The origin parameter controls the placement of the filter. Default is 0.

You can actually mimic the behavior of your implementation using scipy convolve:

from scipy.ndimage.filters import convolve as convolve_sci
from pylab import *

N = 100
start=N//8
end = N-start
A = zeros(N)
A[start:end] = 1
B = zeros(N)
B[start:end] = linspace(1,0,end-start)

figure(figsize=(6,7))
subplot(411); grid(); title('Signals')
plot(A)
plot(B)
subplot(412); grid(); title('A*B numpy')
plot(convolve(A,B, mode='same'))
subplot(413); grid(); title('A*B scipy (zero padding and moved origin)')
plot(convolve_sci(A,B, mode='constant', origin=-N//2))
tight_layout()
show()

Summing it up, doing a convolution you have to decide how to treat values outside the sequence (eq. setting to zero (numpy), reflecting, wrapping, ...) and where you place the origins of the signals.

Please notice that also the defaults of numpy and scipy differ how they treat boundary effects (zero padding vs. reflecting).

Answer 2

First, to match the notation of the documentation, output[i] += signal[i - j] * ir[j] should be output[i] += signal[j] * ir[i - j]

Using the variable names in the documentation to make it simpler to understand:

i = len(signal)
for n in range(i):
    for k in range(n):
        output[n] += signal[k] * ir[n - k]

You get away with this because convolution commutes, so f*g == g*f (see your diagram)

The main difference though is that a "basic" convolution of a length m signal and a length n impulse is length m + n -1 (see np.convolve docs), but both np.convolve( . . . , mode = 'same') and scipy.ndimage.convolve1d return a length m signal by trimming elements off both ends of the signal.

So your problem is you're only trimming your signal from the right, which is why

np.all(
       np.convolve(single_square_wave, single_saw_wave)[:len(single_square_wave)]\
       ==\
       convolve1d(single_square_wave, single_saw_wave)
       )

True

In order to do the same trimming as np.convolve(..., mode = 'same'), you'd need:

def convolve1d_(signal, ir):
    """
    we use the 'same' / 'constant' method for zero padding. 
    """
    pad = len(ir)//2 - 1
    n_ = range(pad, pad + len(signal))
    output = np.zeros(pad + len(signal))

    for n in n_:
        kmin = max(0, n - len(ir) + 1)
        kmax = min(len(ir), n)
        for k in range(kmin, kmax):
            output[n] += signal[k] * ir[n - k]

    return output[pad:]

Testing:

np.all(
       np.convolve(single_square_wave, single_saw_wave, mode = 'same')\
       ==\
       convolve1d_(single_square_wave,single_saw_wave)
       )

True