6.16. Кривые и интегрирование

Страницы по математике охватывали операции, которые принимают массив и производят массив (или скаляр) на выходе – арифметику, редукции, трансляцию. Эта страница охватывает другой класс операций: те, которые трактуют массив как дискретизированную функцию и задают вопросы о самой функции. Интерполяцию между отсчётами, подгонку кривой к ним, интегрирование под ними, свёртку их с другим буфером.

Все они принимают входные данные ndarray и возвращают либо скаляр с плавающей точкой, либо ndarray чисел с плавающей точкой.

6.16.1. Интерполяция

interp() выполняет одномерную линейную интерполяцию. xp – это монотонно возрастающий одномерный массив независимых значений; fp – соответствующие зависимые значения; x – точки, в которых должен вычисляться интерполянт:

xp = np.array([0.0, 1.0, 2.0, 3.0])
fp = np.array([10.0, 20.0, 25.0, 30.0])
x  = np.array([0.5, 1.5, 2.5])

np.interp(x, xp, fp)
# array([15.0, 22.5, 27.5])

За пределами диапазона [xp[0], xp[-1]] результат ограничивается значениями fp[0] и fp[-1] соответственно; ключевые слова left= и right= переопределяют эти конечные точки.

Используется на камере для переотображения калибровочной таблицы в произвольные позиции отсчётов – таблицы температура-напряжение от термистора, нелинейной кривой отклика пикселей, поканальной гамма-таблицы. Один вызов библиотеки, один проход по входам, без цикла Python.

6.16.2. Подгонка и вычисление полиномов

polyfit() подгоняет полином степени deg к точкам данных (x, y) методом наименьших квадратов и возвращает коэффициенты (начиная с наивысшей степени):

x = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y = np.array([0.0, 1.1, 3.9, 9.1, 15.8])

coeffs = np.polyfit(x, y, 2)
# array([1.0, 0.0, 0.0])  (approximately)

Когда x опущен, используется range(len(y)) – полезно для быстрой подгонки полинома к регулярно дискретизированному буферу без связанной оси x:

np.polyfit(y, 2)

polyval() вычисляет полином, коэффициенты которого равны p, в точке x. Вход x может быть скаляром (возвращает float) или ndarray (возвращает ndarray):

p = np.polyfit(x, y, 2)
fitted = np.polyval(p, x)

Естественное сочетание – вызвать polyfit один раз во время калибровки, сохранить коэффициенты и вызывать polyval для вычисления получившейся кривой в каждом кадре. Шаг вычисления полинома – это несколько операций с плавающей точкой на отсчёт, что дёшево даже на самых маленьких камерах.

6.16.3. Свёртка

convolve() возвращает полноразмерную дискретную линейную свёртку двух одномерных массивов. Реализован только режим full; длина вывода равна len(a) + len(v) - 1. Срежьте результат, чтобы получить тот же эффект, что и режимы same и valid, которые предлагает настольный numpy:

a = np.array([1.0, 2.0, 3.0])
v = np.array([0.5, 0.5])

np.convolve(a, v)
# array([0.5, 1.5, 2.5, 1.5])

Полезно для коротких КИХ-фильтров и сглаживающих ядер (прямоугольное, треугольное, гауссово), где настройка цепочки SOS избыточна. Время выполнения пропорционально произведению длин двух массивов – подходит для коротких ядер, но быстро становится дороже свёртки через БПФ для длинных.

6.16.4. Интегрирование методом трапеций

trapz() интегрирует дискретизированную функцию по составной формуле трапеций:

x = np.linspace(0, np.pi, num=128)
y = np.sin(x)
np.trapz(y, x)             # ~ 2.0

Передавайте dx=, когда интервал между отсчётами равномерен и важен только шаг; передавайте x=, когда отсчёты расположены неравномерно. Это подходящий вызов для интегрирования уже захваченных данных датчика, когда аналитическая форма недоступна.

Для данных отсчётов, ограниченных по полосе (например, аудиобуфера после фильтра против наложения спектров), формула трапеций сходится как квадрат числа отсчётов, что означает, что удвоение длины буфера уменьшает ошибку в четыре раза.

Когда подынтегральное выражение – это не буфер отсчётов, а функция Python, которую приложение может вычислять в произвольных точках, подходящим вызовом будет другое семейство решателей.