Jaki byłby idealny sposób na znalezienie średniej i odchylenia standardowego sygnału dla aplikacji w czasie rzeczywistym. Chciałbym być w stanie wyzwolić kontroler, gdy sygnał był większy niż 3 odchylenie standardowe od średniej przez pewien czas.
Zakładam, że dedykowany DSP zrobiłby to dość łatwo, ale czy jest jakiś „skrót”, który może nie wymagać czegoś tak skomplikowanego?
statistics
real-time
measurement
jonsca
źródło
źródło
Odpowiedzi:
W odpowiedzi Jasona R jest błąd, który jest omawiany w „Art of Computer Programming” Knutha, t. 2. Problem pojawia się, gdy masz odchylenie standardowe, które stanowi niewielki ułamek średniej: obliczenie E (x ^ 2) - (E (x) ^ 2) cierpi z powodu dużej wrażliwości na błędy zaokrąglania zmiennoprzecinkowego.
Możesz nawet spróbować tego samodzielnie w skrypcie Python:
Otrzymuję -128,0 jako odpowiedź, co oczywiście nie jest poprawne obliczeniowo, ponieważ matematyka przewiduje, że wynik powinien być nieujemny.
Knuth przytacza podejście (nie pamiętam nazwiska wynalazcy) do obliczania średniej bieżącej i odchylenia standardowego, które wygląda mniej więcej tak:
a następnie po każdym kroku wartość
m
jest średnią, a odchylenie standardowe można obliczyć jakosqrt(S/n)
lub wsqrt(S/n-1)
zależności od ulubionej definicji odchylenia standardowego.Równanie, które piszę powyżej, jest nieco inne niż w Knuth, ale jest obliczeniowo równoważne.
Kiedy mam jeszcze kilka minut, koduję powyższą formułę w Pythonie i pokażę, że otrzymasz nieujemną odpowiedź (która, mam nadzieję, jest zbliżona do poprawnej wartości).
aktualizacja: oto jest.
test1.py:
wynik:
Zauważysz, że nadal występuje błąd zaokrąglania, ale nie jest zły, podczas gdy
naive_stats
po prostu wymiotuje.edycja: Właśnie zauważyłem komentarz Belizariusza cytujący Wikipedię, która wspomina o algorytmie Knuth.
źródło
W dziedzinie częstotliwości „wykładniczo ważona średnia bieżąca” jest po prostu prawdziwym biegunem. Jest prosty do wdrożenia w dziedzinie czasu.
Implementacja w dziedzinie czasu
Niech
mean
imeansq
być obecne szacunki średniej i średniej kwadratu sygnału. W każdym cyklu aktualizuj te szacunki o nową próbkęx
:To, co zostało wyrażone powyżej jako program imperatywny, może być również przedstawione jako diagram przepływu sygnału:
Analiza
Scalając filtry IIR we własne bloki, schemat wygląda teraz tak:
Referencje
źródło
0 > a > 1
0 < a < 1
. Jeśli twój system ma próbkowanie TmieT
i chciałbyś uśrednić stałą czasowątau
, wybierza = 1 - exp (2*pi*T/tau)
.z=1
(DC) doH(z) = a/(1-(1-a)/z)
i masz 1.Metodą, której użyłem wcześniej w aplikacji do przetwarzania wbudowanego, jest utrzymanie akumulatorów sumy i sumy kwadratów sygnału zainteresowania:
lub możesz użyć:
w zależności od preferowanej metody szacowania odchylenia standardowego . Te równania są oparte na definicji wariancji :
Używałem ich z powodzeniem w przeszłości (chociaż byłem zainteresowany jedynie oszacowaniem wariancji, a nie odchyleniem standardowym), chociaż musisz uważać na typy numeryczne, których używasz do przechowywania akumulatorów, jeśli zamierzasz sumować długi okres czasu; nie chcesz przepełnienia.
Edycja: Oprócz powyższego komentarza dotyczącego przepełnienia należy zauważyć, że nie jest to algorytm odporny numerycznie, gdy jest implementowany w arytmetyki zmiennoprzecinkowej, potencjalnie powodując duże błędy w szacowanych statystykach. Spójrz na odpowiedź Jasona S, aby uzyskać lepsze podejście w tej sprawie.
źródło
Podobnie do powyższej preferowanej odpowiedzi (Jason S.), a także wywodzącej się ze wzoru zaczerpniętego z Knuta (Vol. 2, str. 232), można również wyprowadzić formułę, aby zastąpić wartość, tj. Usunąć i dodać wartość w jednym kroku . Według moich testów, replace zapewnia lepszą precyzję niż dwuetapowa wersja „usuń / dodaj”.
Poniższy kod jest w Javie
mean
i jests
aktualizowany („globalne” zmienne składowe), tak samo jakm
is
powyżej w poście Jasona. Wartośćcount
odnosi się do rozmiaru oknan
.źródło
Odpowiedź Jasona i Nibota różni się jednym ważnym aspektem: metoda Jasona oblicza std dev i średnią dla całego sygnału (od y = 0), podczas gdy Nibota jest obliczeniem „bieżącym”, tj. Waży nowsze próbki silniejsze niż próbki z odległa przeszłość.
Ponieważ wydaje się, że aplikacja wymaga std dev i średniej w funkcji czasu, metoda Nibota jest prawdopodobnie bardziej odpowiednia (dla tej konkretnej aplikacji). Jednak prawdziwą trudną częścią będzie prawidłowe wyważenie czasu. Przykład Nibota wykorzystuje prosty filtr jednobiegunowy.
Wybór filtra dolnoprzepustowego może zależeć od tego, co wiesz o swoim sygnale i rozdzielczości czasu potrzebnej do oszacowania. Niższe częstotliwości odcięcia i wyższy porządek będą skutkować lepszą dokładnością, ale wolniejszym czasem reakcji.
Aby jeszcze bardziej skomplikować sprawę, jeden filtr stosuje się w dziedzinie liniowej, a drugi w dziedzinie do kwadratu. Kwadratowanie znacząco zmienia zawartość spektralną sygnału, więc możesz chcieć użyć innego filtra w domenie kwadratu.
Oto przykład, jak oszacować średnią, rms i std dev w funkcji czasu.
źródło
y1 = filter(a,[1 (1-a)],x);
.