Python numpy.convolve 积分限制
Python numpy.convolve integral limits
我有一个卷积积分,其形式与 post 中描述的类似:(link)
我需要从负无穷到 t 积分,而不是从 0 到 t 积分。但是,我无法使用 numpy.convolve
做到这一点,因为它总是 return 从负无穷大到正无穷大的结果。
使用 scipy.integrate.quad
会非常慢,因为我必须遍历每个 t
,并且它仅适用于具有解析表达式的被积函数。
有没有办法指定numpy.convolve
的下限和上限?非常感谢。
这是代码(很抱歉无法在此处输入 LaTeX 方程式):
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
def gaussian(tau, mu, sigma):
pdf = 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (tau - mu)**2 / (2 * sigma**2))
return pdf
def gaussian_deriv(tau, mu, sigma):
"""
derivative of gaussian function
"""
pdf = -(tau - mu)/sigma**2 * gaussian(tau, mu, sigma)
return pdf
def integral_kernel(tau):
return np.cbrt(1/tau)
def integrand(tau, t, mu, sigma):
return gaussian_deriv(tau, mu, sigma) * integral_kernel(t - tau + 1E-28)
tau = np.linspace(-7, 7, 1000)
dtau = tau[1] - tau[0]
lower_lim = tau[0]
g_deriv = gaussian_deriv(tau, mu=0, sigma=1)
result = np.zeros_like(tau)
for idx, t in np.ndenumerate(tau):
result[idx], err = integrate.quad(integrand, lower_lim, t, args=(t, mu, sigma), points=[t])
result_convolve = np.convolve(g_deriv, integral_kernel(tau), mode='same') * dtau
fig, ax = plt.subplots(2,1, figsize=(10, 6))
ax[0].plot(tau, result, 'r-', label='scipy quad')
ax[1].plot(tau, result_convolve, '.', label='numpy convolve')
ax[0].legend()
ax[1].legend()
plt.show()
如果我没理解错你想要积分
int{-inf...T} f(tau) g(t-tau) dtau
写H(t) = 1/2 (sign(t) + 1) 这和
是一样的
int{-inf...inf} f(tau) g(t-tau) [1 - H(tau-T)] dtau
或
int{-inf...inf} f^(tau) g(t-tau) dtau = f^*g(t)
f^(tau) = f(tau) * [1 - H(tau-T)].
所以你需要做的就是将T右边的f归零得到f^,然后计算常规卷积f^*g。
更新
如果t和T是同一个数字就更简单了:
int{-inf...t} f(tau)g(t-tau) dtau
在这种情况下,你可以将积分变量移动 t,所以你得到:
int{-inf...0} f(tau+t}g(-tau) dtau
接下来用 -tau 替换 tau
int{0...inf} g(tau) f(t-tau) dtau
现在你或多或少地像上面那样进行,只有 T 现在是固定的并且等于零,你不会在 T 的右边但在 T 的左边归零。
同样,这是一个相当古老的post,我仍然在这里展示我的分析。
请参阅 this post 我已经回答了 OP 引用的问题。有必要通过此 post 了解此处使用的符号和符号。
np.convolve(A, B)
给你的是
![c_{j+m} \equiv C[j] = \sum_{i = -\infty}^{+\infty} A[i] B[j - i + m]](https://latex.codecogs.com/gif.latex?c_%7Bj+m%7D&space;%5Cequiv&space;C[j]&space;=&space;%5Csum_%7Bi&space;=&space;-%5Cinfty%7D%5E%7B+%5Cinfty%7D&space;A[i]&space;B[j&space;-&space;i&space;+&space;m])
你要找的是
![c^{my}_{j+m} = \sum_{i = -\infty}^{j+m} A[i] B[j - i + m] = c_{j+m} - \sum_{i = j+m+1}^{+\infty} A[i] B[j - i + m]](https://latex.codecogs.com/gif.latex?c%5E%7Bmy%7D_%7Bj+m%7D&space;=&space;%5Csum_%7Bi&space;=&space;-%5Cinfty%7D%5E%7Bj+m%7D&space;A[i]&space;B[j&space;-&space;i&space;+&space;m]&space;=&space;c_%7Bj+m%7D&space;-&space;%5Csum_%7Bi&space;=&space;j+m+1%7D%5E%7B+%5Cinfty%7D&space;A[i]&space;B[j&space;-&space;i&space;+&space;m])
此处第二项将为零,因为 B
将负系数外推为零。
所以你看你不需要做任何事
.
我有一个卷积积分,其形式与 post 中描述的类似:(link)
我需要从负无穷到 t 积分,而不是从 0 到 t 积分。但是,我无法使用 numpy.convolve
做到这一点,因为它总是 return 从负无穷大到正无穷大的结果。
使用 scipy.integrate.quad
会非常慢,因为我必须遍历每个 t
,并且它仅适用于具有解析表达式的被积函数。
有没有办法指定numpy.convolve
的下限和上限?非常感谢。
这是代码(很抱歉无法在此处输入 LaTeX 方程式):
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
def gaussian(tau, mu, sigma):
pdf = 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (tau - mu)**2 / (2 * sigma**2))
return pdf
def gaussian_deriv(tau, mu, sigma):
"""
derivative of gaussian function
"""
pdf = -(tau - mu)/sigma**2 * gaussian(tau, mu, sigma)
return pdf
def integral_kernel(tau):
return np.cbrt(1/tau)
def integrand(tau, t, mu, sigma):
return gaussian_deriv(tau, mu, sigma) * integral_kernel(t - tau + 1E-28)
tau = np.linspace(-7, 7, 1000)
dtau = tau[1] - tau[0]
lower_lim = tau[0]
g_deriv = gaussian_deriv(tau, mu=0, sigma=1)
result = np.zeros_like(tau)
for idx, t in np.ndenumerate(tau):
result[idx], err = integrate.quad(integrand, lower_lim, t, args=(t, mu, sigma), points=[t])
result_convolve = np.convolve(g_deriv, integral_kernel(tau), mode='same') * dtau
fig, ax = plt.subplots(2,1, figsize=(10, 6))
ax[0].plot(tau, result, 'r-', label='scipy quad')
ax[1].plot(tau, result_convolve, '.', label='numpy convolve')
ax[0].legend()
ax[1].legend()
plt.show()
如果我没理解错你想要积分
int{-inf...T} f(tau) g(t-tau) dtau
写H(t) = 1/2 (sign(t) + 1) 这和
是一样的int{-inf...inf} f(tau) g(t-tau) [1 - H(tau-T)] dtau
或
int{-inf...inf} f^(tau) g(t-tau) dtau = f^*g(t)
f^(tau) = f(tau) * [1 - H(tau-T)].
所以你需要做的就是将T右边的f归零得到f^,然后计算常规卷积f^*g。
更新
如果t和T是同一个数字就更简单了:
int{-inf...t} f(tau)g(t-tau) dtau
在这种情况下,你可以将积分变量移动 t,所以你得到:
int{-inf...0} f(tau+t}g(-tau) dtau
接下来用 -tau 替换 tau
int{0...inf} g(tau) f(t-tau) dtau
现在你或多或少地像上面那样进行,只有 T 现在是固定的并且等于零,你不会在 T 的右边但在 T 的左边归零。
同样,这是一个相当古老的post,我仍然在这里展示我的分析。
请参阅 this post 我已经回答了 OP 引用的问题。有必要通过此 post 了解此处使用的符号和符号。
np.convolve(A, B)
给你的是
你要找的是
此处第二项将为零,因为 B
将负系数外推为零。
所以你看你不需要做任何事
.