4. Analiza falkowa

Z zasady działania transformaty Fouriera wynika jasno, że rozdzielczość czasowa analizy dokonanej z jej pomocą, nie jest zadowalająca. Innymi słowy, jeżeli analizujemy sygnał, który zmienia się gwałtownie z poziomu –1V do +1V to do jego opisu musimy użyć bardzo dużej liczby współczynników (teoretycznie – nieskończenie wiele), a ponadto nie jesteśmy w stanie określić chwili czasu, w której ta zmiana nastąpiła. Dla jej określenia możemy oczywiście użyć opisanego algorytmu STFT. Jednakże, znaczącą poprawę jakości analizy otrzymamy dopiero po „umiejętnym” zastosowaniu metody falkowej.

W przypadku transformaty falkowej, najbardziej charakterystyczne jest to, że w odróżnieniu od funkcji sinus i cosinus (charakterystycznych dla transformaty Fouriera), indywidualne funkcje falkowe są dobrze zlokalizowane w czasie (lub przestrzeni – dla obrazów) i jednocześnie podobnie jak sinus i cosinus, indywidualne falki są dobrze zlokalizowane w częstotliwości, ściśle biorąc tzw. skali. Ponadto w odróżnieniu od funkcji sinus i cosinus, które definiują unikalną transformatę Fouriera, nie ma pojedynczego, unikalnego zbioru falkowych funkcji bazowych. Falki różnią się między sobą zwartością lokalizacji czasowej oraz płynnością i gładkością kształtów. Wynikająca stąd zdolność falek do opisu sygnałów „z nieciągłościami”, przy ograniczonej liczbie współczynników oraz z lokalizacją w czasie, stanowi o jej przewadze nad transformatą Fouriera. Zaryzykować można stwierdzenie, że analiza falkowa jest „lokomotywą” w dziedzinie cyfrowego przetwarzania sygnałów.

Na rys. 4.1 zilustrowano ideę czterech znanych metod analizy sygnałów: czasową, częstotliwościową, czasowo-częstotliwościową (STFT) oraz falkową. Widać na nim wyraźnie, że w odróżnieniu od metody STFT, gdzie rozdzielczość czasowo-częstotliwościowa jest ustalona na całej płaszczyźnie t/f, w metodzie falkowej rozmiary okna czasowo-częstotliwościowego są funkcją jego położenia na tej płaszczyźnie.

Analiza czasowa

 

Analiza częstotliwościowa

 

Analiza STFT

Analiza falkowa

Rys. 4.1. Porównanie metod analizy: czasowa (obserwacja właściwości sygnału w dziedzinie czasu),  częstotliwościowa (obserwacja właściwości sygnału w dziedzinie częstotliwości), STFT (obserwacja właściwości sygnału na płaszczyźnie czas-częstotliwość), falkowa (obserwacja właściwości sygnału na płaszczyźnie czas-skala)

 

 

Swoją niezwykłą efektywność, a zarazem popularność w zakresie analizy sygnałów, transformata falkowa zawdzięcza szybkiemu algorytmowi, opracowanemu przez Mallata w roku 1989, zwanemu piramidą Mallata [13]. Algorytm ten wykorzystywany jest do uzyskania dekompozycji sygnału na składowe falkowe z użyciem tzw. kwadraturowych filtrów lustrzanych. Zastosowane w nim podejście wielorozdzielcze przenosi metodę falkową w realia kodowania podpasmowego.

Chociaż zakres pozycji literaturowych poświęconych metodzie falkowej jest niesłychanie rozległy [7, 10,13], to ogromna większość spośród nich pisana jest przez specjalistów dla specjalistów. Istnieją pozycje polskojęzyczne [2, 20]. Na szczególną uwagę zasługuje pozycja [22]. Celem tego rozdziału jest przybliżenie czytelnikowi praktycznych aspektów analizy falkowej, co miejmy nadzieję, będzie pomocne w zrozumieniu idei oraz prowadzeniu eksperymentów w tej dziedzinie [16].

Ciągłą (całkową) transformatę falkową (Continuous Wavelet Transform: CWT) funkcji x(t)ÎL2 dla pewnej falki y(t) definiuje się jako:

 

  W\left(\tau,\sigma\right)=:\int_{-\infty}^{+\infty}{x\left(t\right)\overline{\psi_{\tau,\sigma}\left(t\right)}}dt\begin{matrix},&\psi_{\tau,\sigma}\left(t\right)=\frac{1}{\sqrt\sigma}\psi\frac{t-\tau}{\sigma}\\\end{matrix} (4.1)

 

We wzorze tym parametr s oznacza skalę, zaś t – przesunięcie, co odpowiada ich funkcjom pełnionym w zapisie wzoru falkowego. Po to, aby funkcja y(t) mogła stanowić funkcję okna, a także być wykorzystana do odtworzenia x(t) musi spełniać warunek:

 

  \Psi\left(\omega\right)\left|\omega=0\right.=\int_{-\infty}^{+\infty}\psi\left(t\right)dt=:0 (4.2)

 

Ściśle biorąc, ciągła transformata falkowa osadzona jest więc w przestrzeni czas-skala (t/s), a nie czas-częstotliwość (t/f). Jednakże, w konkretnym przypadku, po dokonaniu odpowiedniej transformacji można przeliczyć skalę na częstotliwość. Miarą częstotliwości jest 1/s. Parametr t symbolizuje lokalizację okna wzdłuż osi czasu. Wyrażenie (4.2) natomiast wymusza oscylacje falek, które nie były wymagane w stosunku do okna STFT. Jednakże definitywnie yt,s(t) zajmuje miejsce jt,x(t) i zachowuje się tak jak funkcja okna. Natomiast widmo tego okna Y(w=0) = 0 i odwzorowuje filtr pasmowy.

Warto przy tym zwrócić uwagę, że iloczyn promieni okien, czasowego i częstotliwościowego, ma wartość stałą na całej płaszczyźnie:

 

  2\sigma\Delta_t\cdot\frac{2}{\sigma}\Delta_\omega=4\Delta_t\Delta_\omega (4.3)

 

Położenie okna czasowo-częstotliwościowego na płaszczyźnie t/f, dla transformaty falkowej, pokazano na rys.4.2.

 

Rys. 4.2. Zobrazowanie idei transformaty falkowej na płaszczyźnie t/f

Podobnie jak w przypadkach DFT oraz STFT definiuje się pojęcie dyskretnej transformaty falkowej (Discrete Wavelet Transform: DWT). W tym celu przyjmuje się, że:

 

  \sigma=2^{-s}\begin{matrix},&\tau=\\\end{matrix}2^{-s}\cdot\ l (4.4)

 

co w konsekwencji, po dodatkowym przyjęciu założenia o dyskretyzacji x(t), daje nową formę zapisu transformaty:

 

 

W\left(l,s\right)=W\left(l2^{-s},2^{-s}\right)=W\left(l,s\right)=

=2^{s/2}\int_{-\infty}^{+\infty}{x\left(t\right)\psi\left(2^st-l\right)dt}=

=2^{s/2}\sum_{n}\ x\left(n\right)\psi\left(2^sn-l\right)

(4.5)

 

Warto pamiętać, że transformata falkowa nie spełnia warunku niezmienności względem przesunięcia, a przesunięcie funkcji w czasie: xm(t)=x(t-tm), objawia się w dziedzinie falkowej raczej w sposób bardzo nieprzyjazny [10]:

 

 

W_{\psi,x}\left(l,s\right)=2^{s/2}\int_{-\infty}^{+\infty}{x_m\left(t\right)\psi\left(2^st-l\right)dt=}

=2^{s/2}\sum_{n}\ x\left(n-m\right)\psi\left(2^sn-l\right)=

=2^{s/2}\sum_{n}{x\left(n\right)\psi\left[2^sn-\left(l-m2^s\right)\right]=}

=W_{\psi,x}\left[\left(l-m2^s\right)2^{-s},2^{-s}\right]

(4.6)

 

Na tym etapie rozważań można się pokusić o sformułowanie wyrażenia na szereg falkowy, który istnieje dla dowolnej funkcji x(t)ÎL2:

 

  x\left(t\right)=\sum_{s}\sum_{l}{w_{l,s}\psi_{l,s}\left(t\right)\begin{matrix},&\psi_{l,s}\left(t\right)\\\end{matrix}=2^{s/2}\psi\left(2^st-l\right)} (4.7)

 

Jeżeli zatem {<i><span style="font-family:Symbol">y</span></i><sub>k,s</sub>(<i>t</i>)} tworzy ortonormalną bazę w przestrzeni L2(Â) to podobnie jak w przypadku szeregu Fouriera:

 

  w_{l,s}=\left\langle x\left(t\right),\psi_{k,s}\left(t\right)\right\rangle=2^{s/2}W_{\psi,x}\left(l,s\right) (4.8)

 

Okazuje się, że w odróżnieniu od transformaty Fouriera szereg falkowy otrzymuje się po spróbkowaniu transformaty ciągłej na płaszczyźnie t/f, w wybranych punktach diadycznych: (l2-s, 2-s). Analiza jest utrudniona, gdy baza {<i><span style="font-family:Symbol">y</span><sub>k,s</sub></i>(<i>t</i>)} nie spełnia warunku ortonormalności. Wymaga to zdefiniowania tzw. falki dualnej i sformułowania bazy: ψl,st  [10, 13].

W tabeli 4.1 zamieszczono przykłady skalowania funkcji dla przypadku sinusoidy oraz pewnej (typowej) falki (a - współczynnik skali). W przypadku sinusoidy istnieje ścisłe odwzorowanie skala-częstotliwość, dla falki nie jest ono tak oczywiste. W związku z tym, jak już wspomniano wcześniej, pozostaje się przy pojęciu skalowania.

 

Tabela 4.1  Przykłady skalowania funkcji

x\left(t\right)=sin{t}

a=1

Wysoka skala

x\left(t\right)=sin{2}t

a=1

x\left(t\right)=sin{2}t

a=2

Średnia skala

x\left(t\right)=\psi\left(2t\right)

a=2

x\left(t\right)=sin{4}t

a=4

Nska skala

x\left(t\right)=\psi\left(4t\right)

a=4

       

 

 

Proces skalowania falki może przebiegać w dwu kierunkach, określa się je mianem kompresji (ściskania) lub rozciągania. W tabeli 4.1 do skalowania falki zastosowano kompresję.

Drugi parametr rozkładu falkowego to przesunięcie. Ilustracja przesuwania w czasie zawarta jest w tabeli 4.2.

 

Tabela 4.2 Ilustracja przesuwania w czasie falki

 

\psi\left(t\right) \psi\left(t-\tau\right) \psi\left(t-2\tau\right)

 

Na tle powyższych rozważań można, w sposób opisowy, zdefiniować proces rozkładu falkowego. Zawiera on 5 charakterystycznych kroków (tabela 4.3).

 

Tabela 4.3 Ilustracja rozkładu falkowego

 

1. Wybraną falkę ustawić na początku fragmentu sygnału przeznaczonego do analizy,Wyznaczyć umowną wartość liczbową odpowiadającą korelacji między bieżącą falką i odpowiadającym jej segmentem sygnału.

Uwaga! – w przypadku unormowania energii sygnału w aspekcie użytej falki, wspomniana liczba będzie równoważna wartości współczynnika korelacji wzajemnej między falką, a wybranym segmentem sygnału.

2 .Przesunąć falkę o jeden cykl w prawo i powtórzyć działanie opisane w kroku 2. Sekwencję kroków 3,2 powtarzać aż do końca trwania sygnału.

3. Rozciągnąć falkę i powtórzyć kroki od 1 do 3.

4. Powtórzyć kroki od 1 do 4 aż do wyczerpania wszystkich skal.

 

Powyższy przykład operuje w obrębie tzw. diadycznego charakteru zmian w obrębie skali i przesunięcia charakterystycznych dla dyskretnej transformaty falkowej (DWT). Pod pojęciem ciągłej transformaty falkowej (CWT) kryje się sposób umożliwiający użycie dowolnej, zmienianej w sposób ciągły skali oraz ciągłego przesunięcia w czasie. Oczywiście, w kontekście sygnałów dyskretnych ciągłość, w obydwu wskazanych przypadkach, oznacza zmiany w obrębie jednej próbki sygnału (co jedną próbkę),

Należy zauważyć, że w analizie falkowej wyższa skala równoważna jest bardziej rozciągniętej falce (patrz tabela 4.1). Im bardziej rozciągnięta falka (wyższa skala) tym większa sekcja sygnału, z którą jest porównywana i tym bardziej zgrubne cechy sygnału wyeksponowane za pomącą odpowiadającego jej współczynnika. Podsumowanie tego spostrzeżenia zawarto w tabeli 4.4

Tabela 4.4  Charakterystyka analizy falkowej

Niska skala analizy falkowej

Ściśnięta falka

Szybkozmienne
detale

Wysoka częstotliwość analizy

Wysoka skala analizy falkowej

Rozciągnięta falka

Wolnozmienne
cechy sygnału

Niska częstotliwość analizy

 

Fakt, że analiza falkowa nie obrazuje cech sygnału na płaszczyźnie czas-częstotliwość, lecz czas-skala nie stanowi o słabości metody, a wręcz przeciwnie – o jej sile. Okazuje się, że jest to naturalna metoda opisu wielu zjawisk fizycznych odbieranych przez zmysły człowieka. Trzeba się z nią pogodzić i do niej przyzwyczaić.

Zjawiskiem ściśle związanym z metodą falkową jest analiza wielorozdzielcza [1]. Jak wskazuje nazwa, umożliwia ona obserwację właściwości sygnału na różnych poziomach aproksymacji lub inaczej rozdzielczości. Stosując analizę wielorozdzielczą dokonuje się podziału złożonej funkcji na kilka (lub więcej) prostych składowych i bada je niezależnie.

Bardzo efektywna metoda implementacji tego algorytmu analizy dokonanej z użyciem filtrów opracowana została w 1988 roku przez Mallata [13]. Nawiązuje ona do, znanej z analizy częstotliwościowej, metody kodowania w podpasmach i realizuje tzw. szybką transformatę falkową (Fast Wavelet Transform: FWT).

Do analizy falkowej wprowadzono dwa bardzo adekwatne pojęcia: aproksymacji i detalu. Pod pojęciem aproksymacji rozumie się wysoką skalę analizy, niskoczęstotliwościowe składowe sygnału. Detale to niska skala i wysokoczęstotliwościowe składowe. Wspomniany proces filtracji, w swoim najniższym poziomie obejmuje dwa filtry: dolnopasmowy (L) i górnopasmowy (H). Filtr dolnopasmowy odtwarza aproksymację, a górnopasmowy detal sygnału (rys.4.3).

 

Rys. 4.3. Schemat filtracyjnej dekompozycji falkowej

Oryginalny sygnał S przechodzi przez parę komplementarnych filtrów, które rozdzielają go na dwie składowe A (Aproksymacja) i D (Detal). W przypadku filtracji cyfrowej podwaja się liczba danych, przeznaczonych do dalszego przetwarzania. Wygodnym sposobem ograniczenia tej liczby w metodzie falkowej jest decymacja, polegająca na odrzuceniu co drugiej próbki danych. Przykład działania filtrów falkowych wraz z interpretacją zawiera tabela 4.5.

Tabela 4.5  Interpretacja działania filtrów falkowych

x(n)

Sygnał oryginalny: sinusoida zaszumiona szumem białym o rozkładzie normalnym amplitud.

 

 

 

xa(n)

Odfiltrowana za pomocą falkowego filtru dolnopasmowego aproksymacja sygnału oryginalnego, w postaci sinusoidy o znacznie mniejszym poziomie szumu.

xd(n)

Odfiltrowane za pomocą falkowego filtru górnopasmowego detale sygnału oryginalnego, w postaci szumu wysokoczęstotli-wościowego.

 

Przykład pierwszego poziomu dekompozycji falkowej pewnego rzeczywistego sygnału pomiarowego zamieszczono na rysunku 4.4.

 

Rys. 4.4. Przykład pierwszego poziomu dekompozycji falkowej sygnału

Pełny proces dekompozycji zawiera szereg członów tworzących tzw. drzewo dekompozycji falkowej. Przykład takiego drzewa zamieszczono na rysunku 4.5.

Rys. 4.5. Drzewo dekompozycji falkowej: S – sygnał oryginalny,
 di – detal w i-tej skali, ai –aproksymacja w i-tej skali

Z powodu tego, że proces dekompozycji ma charakter iteracyjny, teoretycznie może on trwać do nieskończoności. W rzeczywistości, może być kontynuowany tylko do chwili, aż w kolejnym detalu pozostanie pojedyncza próbka. W praktyce zaś, wybiera się odpowiednią liczbę poziomów, dostosowaną do natury sygnału lub innego kryterium.

Z metrologicznego punktu widzenia można by poprzestać na omówieniu procesu dekompozycji falkowej, zwanego inaczej procesem analizy. Nierzadko jednak, np. w zastosowaniu do kompresji sygnałów konieczne jest, możliwie wierne, odtworzenie postaci sygnału oryginalnego. Proces ten, w metodzie falkowej określany jest mianem rekonstrukcji lub syntezy. Matematyczny zapis tego procesu obejmuje definicje odwrotnej dyskretnej transformaty falkowej. Metoda filtracyjna zaś, polega na użyciu tzw. filtrów rekonstrukcyjnych. Znamienna dla metody falkowej jest konieczność tzw. rekonstrukcji perfekcyjnej. Jest ona następstwem systematycznego wprowadzania zniekształceń składowych sygnału występujących w procesie analizy drogą usuwania co drugiej próbki. W procesie rekonstrukcji, brakujące próbki zastępowane są wartościami zerowymi. Filtry rekonstrukcyjne zaś dobrane są w taki sposób, aby „niwelowały” zgubne następstwa aliasingu. Wspólnie z filtrami użytymi w procesie analizy tworzą one system określany mianem kwadraturowych filtrów lustrzanych: H-H’,L-L’ (rys.4.6).

 

Rys. 4.6. Złożenie dekompozycji i rekonstrukcji

Dobór charakterystyk filtrów rozkładu falkowego jest jednoznaczny z doborem kształtu falki tego rozkładu. Wydaje się, że z punktu widzenia analizy metrologicznej ten sposób doboru parametrów jest lepszy niż pierwotne skonstruowanie falki, a potem projektowanie dla niej filtru[1]. Zwłaszcza, że niemożliwe jest użycie całkiem dowolnej ograniczonej w czasie funkcji, o zerowej wartości średniej, nazwanie jej falką i użycie w procesie analizy, jeżeli zamierza się zrekonstruować sygnał w przyszłości. Trzeba to zrobić w kontekście istnienia zespołu kwadraturowych filtrów lustrzanych.

Ściśle rzecz ujmując, kształt falki c(t) jest jednoznacznie związany z charakterystyką filtru górnopasmowego wyodrębniającego detal w rozkładzie falkowym. Istnieje jeszcze jedna bardzo charakterystyczna funkcja, która związana jest z niektórymi (nie wszystkimi) zbiorami falek. Jest to tzw. funkcja skalująca, oznaczana symbolem w(t). Jej kształt związany jest z charakterystykami dolnopasmowych kwadraturowych filtrów lustrzanych odpowiedzialnych za wyodrębnienie aproksymacji sygnału. Kształt funkcji skalującej jest zbliżony do kształtu odpowiadającej jej falki. Definiuje się ją w rekurencyjnym zapisie matematycznym za pomocą równania dylatacyjnego:

 

  \phi\left(t\right)=\sqrt2\sum_{k=0}^{N-1}{h_k\phi\left(2t-k\right)} (4.9)

 

W kontekście funkcji skalującej falka zdefiniowana jest za pomocą tego samego równania dylatacyjnego opisanego za pomocą innego zestawu współczynników.

 

  \psi\left(t\right)=\sqrt2\sum_{k=0}^{N-1}g_k\phi\left(2t-k\right) (4.10)

 

Współczynniki H={<i>h<sub>k</sub></i>}, oraz G={<i>g<sub>k</sub></i>} są rozumiane jako współczynniki pary kwadraturowych filtrów lustrzanych. W przypadku bazy ortonormalnej związane są zależnością wzajemną: gk=(-1)khN-k.

Jak już wspomniano wcześniej, falkowe funkcje bazowe otrzymuje się z falki matki drogą skalowania i przesunięcia:

 

  \psi_{s,l}\left(t\right)=\sqrt{2^s}\psi\left(2^st-l\right) (4.11)

 

Istnieje nieograniczona wręcz liczba możliwych do utworzenia falek jak i tzw. banków filtrów. Która z nich jest najlepsza zależy od konkretnej implementacji, tak jak było to w przypadku doboru okna czasowego w analizie widmowej. Nazwy konkretnych rozwiązań pochodzą zwykle od kształtów lub nazwisk osób, które je po raz pierwszy użyły i opublikowały wyniki. Przykładowe nazwy to: Daubechies, Haar, Coiflets, Symlet, Spline, Battle-Lemarie.

W zasadzie każda dyskusja o falkach może rozpoczynać się od funkcji Haara. Ta falka równoważna jest Daubechies (db1). Ingrid Daubechies wprowadziła rodzinę falek, które są często określane mianem zbioru ortonormalnego o zwartym nośniku [7].

Dla ustanowienia pełnego obrazu metody falkowej poniżej (rys.4.7 i 4.8) przedstawiono wyniki analizy falkowej (CWT) uzyskane z wykorzystaniem wirtualnego przyrządu pomiarowego [11][2]. Znaczenie wykresów, idąc od dołu do góry, jest następujące: przebieg czasowy sygnału, zbiór współczynników w wybranej skali - zaznaczonej kursorem na wykresie górnym, rozkład falkowy na płaszczyźnie czas-skala (t/s). Obydwa prezentowane sygnały zostały wygenerowane w sposób sztuczny i były użyte poprzednio do zobrazowania właściwości STFT (rys.4.7 i 4.8). Konfrontacja obydwu par rysunków jest bardzo interesująca. W szczególności dotyczy to wzajemnej relacji między częstotliwością i skalą.

 

Rys. 4.7. Przykładowy wynik analizy falkowej (CWT): od dołu: przebieg czasowy sygnału, zbiór współczynników w wybranej skali - zaznaczonej kursorem, rozkład falkowy na płaszczyźnie czas-skala (t/s)

 

Rys. 4.8. Przykładowy wynik analizy falkowej (CWT): od dołu: przebieg czasowy sygnału, zbiór współczynników w wybranej skali - zaznaczonej kursorem, rozkład falkowy na płaszczyźnie czas-skala (t/s)

Na zakończenie opisu metody falkowej trzeba wspomnieć o istnieniu uogólnionej metody rozkładu falkowego zwanej pakietami falkowymi, dostarczającej bogatszego zbioru składowych sygnału. W odróżnieniu od standardu w tej metodzie podlegają dalszemu rozkładowi również detale. Schemat drzewa dekompozycji pakietowej pokazano na rysunku 4.9.

Rys. 4.9. Drzewo falkowej dekompozycji pakietowej: S – sygnał oryginalny, d – detal, a – aproksymacja, dd – detal (detalu) w kolejnym rozkładzie, da – aproksymacja (detalu) w kolejnym rozkładzie, ad – detal (aproksymacji) w kolejnym rozkładzie, aa – aproksymacja (aproksymacji) w kolejnym rozkładzie

 

 

 

[1] W szczególnych przypadkach może się okazać, że łatwiej jest dobrać kształt falki do kształtu sygnału. Przykładem może tu być sygnał elektrokardiogramu.

[2] Przyrządy opracowane w Zakładzie Miernictwa Elektrycznego IETiME PW w ramach projektu bad.