为什么我的卷积例程与 numpy 和 scipy 不同？

Question

我想手动编写一维卷积代码，因为我正在使用内核进行时间序列分类，因此我决定制作著名的 Wikipedia 卷积图像，如图所示。

这是我的脚本。我正在使用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()

这是我得到的情节：

如您所见，由于某种原因，我的卷积发生了偏移。曲线中的数字（y 值）相同，但移动了过滤器本身大小的大约一半。

有人知道这里发生了什么吗？

Answer 1

就像您链接到的公式一样，卷积总结了从负到正无穷的索引。对于有限序列，您必须以某种方式处理不可避免的边界效应。 Numpy 和 scipy 提供了不同的方法：

numpy convolve:

模式：{‘full’、‘valid’、‘same’}、可选

scipy convolve:

模式：{'reflect','constant','nearest','mirror','wrap'}，可选

下一点是您放置原点的位置。在您提供的实现中，您从t=0 开始您的信号并丢弃负t 的加法。 Scipy 提供了一个参数origin 来考虑这一点。

origin : array_like，可选 origin 参数控制过滤器的位置。默认为 0。

您实际上可以使用 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()

总结起来，做一个卷积你必须决定如何处理序列之外的值（eq.设置为零（numpy），反射，换行，......）以及你放置的位置信号的来源。

请注意，numpy 和 scipy 的默认值也不同，它们处理边界效果（零填充与反射）的方式不同。

Answer 2

首先，为了匹配文档的符号，output[i] += signal[i - j] * ir[j] 应该是 output[i] += signal[j] * ir[i - j]

使用文档中的变量名更容易理解：

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

因为卷积可以通勤，所以你可以逃脱，所以f*g == g*f（见你的图表）

主要区别在于，长度m 信号和长度n 脉冲的“基本”卷积是长度m + n -1（参见np.convolve docs），但np.convolve( . . . , mode = 'same') 和scipy.ndimage.convolve1d 都返回通过修剪信号两端的元素来获得长度 m 的信号。

所以你的问题是你只是从右边修剪你的信号，这就是为什么

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

True

为了进行与np.convolve(..., mode = 'same') 相同的修剪，您需要：

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:]

测试：

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

True