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

Question

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

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

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

scipy convolve :

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

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

origin : array_like, optional The origin parameter controls the placement of the filter. Default is 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 的默认值也不同，它们处理边界效应的方式(零填充与反射)。