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

Question

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

这是我的脚本。我正在使用数字信号卷积的标准公式。

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 提供了不同的方法来做到这一点：

\n\n

numpy 卷积：

\n\n

\n

模式：{\xe2\x80\x98full\xe2\x80\x99，\xe2\x80\x98valid\xe2\x80\x99，\xe2\x80\x98same\xe2\x80\x99}，可选

\n

\n\n

scipy 卷积：

\n\n

\n

模式 : {\xe2\x80\x98reflect\xe2\x80\x99,\xe2\x80\x99constant\xe2\x80\x99,\xe2\x80\x99nearest\xe2\x80\x99,\xe2\x80\x99mirror\xe2 \x80\x99、\xe2\x80\x98wrap\xe2\x80\x99}，可选

\n

\n\n

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

\n\n

\n

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

\n

\n\n

您实际上可以使用 scipy 模仿实现的行为convolve：

\n\n

from scipy.ndimage.filters import convolve as convolve_sci\nfrom pylab import *\n\nN = 100\nstart=N//8\nend = N-start\nA = zeros(N)\nA[start:end] = 1\nB = zeros(N)\nB[start:end] = linspace(1,0,end-start)\n\nfigure(figsize=(6,7))\nsubplot(411); grid(); title(\'Signals\')\nplot(A)\nplot(B)\nsubplot(412); grid(); title(\'A*B numpy\')\nplot(convolve(A,B, mode=\'same\'))\nsubplot(413); grid(); title(\'A*B scipy (zero padding and moved origin)\')\nplot(convolve_sci(A,B, mode=\'constant\', origin=-N//2))\ntight_layout()\nshow()\n

\n\n

\n\n

总而言之，进行卷积时，您必须决定如何处理序列之外的值（例如设置为零（numpy）、反射、环绕……）以及放置信号源的位置。

\n\n

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

\n\n

\n