scipy.signal.convolve 中黎曼求和的工件

Question

小结：如何快速计算两个数组的有限卷积？

问题描述

我试图获得由

定义的两个函数 f(x), g(x) 的有限卷积为了实现这一目标，我提取了函数的离散样本，并将它们转换为长度为 steps 的数组：

xarray = [x * i / steps for i in range(steps)]
farray = [f(x) for x in xarray]
garray = [g(x) for x in xarray]

然后，我尝试使用 scipy.signal.convolve 函数计算卷积。此函数给出的结果与conv相同>这里。然而，结果与解析解有很大不同。修改算法 conv 以使用梯形规则可以得到所需的结果。

为了说明这一点，我让

f(x) = exp(-x)
g(x) = 2 * exp(-2 * x)

结果为：

这里 Riemann 代表一个简单的黎曼求和，trapezoidal 是黎曼算法的修改版本，使用梯形规则，scipy.signal.convolve 是 scipy 函数， analytical 是分析卷积。

现在让 g(x) = x^2 * exp(-x) ，结果变为：

这里“比率”是从 scipy 获得的值与分析值的比率。上面证明了这个问题不能通过重新归一化积分来解决。

问题

是否可以使用 scipy 的速度但保留梯形规则的更好结果，或者我是否必须编写 C 扩展才能达到所需的结果？

示例

只需复制并粘贴下面的代码即可查看我遇到的问题。通过增加 steps 变量可以使两个结果更加一致。我认为这个问题是由于右手黎曼和的人为因素造成的，因为积分在增加时被高估，并且在减少时再次接近解析解。

编辑：我现在已经包含了原始算法2 作为比较，它给出与 scipy.signal.convolve 函数相同的结果。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import math

def convolveoriginal(x, y):
    '''
    The original algorithm from http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
    '''
    P, Q, N = len(x), len(y), len(x) + len(y) - 1
    z = []
    for k in range(N):
        t, lower, upper = 0, max(0, k - (Q - 1)), min(P - 1, k)
        for i in range(lower, upper + 1):
            t = t + x[i] * y[k - i]
        z.append(t)
    return np.array(z) #Modified to include conversion to numpy array

def convolve(y1, y2, dx = None):
    '''
    Compute the finite convolution of two signals of equal length.
    @param y1: First signal.
    @param y2: Second signal.
    @param dx: [optional] Integration step width.
    @note: Based on the algorithm at http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
    '''
    P = len(y1) #Determine the length of the signal
    z = [] #Create a list of convolution values
    for k in range(P):
        t = 0
        lower = max(0, k - (P - 1))
        upper = min(P - 1, k)
        for i in range(lower, upper):
            t += (y1[i] * y2[k - i] + y1[i + 1] * y2[k - (i + 1)]) / 2
        z.append(t)
    z = np.array(z) #Convert to a numpy array
    if dx != None: #Is a step width specified?
        z *= dx
    return z

steps = 50 #Number of integration steps
maxtime = 5 #Maximum time
dt = float(maxtime) / steps #Obtain the width of a time step
time = [dt * i for i in range (steps)] #Create an array of times
exp1 = [math.exp(-t) for t in time] #Create an array of function values
exp2 = [2 * math.exp(-2 * t) for t in time]
#Calculate the analytical expression
analytical = [2 * math.exp(-2 * t) * (-1 + math.exp(t)) for t in time]
#Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
#Calculate the scipy convolution
sci = signal.convolve(exp1, exp2, mode = 'full')
#Slice the first half to obtain the causal convolution and multiply by dt
#to account for the step width
sci = sci[0:steps] * dt
#Calculate the convolution using the original Riemann sum algorithm
riemann = convolveoriginal(exp1, exp2)
riemann = riemann[0:steps] * dt

#Plot
plt.plot(time, analytical, label = 'analytical')
plt.plot(time, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(time, riemann, 'o', label = 'Riemann')
plt.plot(time, sci, '.', label = 'scipy.signal.convolve')
plt.legend()
plt.show()

谢谢您的宝贵时间！

Answer 1

或者，对于那些更喜欢 numpy 而不是 C 的人。它会比 C 实现慢，但只有几行。

>>> t = np.linspace(0, maxtime-dt, 50)
>>> fx = np.exp(-np.array(t))
>>> gx = 2*np.exp(-2*np.array(t))
>>> analytical = 2 * np.exp(-2 * t) * (-1 + np.exp(t))

在这种情况下，这看起来像梯形（但我没有检查数学）

>>> s2a = signal.convolve(fx[1:], gx, 'full')*dt
>>> s2b = signal.convolve(fx, gx[1:], 'full')*dt
>>> s = (s2a+s2b)/2
>>> s[:10]
array([ 0.17235682,  0.29706872,  0.38433313,  0.44235042,  0.47770012,
        0.49564748,  0.50039326,  0.49527721,  0.48294359,  0.46547582])
>>> analytical[:10]
array([ 0.        ,  0.17221333,  0.29682141,  0.38401317,  0.44198216,
        0.47730244,  0.49523485,  0.49997668,  0.49486489,  0.48254154])

最大绝对误差：

>>> np.max(np.abs(s[:len(analytical)-1] - analytical[1:]))
0.00041657780840698155
>>> np.argmax(np.abs(s[:len(analytical)-1] - analytical[1:]))
6

Answer 2

简答：用C编写！

长答案

使用关于numpy数组的食谱我用C重写了梯形卷积方法。为了使用 C 代码，需要三个文件 (https://gist.github.com/1626919）

• C 代码（performancemodule.c）。
• 用于构建代码并使其可从 python 调用的安装文件 (performancemodulesetup.py)。
• 使用 C 扩展名的 python 文件 (performancetest.py) 通过

执行以下操作，代码应在下载后运行

• 调整 performancemodule.c 中的包含路径。

• 运行以下命令

  python PerformanceModuleSetup.py 构建
  python Performancetest.py

您可能需要将库文件 performancemodule.so 或 performancemodule.dll 复制到与 performancetest.py 相同的目录中。

结果和性能

结果彼此完全一致，如下所示：

C 方法的性能甚至更好比 scipy 的卷积方法。运行数组长度为 50 的 10k 卷积

convolve (seconds, microseconds) 81 349969
scipy.signal.convolve (seconds, microseconds) 1 962599
convolve in C (seconds, microseconds) 0 87024

因此，C 实现比 python 实现快大约 1000 倍，比 scipy 实现快 20 倍多（诚然，scipy 实现比多才多艺的）。

编辑：这并不能完全解决原来的问题，但足以满足我的目的。