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

模式：{ '全'， '有效'， '同'}，可选

scipy convolve：

mode：{'reflect'，'constant'，'nearest'，'mirror'，'wrap'}，可选

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

origin：array_like，optional 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() 

总结起来，做一个卷积，你必须决定如何处理序列之外的值（例如设置为零（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信号和长度为0的“基本”卷积脉冲长度为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