python-3.x - 查找大型准周期阵列中的所有正向零交叉点

Question

我需要在大致周期性函数的一维数组中找到零交叉点。这将是轨道卫星向北穿过地球赤道的点。

我已经制定了一个简单的解决方案，该解决方案基于找到一个值为零或负值而下一个值为正值的点，然后使用二次或三次插值器scipy.optimize.brentq找到附近的零点。

插值器不会超出三次，在我学会使用更好的插值器之前，我首先想检查在 numpy 或 scipy 中是否已经存在一种快速方法来查找大型数组中的所有过零（n = 1E +06 至 1E+09）。

问题：所以我问在 numpy 或 scipy 中是否已经存在一种比我在这里完成的方式更快的方法来查找大型数组（n = 1E+06 到 1E+09）中的所有零交叉？

该图显示了插值零点与函数实际值之间的误差，较小的线是三次插值，较大的是二次插值。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import brentq

def f(x):
    return np.sin(x + np.sin(x*e)/e) # roughly periodic function

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
e = np.exp(1)

x = np.arange(0, 10000, 0.1)
y = np.sin(x + np.sin(x*e)/e)
crossings = np.where((y[1:] > 0) * (y[:-1] <= 0))[0]

Qd = interp1d(x, y, kind='quadratic', assume_sorted=True)
Cu = interp1d(x, y, kind='cubic', assume_sorted=True)

x0sQd = [brentq(Qd, x[i-1], x[i+1]) for i in crossings[1:-1]]
x0sCu = [brentq(Cu, x[i-1], x[i+1]) for i in crossings[1:-1]]

y0sQd = [f(x0) for x0 in x0sQd]
y0sCu = [f(x0) for x0 in x0sCu]

if True:
    plt.figure()
    plt.plot(x0sQd, y0sQd)
    plt.plot(x0sCu, y0sCu)
    plt.show()

标签： python-3.xnumpyscipysignal-processinginterpolation