Численное интегрирование сигналов с абсолютным значением ⇐ Python

[Anonymous]

Предположим, у меня есть массив значений ускорения numpy s, представляющий некоторый сигнал, дискретизированный с фиксированной частотой dt. Я хочу вычислить совокупную абсолютную скорость, т.е. np.trapz(np.abs(s), dx=dt).
Это здорово, за исключением случаев, когда 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[i] += 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[i] += x0 / 2 * np.abs(waveform[i, j]) + (dt - x0) / 2 * np.abs(
waveform[i, j + 1]
)
return cav

Пример данных
Я загрузил небольшую широкополосную симуляцию движения грунта по ссылке Dropbox (около 91 МБ) для тестирования. Эти данные получены в результате конечно-разностного моделирования недавнего землетрясения недалеко от Веллингтона, Новая Зеландия, а также некоторого высокочастотного шума, полученного эмпирическим путем. Файл представляет собой HDF5, содержащий некоторые данные станции (не имеющие отношения к нашим целям) и моделируемые сигналы в ключе «waveforms». Массив имеет форму (количество станций, временные шаги, компоненты) = (1433, 5876, 3). 1d массив numpy waveform[i, :, j] — это смоделированное ускорение для i-й станции в j-м компоненте. Нам необходимо вычислить совокупную абсолютную скорость (CAV) для каждого компонента и каждой станции независимо. Тестовый код для этого можно найти ниже:

Код: Выделить всё

import time

import h5py
import numba
import numpy as np
import numpy.typing as npt

broadband_input_file = h5py.File("broadband.h5", "r")
# Load the entire dataset into memory so that the first method is not arbitrarily slowed down by file I/O
waveforms = np.array(broadband_input_file["waveforms"])
dt = 0.01

start = time.process_time()
cav_naive = np.trapz(np.abs(waveforms), dx=dt, axis=1)
print(f"CAV naive time: {time.process_time() - start}")

@numba.njit
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], waveform.shape[-1]), dtype=np.float32)
for c in range(waveform.shape[-1]):
for i in range(waveform.shape[0]):
for j in range(waveform.shape[1] - 1):
if np.sign(waveform[i, j, c]) * np.sign(waveform[i, j + 1, c]) >= 0:
cav[i, c] += (
dt
/ 2
* (np.abs(waveform[i, j, c]) + np.abs(waveform[i, j + 1, c]))
)
else:
slope = (waveform[i, j + 1, c] - waveform[i, j, c]) / dt
x0 = -waveform[i, j, c] / slope
cav[i, c] += x0 / 2 * np.abs(waveform[i, j, c]) + (
dt - x0
) / 2 * np.abs(waveform[i, j + 1, c])
return cav

# Warm up the numba compilation cache
_ = cav_integrate(waveforms, dt)

start = time.process_time()
cav_bespoke = cav_integrate(waveforms, dt)
print(f"Custom CAV time: {time.process_time() - start}")

print(
f"Ratio naive CAV / custom CAV (0, 25, 50, 75, 100% quartiles): {np.percentile(cav_naive / cav_bespoke, [0, 25, 50, 75, 100])}"
)

Что дает следующий результат

Код: Выделить всё

CAV naive time: 0.14353649699999993
Custom CAV time: 0.11182449700000019
Ratio naive CAV / custom CAV (0, 25, 50, 75, 100% quartiles): [1.00607312 1.00999796 1.01163089 1.01318455 1.02221394]

Эти различия достаточно малы, лучшие примеры более крупных различий показаны в комментариях. Некоторые из наблюдаемых сигналов имеют различия между методами на 20-40%. Даже разница в 2% может быть важна для некоторых исследователей, которых я поддерживаю. Также обратите внимание, что расчет CAV выполняется в одном потоке для сравнения, но я бы распараллелил оба метода на самом деле для самых больших массивов сигналов (с 6 или 7-кратным увеличением станций и 10-20-кратным увеличением временных шагов в зависимости от временного разрешения моделирования). . Как ни странно, параллельные накладные расходы для этого небольшого файла делают cav_integrate медленнее, чем наивный подход, если он включен.
На самом деле мы выполняем расчет CAV для всех линейных комбинаций cos(theta) * waveform [i, :, 0] + sin(theta) * waveform[i, :, 1] где theta = 0, 1,...180° для получения независимых от ориентации измерений CAV. Это одна из причин, по которой он должен быть быстрым.

Подробнее здесь: https://stackoverflow.com/questions/791 ... lute-value