Podręcznik: Krótkoczasowa transformata Fouriera

3. Analiza czasowo-częstotliwościowa

3.1. Krótkoczasowa transformata Fouriera

Jeżeli mamy do czynienia z sygnałem nieokresowym, sumowanie funkcji okresowych typu sinus i cosinus, nie może dać wiarygodnej reprezentacji sygnału. Sygnał ten można w sposób sztuczny doprowadzić do okresowości, ale wymaga to zadbania o zachowanie ciągłości na krańcach przedziałów. Tak zwana Krótkoczasowa Transformata Fouriera (Short Time Fourier Transform – STFT) jest skutecznym narzędziem rozwiązania problemu budowania reprezentacji widmowej funkcji nieokresowej.

Krótkoczasowa Transformata Fouriera umożliwia wydobycie z sygnału informacji o tym, jak zmienia się jego widmo w czasie, czyli jednoczesną obserwację jego właściwości zarówno w dziedzinie czasu jak i częstotliwości. Wycinek sygnału (blok próbek o rozmiarze L) przeznaczony do analizy zostaje podzielony na segmenty, a każdy segment podlega analizie widmowej niezależnie. Podobnie jak w przypadku tradycyjnym, aby usunąć gwałtowne zmiany (cięcia) sygnału na krańcach przedziałów, stosuje się wymyślne okna czasowe w odniesieniu do wspomnianych sekcji. Sukcesywne „przesuwanie okna czasowego” umożliwia lokalizację parametrów widmowych sygnału w czasie.

Ten proces daje się zapisać za pomocą prostego równania:

$X\left(t,\omega\right)=\int_{-\infty}^{\infty}{x\left(t\right)\phi\left(t-\tau\right)e^{-j\omega\tau}d\tau}$

(3.1)

w którym j () opisuje funkcję okna czasowego. Przesuwając okno j () w czasie, wzdłuż sygnału, wyznacza się jego zawartość widmową wewnątrz przedziału czasowego, którego długość jest określona szerokością okna. Kształt okna czasowego odgrywa kluczową rolę w przypadku STFT. Iloczyn szerokości okna w dziedzinie czasu i szerokości okna w dziedzinie częstotliwości jest wielkością stałą dla danego okna. Stąd też, poprawiając rozdzielczość w dziedzinie czasu, będziemy ją pogarszać w dziedzinie częstotliwości i odwrotnie. Zatem szerokość okna wybierana jest na drodze kompromisu.

Przechodząc do konkretów, załóżmy, że j (t)ÎL² (Â), jest funkcją rzeczywistą i opisuje kształt okna czasowego. Efekt zastosowania okna czasowego opisuje zależność:

$x_\tau\left(t\right)=:x\left(t\right)\phi\left(t-\tau\right)$

(3.2)

Dla przypadku okna prostokątnego (rys.3.1): $\phi_r\left(t\right)=\chi_{\left[-\varepsilon,\varepsilon\right)}\left(t\right)$

Rys. 3.1. Prostokątne okno czasowe

Efekt działania prostokątnego okna czasowego opisuje zależność (3.3).

$x_\tau\left(t\right)=\left\{\begin{matrix}x\left(t\right),\begin{matrix}&t\in\left[\tau-\varepsilon,\tau+\varepsilon\right)\\\end{matrix}\\0\begin{matrix}&\\\end{matrix}\mathrm{dla}\ \buildrel\over{}\mathrm{pozostalych}\ \\\end{matrix}\right.$

(3.3)

Teoretycznie, funkcja nie może być jednocześnie ograniczona w dziedzinie czasu i w dziedzinie częstotliwości. Jednak, przy tak zdefiniowanych promieniach okien, można osiągnąć jednoczesną ich ograniczoność w dziedzinie czasu i częstotliwości. W takim przypadku j(t) określa się mianem okna czasowo-częstotliwościowego. Niestety, okno prostokątne nie jest tutaj dobrym przykładem, gdyż D_j= 3^-1/2, a D_F= ¥. Wynika stąd jasno, że jest to idealne okno czasowe i jednocześnie najgorsze, wręcz nieakceptowalne okno częstotliwościowe.

Miarą jakości okna czasowo-częstotliwościowego jest iloczyn promienia w dziedzinie czasu i częstotliwości, czyli „pole” okna położonego na płaszczyźnie t/f. Przy czym zgodnie z zasadą nieoznaczoności:

$\Delta_t\Delta_\omega\geq\frac{1}{2}$

(3.4)

W kontekście powyższych rozważań, krótkoczasową transformatę Fouriera (STFT) sygnału x(t), w odniesieniu do okna j(t) rozmieszczonego w pozycji (t,x) na płaszczyźnie t/f zdefiniować można jako:

$G\left(\tau,\xi\right)=:\int_{-\infty}^{+\infty}{x\left(t\right)\overline{\phi_{\tau,\xi}\left(t\right)}dt}$

(3.5)

gdzie: $\phi_{\tau,\xi}\left(t\right)=:\phi\left(t-\tau\right)e^{j\xi.t}$

Dopuszcza się, aby j(t) było funkcją zespoloną i w związku z tym musi spełniać uogólniony warunek zapewniający, że F(w) stanowi filtr dolnopasmowy:

$\Phi\left(0\right)=\int_{-\infty}^{+\infty}{\phi\left(t\right)dt\neq0}$

(3.6)

W odróżnieniu od tradycyjnej transformaty Fouriera, dla której, do wyznaczenia pojedynczej składowej, konieczna jest znajomość funkcji x(t) na całej osi czasu, w tym przypadku wymagana jest znajomość x(t) tylko w przedziale określonym przez położenie j(t-t).

Niezwykle interesujące jest podejście przyjmujące, że w STFT okno j(t) użyte jest do modulacji sygnału sinusoidalnego e^-j^xt. Wyrażenie opisujące STFT przyjmuje wtedy postać:

$G\left(\tau,\xi\right)=:\left\langle x\left(t\right)\phi\left(t-\tau\right)e^{j\xi.t}\right\rangle=\left\langle x\left(t\right)\phi_{\tau,\xi}\left(t\right)\right\rangle$

(3.7)

Można powiedzieć, że funkcja:

$\phi_{\tau,\xi}\left(t\right)=:\phi\left(t-\tau\right)e^{j\xi\begin{matrix}t&\\\end{matrix}}$

(3.8)

zachowuje się tak jak pakiet fal, gdzie sinusoida oscyluje wewnątrz okna j(t). Każdy pakiet traktować można jako pewną funkcję bazową. Podejście to pozwala traktować STFT jako rozkład x(t) na składowe należące do tak zdefiniowanej bazy na płaszczyźnie t/f.

Funkcję x_t(t) można odtworzyć wyznaczając odwrotną transformatę Fouriera o postaci:

$x_\tau\left(t\right)=\phi\left(t-\tau\right)x\left(t\right)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}{G\left(\tau,\xi\right)e^{j\xi t}d\xi}$

(3.9)

Funkcję oryginalną x(t), natomiast otrzyma się dopiero po przebiegnięciu wszystkich pozycji położenia okna:

$x( t) =\frac{1}{2\pi \Vert \phi ( t)\Vert }\int ^{+\infty }_{-\infty } d\xi \ e^{j\xi t}\int ^{+\infty }_{-\infty } G( \tau ,\xi )\overline{\phi ( t-\tau )} dt$

(3.10)

Transformata Gabora

Znany jest pewien przypadek szczególny dla STFT, w którym okno czasowe ma kształt gaussowski:

$\phi\left(t\right)=g_\alpha\left(t\right)=\frac{1}{2\pi\alpha}e^{-\frac{t^2}{4\alpha}},\begin{matrix}&\alpha>0\\\end{matrix}$

(3.11)

Okazuje się, że widmo tego okna ma również kształt gaussowski:

$\Phi\left(\omega\right)=G_\alpha\left(\omega\right)=e^{-\alpha\omega^2}\begin{matrix},&\alpha>0\\\end{matrix}$

(3.12)

Przypadek ten spełnia warunki okna czasowo-częstotliwościowego o parametrach:

okno czasowe: $t^*=0\begin{matrix},&\Delta_t=\sqrt\alpha\\\end{matrix}$
okno częstotliwościowe: $\omega^*=0\begin{matrix},&\Delta_\omega=\frac{1}{2\sqrt\alpha}\\\end{matrix}$

Okno to, od nazwiska pierwszego użytkownika, nosi miano okna Gabora [10]. Miara jakości okna Gabora wynosi:

$\Delta_t\cdot\Delta_\omega=\frac{1}{2}$

(3.13)

i osiąga dolną granicę zasady niepewności.

Idea analizy czasowo-częstotliwościowej z użyciem transformaty Gabora, przedstawiona jest na rys.3.2.

Rys. 3.2. Zobrazowanie idei STFT na płaszczyźnie t/f

Analityczny zapis tej transformaty przedstawia zależność (3.14):

$G\left(b,\xi\right)=\int_{-\infty}^{+\infty}{x\left(t\right)\overline{g_\alpha\left(t-\tau\right)} e^{-j\xi \ t}dt}=\frac{1}{2\pi\alpha}\int_{-\infty}^{+\infty}{f\left(t\right) e^\frac{-(t-\tau)^2}{4\alpha}e^{-j\xi \ t}dt}$

(3.14)

Właściwości STFT

Właściwości krótkoczasowej transformaty Fouriera ułatwiają zrozumienie istoty analizy czasowo-częstotliwościowej. Poniżej przedstawiono trzy najbardziej istotne:

liniowość: $x\left(t\right)=ax_1\left(t\right)+bx_2\left(t\right)$

$G_=aG_1\left(\tau,\xi\right)+bG_2\left(\tau,\xi\right)$

(3.15)

przesunięcie w czasie: $x_0\left(t\right)=x\left(t-t_0\right)$

$G_0\left(b,\xi\right)=\int_{-\infty}^{+\infty}{x\left(t-t_0\right)\phi\left(t-\tau\right)e^{-j\xi\cdot t}dt=}$

$=\int_{-\infty}^{+\infty}f\left(t\right)\phi\left[t-\left(\tau-t_0\right)\right]e^{-j\xi\cdot t}e^{-j\xi\cdot t_0}dt=$

$=e^{-j\xi\cdot t_0}G\left(\tau-t_0,\xi\right)$

(3.16)

Oznacza to, że jeżeli położenie oryginalnej funkcji x(t) zostanie przesunięte w czasie o t₀, to lokalizacja widma na płaszczyźnie t/f zostanie również przesunięta o t₀ w czasie, a lokalizacja częstotliwościowa pozostanie bez zmian. Dodatkowo pojawi się przesunięcie fazowe proporcjonalne do t₀.

przesunięcie w częstotliwości: $x_0\left(t\right)=x\left(t\right)e^{j\omega_0t}$

$G\left(\tau,\xi\right)=\int_{-\infty}^{+\infty}{x\left(t\right)e^{j\omega{_0}t}}\phi\left(t-\tau\right)e^{-j\xi t}dt=$

$=\int_{-\infty}^{+\infty}{x\left(t\right)\phi\left(t-\tau\right)e^{-j\left(\xi-\omega_0\right)t}}dt=$

$=\left\langle x\left(t\right),\phi\left(t-\tau\right)e^{-j\left(\xi-\omega_0\right)t}\right\rangle=G\left(\tau,\xi-\omega_0\right)$

(3.17)

Powyższa zależność potwierdza, że w przypadku „przesunięcia” funkcji x(t) w dziedzinie częstotliwości, widmo STFT zostaje również przesunięte na osi częstotliwości bez zmiany amplitudy i fazy.

Do celów obliczeń numerycznych, jak zwykle wprowadza się pojęcie Dyskretnej Transformaty Gabora, określanej w skończonym zbiorze punktów na płaszczyźnie t/f. Wtedy:

$G(\tau_n,\xi_k)=\int_{-\infty}^{+\infty}x(t) \overline {g_{\alpha}(t-\tau_n)}e^{-j\xi_kt}dt=$

$=\left\langle x\left(t\right),g_\alpha\left(t-\tau_n\right) \ e^{j\xi_kt}\right\rangle=\left\langle x\left(t\right),\phi_{n,k}\left(t\right)\right\rangle$

(3.18)

Funkcja x(t) może zostać odtworzona na podstawie zależności określanej mianem rozkładu Gabora:

$f\left(t\right)=\sum_{n}\sum_{k}{G\left(\tau_n,\xi_k\right)g_\alpha\left(t-\tau_n\right)e^{j\xi_kt}}$

(3.19)

Współczynnikami tego rozkładu są liczby G_n,k= G(t_n,x_k). W kontekście powyższego warto podać uogólnioną formułę rozkładu Gabora, w której występuje pakiet fal sinusoidalnych zmodulowanych za pomocą okna gaussowskiego o symbolu j_n,k(t):

$x\left(t\right)=\sum_{n}\sum_{k}{G_{n,k}\phi_{n,k}\left(t\right)}$

(3.20)

Pakiet ten określany jest mianem bazy rozkładu Gabora. Baza ta, z uwagi na gaussowski kształt okna o minimalnym rozmiarze na płaszczyźnie t/f, uchodzi za optymalną w sensie budowy reprezentacji sygnału.

Poniżej (rys.3.3) podany jest przykład zastosowania rozkładu STFT do sygnału o postaci [10]:

$x\left(t\right)=sin{2}\pi\ f_1t+sin{2}\pi\ f_2t+\alpha\left[\delta\left(t-t_1\right)+\delta\left(t-t_2\right)\right]$

(3.21)

dla różnych szerokości (promieni) D_t okna czasowego j(t).

W pierwszym przypadku, dla szerokiego okna czasowego, otrzymuje się bardzo dobrą rozdzielczość częstotliwościową i bardzo słabą rozdzielczość czasową. W miarę zwężania okna czasowego poprawia się rozdzielczość czasowa kosztem częstotliwościowej. Dla dobrego zobrazowania wszystkich elementów składowych sygnału, trzeba zmieniać parametry okna czasowo-częstotliwościowego i za każdym razem, od nowa wyznaczać STFT. Wymaga to dużych nakładów obliczeniowych.

Na zakończenie tego punktu na kolejnych rysunkach pokazano praktyczne wyniki analizy t/f, przeprowadzonej z użyciem transformaty Gabora, za pomocą wirtualnych przyrządów pomiarowych [11][1]. W obydwu przypadkach użyte zostały bardzo charakterystyczne sygnały, ułatwiające interpretację rozkładów widmowych [16]. W pierwszym przypadku (rys.3.4) jest to sygnał złożony z segmentów powstałych w wyniku skokowej zmiany częstotliwości, a w drugim (rys. 3.5) sygnał o zmodulowanej częstotliwości ze zmiennością wykładniczą.

Rys. 3.3. Przykład zastosowania rozkładu STFT dla czterech różnych szerokości okna czasowego

Rys. 3.4. Praktyczny wynik analizy t/f, przeprowadzonej z użyciem transformaty Gabora

Rys. 3.5. Praktyczny wynik analizy t/f, przeprowadzonej z użyciem transformaty Gabora

Podsumowanie właściwości aplikacyjnych STFT:

D_t i D_w - nie zależą od położenia na płaszczyźnie (t/f).
ustalona rozdzielczość czasowo – częstotliwościowa
W analizie czasowo-częstotliwościowej wskazane byłoby posługiwanie się oknem strojonym:

- szerokie w czasie (wąskie w częstotliwości) dla małych częstotliwości

- wąskie w czasie (szerokie w częstotliwości) dla dużych częstotliwości.

[1] Przyrządy opracowane w Zakładzie Miernictwa Elektrycznego IETiME PW w ramach projektu badawczego