Это здорово, за исключением случаев, когда dt является «большим» (например, 0,01), а сигнал s одновременно длинный и часто пересекает положительные и отрицательные значения (к сожалению, частое явление), накапливается ошибка из-за того, что принимая |s| удаляет информацию об исходном знаке s. Посмотрите на картинку, чтобы лучше понять, как на самом деле выглядит эта ошибка.
У меня есть специальный код, который может правильно объяснить эту ошибку, создав модифицированное правило трапеции с numba, но есть другие очень похожие функции, которые мне нужно реализовать, например, np.trapz(np.square(s), dx=dt). Существует ли готовое решение для такого рода числового интеграла, которое:
- Можно повторно использовать для обоих интегралов np.trapz(np.square(s), dx=dt) и np.trapz(np.abs(s), dx=dt) и т. д.
- Идеально векторизован, чтобы обеспечить интеграцию можно сделать для десятков тысяч сигналов одновременно за разумное время?
@numba.njit(parallel=True)
def cav_integrate(
waveform: npt.NDArray[np.float32], dt: float
) -> npt.NDArray[np.float32]:
"""Compute the Cumulative Absolute Velocity (CAV) of a waveform."""
cav = np.zeros((waveform.shape[0],), dtype=np.float32)
for i in range(waveform.shape[0]):
for j in range(waveform.shape[1] - 1):
if np.sign(waveform[i, j]) * np.sign(waveform[i, j + 1]) >= 0:
cav += dt / 2 * (np.abs(waveform[i, j]) + np.abs(waveform[i, j + 1]))
else:
slope = (waveform[i, j + 1] - waveform[i, j]) / dt
x0 = -waveform[i, j] / slope
cav += x0 / 2 * np.abs(waveform[i, j]) + (dt - x0) / 2 * np.abs(
waveform[i, j + 1]
)
return cav
Подробнее здесь: https://stackoverflow.com/questions/791 ... lute-value