为什么我的卷积例程与numpy&卑鄙的? [英] why does my convolution routine differ from numpy & scipy's?

问题描述

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

这是我的剧本.我正在使用

如您所见，由于某种原因，我的卷积发生了变化.曲线中的数字(y值)是相同的，但偏移了过滤器本身大小的一半.

有人知道这是怎么回事吗?

解决方案

就像您要链接的公式一样，卷积会将索引从负到正无穷大.对于有限序列，您必须以某种方式处理不可避免发生的边界效应.Numpy和scipy提供了不同的方法:

总结，进行卷积时，您必须决定如何处理序列外的值(例如，将值设置为零(numpy)，反射，换行等)以及放置位置信号的起源.

请注意，numpy和scipy的默认设置也不同，它们如何处理边界效果(零填充与反射).

I wanted to manually code a 1D convolution because I was playing around with kernels for time series classification, and I decided to make the famous Wikipedia convolution image, as seen here.

Here's my script. I'm using the 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()

And here's the plot I get:

As you can see, for some reason, my convolution is shifted. The numbers in the curve (y values) are identical, but shifted about half the size of the filter itself.

Anyone know what's going on here?

解决方案

Like in the formula you are linking to, convolution summs up the indices from minus to plus infinity. For finite sequences you somehow have to deal with the boundary effects that will inevitable ocure. Numpy and scipy offer different ways to do so:

numpy convolve:

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

scipy convolve:

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

The next point is where you place your origin. In the implementation you provide, you start your signal at t=0 and discard summands for negative t. Scipy offers a parameter origin to take account for that.

origin : array_like, optional The origin parameter controls the placement of the filter. Default is 0.

You can actually mimic the behavior of your implementation using 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()

Summing it up, doing a convolution you have to decide how to treat values outside the sequence (eq. setting to zero (numpy), reflecting, wrapping, ...) and where you place the origins of the signals.

Please notice that also the defaults of numpy and scipy differ how they treat boundary effects (zero padding vs. reflecting).

这篇关于为什么我的卷积例程与numpy&卑鄙的?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！