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)，结果变为：

在这里输入图片描述

这里的 'ratio' 是从 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()

感谢你的时间！

scipy 数值计算数值积分伪影梯形法则有限卷积离散采样黎曼和

scipy.signal.convolve中的黎曼和伪影

问题描述

问题

示例

2 个回答

详细回答

结果和性能

撰写回答