Jak uśredniać złożone odpowiedzi (i uzasadnienie)?

Zajmuję się tworzeniem oprogramowania, które oblicza odpowiedź systemu poprzez porównanie FFT sygnałów wejściowych i wyjściowych. Sygnały wejściowy i wyjściowy są podzielone na okna i dla każdego okna sygnały są odejmowane mediany i mnożone przez funkcję Hann. Reakcja instrumentu dla tego okna jest wówczas stosunkiem FFT przetwarzanych danych.

Uważam, że powyższe jest standardową procedurą, chociaż być może źle ją opisuję. Mój problem polega na tym, jak połączyć odpowiedzi z wielu okien.

O ile widzę, poprawne podejście polega na uśrednieniu wartości złożonych we wszystkich oknach. Amplituda i odpowiedź fazowa są wówczas amplitudą i fazą średniej wartości zespolonej dla każdej częstotliwości:

av_response = sum_windows(response) / n
av_amplitude = sqrt(real(av_response)**2 + imag(av_response)**2)
av_phase = atan2(imag(av_response), real(av_response))

z niejawnymi pętlami nad przedziałami częstotliwości.

Ale zostałem poproszony o zmianę tego, aby najpierw obliczyć amplitudę i fazę w każdym oknie , a następnie uśrednić amplitudy i fazy we wszystkich oknach:

amplitude = sqrt(real(response)**2 + imag(response)**2)
av_amplitude = sum_windows(amplitude) / n
phase = atan2(imag(response), real(response))
av_phase = sum_windows(phase) / n

Argumentowałem, że jest to niepoprawne, ponieważ uśrednianie kątów jest „po prostu złe” - na przykład średnia 0 i 360 stopni wynosi 180, ale ludzie, z którymi pracuję, odpowiedzieli „OK, pokażemy tylko amplitudę”.

Więc moje pytania to:

Czy mam rację sądząc, że drugie podejście jest generalnie niewłaściwe również dla amplitud?
Jeśli tak, to czy są jakieś wyjątki, które mogą być istotne i które mogą wyjaśnić, dlaczego ludzie, z którymi pracuję, wolą drugą metodę? Na przykład wygląda na to, że oba podejścia się ze sobą zgadzają, ponieważ hałas staje się mały, więc może jest to przyjęte przybliżenie niskiego poziomu hałasu?
Jeśli drugie podejście jest nieprawidłowe, czy istnieją jakieś przekonujące, autorytatywne odniesienia, których mogę użyć, aby to pokazać?
Jeśli drugie podejście jest niepoprawne, czy istnieją jakieś dobre, łatwe do zrozumienia przykłady, które pokazują to dla amplitudy (podobnie jak dla fazy 0 i 360 stopni)?
Ewentualnie, jeśli się mylę , jaka byłaby dobra książka dla mnie, aby lepiej się kształcić?

Próbowałem argumentować, że średnia -1 1 1 -1 -1 1 -1 -1 powinna wynosić zero zamiast 1, ale to nie było przekonujące. I chociaż wydaje mi się, że z czasem mogłem skonstruować argument oparty na oszacowaniu maksymalnego prawdopodobieństwa, biorąc pod uwagę konkretny model hałasu, nie jest to rodzaj rozumowania, którego słuchają ludzie, z którymi pracuję. Więc jeśli się nie mylę, potrzebuję mocnego argumentu ze strony autorytetu lub „oczywistej” demonstracji.

[Próbowałem dodać więcej tagów, ale nie mogę znaleźć odpowiednich i nie mogę zdefiniować nowych jako nowy użytkownik - przepraszam]

fft statistics frequency-response

— Andrzej Cooke
źródło

Jaki podają powód, dla którego nie podoba ci się twoja metoda?

— nibot

odpowiedź wygląda na płynniejszą, gdy wykreślona jest drugą metodą. myślę, że dzieje się tak, ponieważ w analizowanych przypadkach nie ma znaczącego sygnału (przy wyższej f), podczas gdy drugie podejście wymusza sygnał „pojawiania się” z szumu. a także różne kwestie polityczne / komunikacyjne, jak można się domyślić.

— Andrew Cooke

Czy próbowałeś podać jakieś przypadki testowe? Weź losowe dane i przefiltruj je przez niektóre filtry o znanej odpowiedzi częstotliwościowej. Sprawdź, czy oszacowanie funkcji przesyłania jest zbieżne ze znaną funkcją przesyłania.

— nibot

Nie. nie mam. to dobra sugestia. dzięki. jeśli dobrze zaprezentowane, przekonam się, że to przekonujące.

— Andrew Cooke

Odpowiedzi:

Oszacowanie funkcji przenoszenia jest zwykle realizowane nieco inaczej niż opisana metoda.

Twoja metoda jest obliczana

⟨ \frac{F [y]}{F [x]} ⟩

$\left\langle \frac{\mathcal{F}[y]}{\mathcal{F}[x]} \right\rangle$

$\langle$ $\rangle$ $\mathcal{F}$

Bardziej typowa implementacja obliczy gęstość widmową krzyża xiy podzieloną przez gęstość widmową mocy x:

\frac{⟨ F [y] \cdot F [x]^{*} ⟩}{⟨ | F [x] |^{2} ⟩} = \frac{⟨ F [y] \cdot F [x]^{*} ⟩}{⟨ F [x] \cdot F [x]^{*} ⟩}

$\frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle|\mathcal{F}[x]|^2\rangle} = \frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle\mathcal{F}[x]\cdot\mathcal{F}[x]^*\rangle}$

$\cdot$ $*$

$\mathcal{F}[x]$

Niespójne oszacowanie

Twój pracodawca zasugerował, abyś oszacował za pomocą funkcji przeniesienia

\frac{| ⟨ F [y] | ⟩}{⟨ | F [x] | ⟩}

$\frac{|\langle \mathcal{F}[y]|\rangle}{\langle |\mathcal{F}[x]|\rangle}$

To zadziała , ale ma dwie duże wady:

Nie otrzymujesz żadnych informacji o fazie.
$x$ $y$

Twoja metoda i metoda, którą opisałem, omijają te problemy, stosując spójne uśrednianie .

Bibliografia

Ogólna idea stosowania nakładających się, uśrednionych segmentów do obliczania gęstości widmowych mocy jest znana jako metoda Welcha . Uważam, że rozszerzenie zastosowania tego do szacowania funkcji przenoszenia jest często znane jako metoda Welcha, chociaż nie jestem pewien, czy wspomniano o tym w pracy Welcha. Przyglądanie się pracy Welcha może być cennym zasobem. Przydatną monografią na ten temat jest książka Bendata i Piersola „ Random Data: Analysis and Measurement Procedures” .

Uprawomocnienie

Aby sprawdzić poprawność oprogramowania, sugeruję zastosowanie kilku przypadków testowych, w których generujesz biały szum Gaussa i przepuszczasz go przez filtr cyfrowy o znanej funkcji przenoszenia. Wprowadź dane wejściowe i wyjściowe do procedury szacowania funkcji przenoszenia i sprawdź, czy oszacowanie jest zbieżne ze znaną wartością funkcji przenoszenia.

— nibot
źródło

ah! Dziękuję Ci. zbadam / spróbuję tego.

— Andrew Cooke

@nibot Jakie dokładnie długości FFT są tutaj stosowane?

— Spacey

Możesz użyć dowolnej długości. Długość określa rozdzielczość, a domyślnie (przy ustalonej ilości danych do pracy) liczbę średnich. Dłuższe fft = lepsza rozdzielczość, ale także większe błędy z powodu mniejszej liczby średnich.

— nibot

ok, inna różnica polega na tym, że masz <F (y) F * (x)> / <F (x) F * (x)> podczas gdy Phonon ma <F (y)> <F * (x)> / (< F (x)> <F * (x)>) afaict: o (

— andrew cooke

Obliczanie <F (y)> <F * (x)> / (<F (x)> <F * (x)>) nie ma sensu, ponieważ <F * (x)> zostaną natychmiast anulowane. Myślę, że to poprawne, jak to napisałem.

— nibot

Witamy w Signal Processing!

Masz całkowitą rację. Nie można po prostu uśrednić wielkości DFT i faz oddzielnie, szczególnie faz. Oto prosta demonstracja:

$z = a+bi$ $|z|$ $\angle z$ $z$

| z | = \sqrt{a^{2} + b^{2}}

$|z| = \sqrt{a^2 + b^2}$

∠ z = \tan^{- 1} (\frac{b}{a})

$\angle z = \tan ^{-1} \left( \frac{b}{a} \right)$

$z$ $z_1$ $z_2$

z = \frac{z_{1} + z_{2}}{2} = \frac{a_{1} + b_{1} i + a_{2} + b_{2} i}{2} = \frac{(a_{1} + a_{2}) + (b_{1} + b_{2}) i}{2}

$z = \frac{ z_1 + z_2 } {2} = \frac{ a_1+b_1i + a_2 + b_2i } {2} = \frac{ (a_1+a_2) + (b_1 + b_2)i } {2}$

W tym przypadku,

| z | = \sqrt{\frac{(a_{1} + a_{2})^{2}}{4} + \frac{(b_{1} + b_{2})^{2}}{4}} = \frac{1}{2} \sqrt{(a_{1} + a_{2})^{2} + (b_{1} + b_{2})^{2}} \neq \frac{\sqrt{a_{1}^{2} + b_{1}^{2}} + \sqrt{a_{2}^{2} + b_{2}^{2}}}{2}

$|z| = \sqrt{\frac{(a_1+a_2)^2}{4} + \frac{(b_1+b_2)^2}{4}} = \frac{1}{2}\sqrt{(a_1+a_2)^2 + (b_1+b_2)^2}\neq \frac{\sqrt{a_1^2 + b_1^2}+\sqrt{a_2^2 + b_2^2}}{2}$

Również,

∠ z = \frac{\tan^{- 1} (\frac{b_{1}}{a_{1}}) + \tan^{- 1} (\frac{b_{2}}{a_{2}})}{2} \neq \tan^{- 1} (\frac{2 (b_{1} + b_{2})}{2 (a_{1} + a_{2})})

$\angle z = \frac{\tan ^{-1} \left( \frac{b_1}{a_1} \right) + \tan ^{-1} \left( \frac{b_2}{a_2} \right)}{2} \neq \tan ^{-1} \left( \frac{2(b_1+b_2)}{2(a_1+a_2)} \right)$

$|z|$ $\angle z$

Teraz, aby zrobić to, co próbujesz zrobić, proponuję następujące. Teoretycznie odpowiedź impulsową systemu można znaleźć, dzieląc DFT wyniku przez DFT wejścia. Jednak w obecności hałasu uzyskasz bardzo dziwne wyniki. Nieco lepszym sposobem na to byłoby użycie dwukanałowej estymacji odpowiedzi impulsowej FFT, która wygląda następująco (nie podano tutaj pochodnej, ale można ją znaleźć online).

$G_i(f) = \dfrac{ F^1_i(f) + F^2_i(f) + \cdots + F^N_i(f) }{N}$ $F^k_i(f)$ $k$ $k$ $i$ $G_o(f) = \dfrac{ F^1_o(f) + F^2_o(f) + \cdots + F^N_o(f) }{N}$ $G$ $\hat{H}(f)$ $H(f)$

\hat{H} (f) = \frac{G_{o} (f) G_{i}^{*} (f)}{| G_{i} (f) |^{2}}

$\hat{H}(f) = \frac{G_o(f)G_i^*(f)}{|G_i(f)|^2}$

$(\cdot)^*$

— Phonon
źródło

dzięki; Nie byłem pewien, czy do głosowania ten jeden lub nibot jest jako najlepsza odpowiedź - myślę, że są one opowiadają tę samą procedurę, więc poszedł z zaleceniem książki, ale gdybym miał dwa głosy byłyby włączone to zbyt ...

— Andrew Cooke

@andrewcooke Tak, oboje opowiadają się dokładnie za tym samym. Mam nadzieję, że to wyjaśni sprawę tobie i twoim kolegom.

— Phonon

to była dla mnie ogromna pomoc (jeszcze raz dziękuję). w poniedziałek zasugeruję, aby (1) wdrożyć sugerowaną metodę i (2) dokonać porównań ze znanymi (syntetycznymi) danymi dla wszystkich trzech. wtedy mam nadzieję, że najlepsze podejście wygra: o)

— Andrew Cooke

@Phonon Jakich długości FFT używamy do obliczenia tutaj FFT? length_of_signal + max_length_of_channel + 1?

— Spacey

@Mohammad Musi być co najmniej dwa razy dłuższy niż oczekiwane opóźnienie. Wynika to z okrągłej symetrii DFT, więc otrzymasz wynik zarówno przyczynowy, jak i nie przyczynowy.

— Phonon

Jest to różnica między spójnym a niespójnym uśrednianiem widm FFT. Spójne uśrednianie jest bardziej prawdopodobne, aby odrzucić przypadkowy szum w analizie. Niespójne jest bardziej prawdopodobne, że zaakcentuje przypadkowe wielkości hałasu. Który z nich jest ważniejszy dla raportu z wyników?

— hotpaw2
źródło

jeśli dają inne wyniki, to chyba potrzebuję obiektywnej oceny. jest albo obiektywny?

— Andrew Cooke