6.16. Krzywe i całkowanie

Strony poświęcone matematyce omawiały operacje, które przyjmują tablicę i produkują tablicę (lub skalar) na wyjściu – arytmetykę, redukcje, rozpropagowanie. Ta strona obejmuje inną klasę operacji: takie, które traktują tablicę jako spróbkowaną funkcję i zadają pytania o samą funkcję. Interpolowanie między próbkami, dopasowywanie do nich krzywej, całkowanie pod nimi, splatanie ich z innym buforem.

Wszystkie te operacje przyjmują dane wejściowe typu ndarray i zwracają albo skalar zmiennoprzecinkowy, albo zmiennoprzecinkową ndarray.

6.16.1. Interpolacja

interp() wykonuje jednowymiarową interpolację liniową. xp to monotonicznie rosnąca jednowymiarowa tablica wartości niezależnych; fp to odpowiadające im wartości zależne; x to punkty, w których ma być obliczony interpolant:

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])

Poza zakresem [xp[0], xp[-1]] wynik jest odpowiednio przycinany do fp[0] oraz fp[-1]; słowa kluczowe left= i right= nadpisują te punkty końcowe.

Na kamerze używana do przemapowania tablicy kalibracyjnej na dowolne pozycje próbek – tablicy temperatura-napięcie z termistora, nieliniowej krzywej odpowiedzi pikseli, tablicy gamma dla poszczególnych kanałów. Jedno wywołanie biblioteki, jedno przejście przez dane wejściowe, bez pętli w Pythonie.

6.16.2. Dopasowanie i wyliczanie wielomianów

polyfit() dopasowuje wielomian stopnia deg do punktów danych (x, y) metodą najmniejszych kwadratów i zwraca współczynniki (od najwyższego stopnia):

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)

Gdy x jest pominięte, używane jest range(len(y)) – przydatne do szybkiego dopasowania wielomianu do regularnie próbkowanego bufora bez powiązanej osi x:

np.polyfit(y, 2)

polyval() wylicza wielomian o współczynnikach p w punkcie x. Dane wejściowe x mogą być skalarem (zwraca float) lub ndarray (zwraca ndarray):

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

Naturalnym połączeniem jest jednorazowe wywołanie polyfit w czasie kalibracji, zapamiętanie współczynników i wywoływanie polyval do wyliczania powstałej krzywej w każdej ramce. Krok wyliczania wielomianu to garść operacji zmiennoprzecinkowych na próbkę, co jest tanie nawet na najmniejszych kamerach.

6.16.3. Splot

convolve() zwraca pełnej długości dyskretny splot liniowy dwóch jednowymiarowych tablic. Zaimplementowany jest tylko tryb full; długość wyjścia wynosi len(a) + len(v) - 1. Wytnij wynik, aby uzyskać efekt trybów same i valid, które oferuje desktopowe 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])

Przydatne dla krótkich filtrów FIR i jąder wygładzających (prostokątne, trójkątne, gaussowskie), gdzie konfigurowanie łańcucha SOS to przesada. Czas wykonania jest proporcjonalny do iloczynu długości obu tablic – w porządku dla krótkich jąder, ale szybko staje się droższy niż splot przez FFT dla długich.

6.16.4. Całkowanie metodą trapezów

trapz() całkuje spróbkowaną funkcję złożoną metodą trapezów:

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

Przekaż dx=, gdy odstęp między próbkami jest jednorodny i liczy się tylko krok; przekaż x=, gdy próbki nie są równomiernie rozmieszczone. Właściwe wywołanie do całkowania już przechwyconych danych z sensora, gdy postać analityczna nie jest dostępna.

Dla danych próbek, które zostały ograniczone pasmowo (na przykład bufora audio po filtrze antyaliasingowym), reguła trapezów zbiega się jak kwadrat liczby próbek, co oznacza, że podwojenie długości bufora zmniejsza błąd czterokrotnie.

Gdy całkowana wielkość nie jest buforem próbek, lecz funkcją Pythona, którą aplikacja może obliczać w dowolnych punktach, właściwym wyborem jest inna rodzina solverów.