为什么我的卷积例程与 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

首先，为了匹配文档的符号，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 脉冲的 "basic" 卷积是长度 m + n -1（参见 np.convolve docs ), 但 np.convolve( . . . , mode = 'same') 和 scipy.ndimage.convolve1d return 长度 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

Answer 2

就像您链接到的公式中一样，卷积对从负到正无穷大的索引求和。对于有限序列，您必须以某种方式处理不可避免会出现的边界效应。 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 的默认值也不同，它们处理边界效应的方式不同（零填充与反射）。

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

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

python

numpy

matplotlib

convolution

scipy