爲什麼我的卷積程序與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