Jak obliczyć widmo mocy rozmieszczone w logach?

Chciałbym obliczyć widmo mocy, w którym częstotliwości są logarytmicznie rozmieszczone.

W metodzie Welcha występuje kompromis między rozdzielczością częstotliwości wynikowego widma mocy a liczbą średnich (tj. Błąd wyniku). Chciałbym, aby ten kompromis był dynamiczny, tzn. Robiłem mniej średnich dla punktów o niskiej częstotliwości, aby uzyskać lepszą rozdzielczość przy niskiej częstotliwości.

Czy istnieje standardowy sposób to zrobić?

Przypuszczam, że jednym ze sposobów byłoby początkowo pwelchuzyskanie bardzo wysokiej rozdzielczości (małej liczby średnich), a następnie ponowne zestawienie uzyskanego widma przy użyciu binowania logarytmicznego.

Znalazłem artykuł, który odnosi się bezpośrednio do tego pytania:

• Michael Tröbs i Gerhard Heinzel, „ Ulepszone oszacowanie widma na podstawie cyfrowych szeregów czasowych na logarytmicznej osi częstotliwości ”, w pomiarach , tom 39 (2006), s. 120–129. (Lub bezpłatny wstępny wydruk ).

Pierwsze kilka liczb w dokumencie dobrze ilustruje problem, który rozwiązuje ten algorytm, a odnośniki zawierają przydatną bibliografię innych podejść (stała transformata Q, temperowana transformata Fouriera, artykuł z badania itp.).

Ich podejście nie polega na ponownym binowaniu danych wyjściowych istniejącej estymacji widma mocy opartej na FFT, ale jedynie na obliczeniu dyskretnej transformaty Fouriera przy (logarytmicznie rozmieszczonych) interesujących częstotliwościach. Dla każdej częstotliwości, która ma być oszacowana, w zasadzie implementują algorytm Welcha, ale z długością transformacji (a zatem także liczbą średnich) specjalnie wybraną dla każdej częstotliwości. Obliczenia każdego przedziału częstotliwości wykorzystują całe szeregi czasowe, ale segmentują inaczej. Wyniki mają pożądaną właściwość polegającą na tym, że rozdzielczość (szerokość przedziału) jest płynną funkcją częstotliwości, a wyniki można skalibrować jako gęstość widmową mocy lub widmo mocy.

Implementacja Matlab tutaj: https://github.com/tobin/lpsd

Ujawnienie: Autorzy tego artykułu są w tej samej instytucji co ja.

W tym przypadku użyłbym metody najmniejszych kwadratów do obliczenia częstotliwości pewnej znanej listy wartości. Najpopularniejszą metodą jest metoda Lomb. Działa dość podobnie do FFT lub DFT, ale będzie obliczać częstotliwość tylko dla określonych częstotliwości i może obsłużyć brakujące dane, gdyby to był problem. Pomysł jest następujący:

1. Wybierz listę częstotliwości ($w$
2. $w$$t_j$$X_j$ znajdziesz moc częstotliwości w następujący sposób:

$P_x(\omega) = \frac{1}{2} \left( \frac { \left[ \sum_j X_j \cos \omega ( t_j - \tau ) \right] ^ 2} { \sum_j \cos^2 \omega ( t_j - \tau ) } + \frac {\left[ \sum_j X_j \sin \omega ( t_j - \tau ) \right] ^ 2} { \sum_j \sin^2 \omega ( t_j - \tau ) } \right)$

Zauważ, że to nie będzie skalować się tak ładnie jak FFT, więc zrobiłbym to tylko wtedy, gdy liczba pożądanych częstotliwości jest znacznie mniejsza niż FFT, która byłaby wymagana do zebrania wszystkich danych.

W przeciwnym razie można wykonać metodę interpolacji lub dowolne inne ponowne próbkowanie FFT lub DFT.