Oblicz i interpretuj częstotliwość chwilową

Znam zasadę obliczania częstotliwości chwilowej i wpadłem na wiele pytań na jej temat. Znajdziesz je wszystkie na liście punktorów na końcu tego tekstu. Tekst może być trochę długi, przepraszam za to, ale naprawdę próbowałem samodzielnie rozwiązać ten problem.

Więc interesuje mnie częstotliwość chwilowa $f(t)$ sygnału o wartościach rzeczywistych $x(t)$. Obliczenia wykonuje się za pomocą sygnału analitycznego$z(t) = x(t) + j y(t)$, gdzie $y(t)$ jest transformacją Hilberta $x(t)$.

Aby obliczyć chwilowe częstotliwości z sygnału analitycznego $z(t)$ Śledziłem artykuł:

Obliczanie częstotliwości chwilowej i chwilowej przepustowości przez Arthura E. Barnsa z 1992 r. W tym artykule wprowadza wiele metod obliczania częstotliwości chwilowej. Zapisuję wszystkie formuły, które zaproponował (i użyłem) za chwilę.

W celu „uczenia się” bawiłem się w MATLAB bardzo prostymi i dwoma nieco bardziej złożonymi sygnałami i chciałem uzyskać ich chwilowe częstotliwości.

Fs = 1000;                                            % sampling-rate = 1kHz
t = 0:1/Fs:10-1/Fs;                                    % 10s 'Timevector'
chirp_signal = chirp(t,0,1,2);                         % 10s long chirp-signal, signal 1
added_sinusoid = chirp_signal + sin(2*pi*t*10);        % chirp + sin(10Hz), signal 2
modulated_sinusoid = chirp_signal .* sin(2*pi*t*10);   % chirp * sin(10Hz), signal 3

Wykresy w dziedzinie czasu tych trzech sygnałów wyglądają następująco:

Wykresy wszystkich chwilowych częstotliwości, które otrzymałem po zastosowaniu wszystkich metod z papieru, są następujące:

Chwilowe częstotliwości czystego sygnału ćwierkającego: Chwilowe częstotliwości sygnału ćwierkającego z dodanym sinusoidą: Chwilowe częstotliwości modulowanego sygnału ćwierkającego: Należy pamiętać, że na wszystkich trzech obrazach oś y wykresu 3 i 4 jest powiększona, więc amplitudy tych sygnały są bardzo małe!

Pierwszą możliwością przejścia z sygnału analitycznego na częstotliwość chwilową jest:

$$\begin{equation} f_{2}(t) = \frac{1}{2\pi} \frac{d}{dt}\theta(t) \end{equation}$$ gdzie $\theta(t)$jest chwilową fazą. Myślę, że jest to obecnie najpopularniejsza metoda, przynajmniej tak jest obliczana na stronie MATLAB . Kod wygląda następująco:

function [instantaneous_frequency] = f2(analytic_signal,Fs)
    factor =  Fs/(2*pi);
    instantaneous_frequency = factor * diff(unwrap(angle(analytic_signal)));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

W artykule Barns sugeruje teraz (a raczej wspomnianą kompilację) cztery inne sposoby obliczania chwilowych częstotliwości z sygnału analitycznego. Wspomina także górną formułę, ale jest zdania, że ​​jest ona niepraktyczna ze względu na niejasności w fazie. Wydaje mi się, że nie znał unwrap()metody ani, mówiąc ściślej, matematyki za nią. (Ja sam dowiedziałem się o tej metodzie dopiero dziś, kiedy patrzę na inne kody źródłowe na chwilowych częstotliwościach)

W jego pracy formuła ma etykietę Numer (2), dlatego nadałem f (t) indeks 2. Wszystkie pozostałe indeksy odpowiadają w ten sam sposób ich liczbom w pracy.

Z powodu niejasności fazowych sugeruje raczej:

$$\begin{equation} f_{3}(t) = \frac{1}{2\pi}\cdot\frac{\overbrace{x(t)y'(t)}^\text{a}-\overbrace{x'(t)y(t)}^{\text{b}}}{\underbrace{x(t)^2}_{\text{c}}+\underbrace{y(t)^2}_{\text{d}}} \end{equation}$$ Wprowadziłem symbole „a”, „b”, „c” i „d”, aby ułatwić programowanie:

function [instantaneous_frequency] = f3(analytic_signal,Fs,T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    diff_x = diff(x);
    diff_y = diff(y);
    factor = Fs/(2*pi);
    a = x(2:end).*diff_y;
    b = y(2:end).*diff_x;
    c = x(2:end).^2;
    d = y(2:end).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Następnie Barner podaje jeszcze trzy formuły, które nazywa „chwilowymi przybliżeniami częstotliwości”:

$$\begin{equation} f_{9}(t) = \frac{1}{2\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{x(t)x(t+T)}_{\text{c}}+\underbrace{y(t)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function[instantaneous_frequency] = f9(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(2*pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = x(1:end-T).*x(1+T:end);
    d = y(1:end-T).*y(1+T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append 0 to return-vector to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{11}(t) = \frac{1}{4\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t-T)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t-T)}^{\text{b}}}{\underbrace{x(t-T)x(t+T)}_{\text{c}}+\underbrace{y(t-T)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = f11(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(4*pi*T);
    a = x(1:end-2*T).*y(1+2*T:end);
    b = x(1+2*T:end).*y(1:end-2*T);
    c = x(1:end-2*T).*x(1+2*T:end);
    d = y(1:end-2*T).*y(1+2*T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [zeros(1,T) instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{14}(t)=\frac{2}{\pi T}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{(x(t)+x(t+T))^2}_{\text{c}}+\underbrace{(y(t)+y(t+T))^2}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = formula14(analytic_signal, Fs, T);
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = 2*Fs/(pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = (x(1:end-T)+x(1+T:end)).^2;
    d = (y(1:end-T)+y(1+T:end)).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

We wszystkich 3 przybliżeniach wzory T zostały ustawione na Fs (T = Fs = 1000 = 1s), jak sugerowano w pracy.

Teraz moje pytania to:

• Formuły f2 i f3 zwracają ten sam wynik dla czystego sygnału ćwierkającego. Myślę, że to dobrze, ponieważ obliczają to samo. Trzy metody aproksymacji nie zwracają tego samego, nawet czegoś, co jest blisko! Dlaczego tak jest? (Mam nadzieję, że to nie tylko błąd programowy ...)
• Chociaż wrócą same, zwłaszcza na końcu działki zaczną „Wiggle” dużo . Jakie jest tego wytłumaczenie? Najpierw pomyślałem o aliasingu, ale moja częstotliwość próbkowania jest dość wysoka w porównaniu z częstotliwością sygnałów, więc myślę, że można to wykluczyć.

• Wydaje się, że przynajmniej f2 i f3 działają poprawnie na czysty sygnał ćwierkający, ale wszystkie metody, w tym f2 i f3, wydają się okropnie zawodzić, jeśli chodzi o więcej niż jedną częstotliwość w sygnale. W rzeczywistości zawsze występuje więcej niż jedna częstotliwość w sygnale. Jak więc uzyskać (mniej więcej) prawidłową częstotliwość chwilową?

  • Właściwie nawet nie wiem, czego się spodziewać, gdy w sygnale jest więcej niż jedna częstotliwość. Obliczenia zwracają jedną liczbę dla danego punktu w czasie, więc co należy zrobić, gdy, tak jak tutaj, jest więcej częstotliwości? Zwraca średnią wszystkich częstotliwości czy coś takiego?

• I moim prawdopodobnie najważniejszym pytaniem jest, jak sobie z tym poradzić w prawdziwym i rozbudowanym oprogramowaniu? Powiedzmy, że chcę znać chwilową częstotliwość modulowanego sygnału na 1,75 s, i wybrałem metodę f2, niż mogę być „szczęśliwy” i uzyskać liczbę bliską 6 [Hz], co jest najprawdopodobniej poprawną odpowiedzią, lub ja wybierz moje wyniki obok kilku próbek i nagle dostaję trochę przewodowego wyniku, zbyt wysokiego, ponieważ niestety wybrałem wartość w piku. Jak sobie z tym poradzić? Czy przetwarzając go ze średnim, a nawet lepszym filtrem mediany? Myślę, że nawet to może stać się naprawdę trudne, szczególnie w regionach, w których wiele skoków znajduje się obok siebie.

I ostatnie, nie tak ważne pytanie, dlaczego większość dokumentów, które znajduję na temat częstotliwości chwilowych, pochodzi z obszaru geografii, szczególnie w obliczaniu zdarzeń sejsmograficznych, takich jak trzęsienia ziemi. Papier Barne'a również bierze to za przykład. Czy częstotliwość chwilowa nie jest interesująca w wielu obszarach?

Na razie tyle, jestem bardzo wdzięczny za każdą odpowiedź, szczególnie gdy ktoś daje mi wskazówki, jak zaimplementować ją w prawdziwym projekcie oprogramowania ;)

Z pozdrowieniami, Patrick

Naprawdę nie odpowiedź, ale może pomocna: osobiście przekonałem się, że koncepcja częstotliwości chwilowej jest użyteczna tylko w przypadku wystarczająco wąskich pasm.

Rozważ prosty przykład dwóch stałych fal sinusoidalnych, powiedzmy 100 Hz i 934 Hz. W takim przypadku możesz z pewnością zdefiniować i obliczyć częstotliwość chwilową (w dowolny sposób), ale jaki powinien być wynik? Jaki możliwy wgląd lub właściwości może mieć chwilowa częstotliwość, która mówi coś znaczącego na temat sygnału? Zastosowanie pojęcia częstotliwości chwilowej do sygnałów o wielu częstotliwościach w tym samym czasie po prostu nie ma większego sensu.

Właśnie dlatego uzyskujesz przyzwoite wyniki dla Sweepów, ale dziwne krzywe dla Sweep + sinus. Jest to również powód, dla którego widzisz wiggles wysoką część wymiata. Przepustowość sygnału staje się zbyt wysoka, aby przypisać mu pojedynczy numer częstotliwości, więc wyniki podskakują.

Wydaje się, że przynajmniej f2 i f3 działają poprawnie na czysty sygnał ćwierkający, ale wszystkie metody, w tym f2 i f3, wydają się okropnie zawodzić, jeśli chodzi o więcej niż jedną częstotliwość w sygnale. W rzeczywistości zawsze występuje więcej niż jedna częstotliwość w sygnale. Jak więc uzyskać (mniej więcej) prawidłową częstotliwość chwilową?

jak sugeruje Hilmar, metoda transformacji Hilberta (lub „sygnału analitycznego”) nie działa na szerokopasmowym, ponieważ istnieje więcej niż jeden składnik częstotliwości. możesz wykonać tę metodę tylko dla jednego komponentu sinusoidalnego.

więc dzięki podejściu Analitycznego Sygnału chcesz użyć tej tożsamości:

$$\arctan u - \arctan v = \arctan \left( \frac{u-v}{1+uv} \right)$$

gdyby $|u-v|$ jest wystarczająco mały, co można uzyskać Barnera "$f_9$„formuła z.

ale w obliczeniach transformacji Hilberta musi być tylko jedna sinusoida zmieniająca się w czasie, aby wykonać to poprawnie. i lepiej dopasuj komponent „w fazie” do wyjścia transformacji Hilberta (która jest opóźniona przez przyczynowy filtr FIR). inaczej dostaniesz bzdury.

Wow, co za wielkie pytanie. Najpierw odpowiem na nie tak ważne pytanie:

I ostatnie, nie tak ważne pytanie, dlaczego większość dokumentów, które znajduję na temat częstotliwości chwilowych, pochodzi z obszaru geografii, szczególnie w obliczaniu zdarzeń sejsmograficznych, takich jak trzęsienia ziemi. Papier Barne'a również bierze to za przykład. Czy częstotliwość chwilowa nie jest interesująca w wielu obszarach?

Powodem jest to, że system sejsmograficzny „vibroseis” jest wykorzystywany w przemyśle naftowym do wykonywania badań sejsmicznych. Ciężarówki, które podłączyłem, wibrują od około 5 Hz do około 90 Hz i mogą być wytwarzane do generowania sygnałów ćwierkających. W przemyśle naftowym jest dużo pieniędzy, a przetwarzanie zwrotów z tych sygnałów może być bardzo, bardzo lukratywne. Dlatego wiele osób spędzało wiele godzin analizując takie sygnały, w tym patrząc na techniki częstotliwości chwilowych.

---

Co do waszych ważniejszych pytań: generalnie robienie różnic arytmetycznych i obliczanie arcus tangentów na sygnałach w czasie dyskretnym jest złą rzeczą$^{\rm TM}$. Wynika to z faktu, że estymaty częstotliwości z czasem dyskretnym należy obliczać przy użyciu „arytmetyki kołowej” (arytmetyki wektorów AKA).

Sprawdź ten artykuł.

Lepsze podejścia zwykle wykorzystują „uśrednione fazy ważone”, jak tu wdrożono . Lub tutaj, aby uzyskać bezpośredni link do Matlaba .

Przepraszam, że odpowiedziałem rok po fakcie, ale natknąłem się na ten post, szukając artykułów na ten właśnie temat. Twoje pytania odzwierciedlają powszechne nieporozumienia i interpretacje „częstotliwości chwilowej”, które nękały to pole od samego początku. Wiele osób powie, jak niektóre odpowiedzi tutaj, że IF ma zastosowanie tylko do sygnałów „wąskopasmowych” lub „mono-komponentowych”. W rzeczywistości nie jest to prawdą: czasami IF uzyskany przez transformację Hilberta doskonale zachowuje się w przypadku sygnałów szerokopasmowych i / lub sygnałów „wieloskładnikowych”. Jedną z proponowanych wielkości, która pozwala uniknąć wielu z tych trudności, jest „średnia ważona częstotliwość chwilowa (WAIF)”, którą można zmierzyć za pomocą spektrogramu.

Zobacz Loughlin w J. Acoust. Soc. Popr., Styczeń 1999. Inne dobre prace na temat IF i powszechnych nieporozumień są autorstwa Picinbono (IEEE Trans. Sig. Proc., Marzec 1997) i Vakman (IEEE Trans. Sig. Proc., Kwiecień 1996).