Interpolacja

Transformacja Fouriera (ang. Fourier Transformation, FT) jest narzędziem analitycznym, które pozwala rozłożyć skomplikowane "obiekty" na proste elementy składowe. Jako przykład, wyobraź sobie, że w kuchni w naczyniu znajduje się pyszny sok. Gdy poddasz go analizie, to otrzymasz informację z czego ten sok się składa oraz w jakiej ilości zawiera swoje składniki: np. 25% borówek, 40% bananów i 35% mleka. Coś podobnego robi transformacja Fouriera z sygnałami: bierze skomplikowany sygnał, np. dźwięk z mikrofonu, i rozbija go na pojedyncze fale sinusoidalne, z których ten dźwięk jest złożony:

Mówiąc językiem fizyki i matematyki: FT to narzędzie, które bierze skomplikowany, „poszarpany” sygnał (np. dźwięk) i rozbija go na pojedyncze fale sinusoidalne, z których ten sygnał został złożony. Przenosi nas z dziedziny czasu \(t\) (gdzie widzimy tylko sumaryczny hałas) do dziedziny częstotliwości \(f\) (gdzie widzimy konkretne nuty/składniki):

W dziedzinie czasu \(t\) (na zwykłym wykresie fali dźwiękowej) sygnały są nieczytelne. Jeśli w tym samym momencie gra gitara, perkusja i ktoś krzyczy, mikrofon rejestruje jedną, chaotyczną, falującą linię.

FT rozwiązuje problem izolacji i analizy danych:

• Separacja: Pozwala komputerowi „usłyszeć” poszczególne instrumenty osobno z jednego nagrania.
• Kompresja: Zamiast zapisywać miliard punktów skomplikowanej fali, FT pozwala zapisać tylko informację typu: „gra dźwięk A przez 5 sekund i dźwięk B przez 3 sekundy”. To drastycznie zmniejsza wagę plików.
• Filtracja: Rozwiązuje problem szumu. Jeśli w nagraniu piszczy stary telewizor (wysoka częstotliwość), FT pokazuje ten pisk jako pojedynczy, wyraźny słupek na wykresie w dziedzinie częstotliwości \(f\). Wystarczy go „skasować” i cofnąć transformację, aby uzyskać czysty dźwięk.

Transformacja Fouriera zyskała na znaczeniu po opracowaniu szybkiego algorytmu dokonującego tej transformacji, tzw. FFT (ang. Fast Fourier Transformation). Tym zajmiemy się w dalszej części rozdziału. Transformację Fouriera będziemy oznaczać skrótem FT.

1805 rok – Carl Friedrich Gauss (Pierwszy, ukryty algorytm FFT)

Wielu uważa, że szybki algorytm obliczania transformacji Fouriera (FFT) powstał w latach 60. XX wieku. W rzeczywistości Carl Friedrich Gauss obliczył go jako pierwszy w 1805 roku, próbując wyznaczyć orbity asteroid Pallas i Juno.

• Co zrobił: Opracował matematyczny trik pozwalający drastycznie skrócić liczbę mnożeń potrzebnych do analizy okresowych orbit.
• Dlaczego świat o tym zapomniał: Gauss nigdy nie opublikował tej pracy za życia. Została ona wydana dopiero pośmiertnie w 1866 roku w jego Dziełach zebranych, napisana po łacinie, przez co świat nauki kompletnie ją przeoczył na kolejne sto lat.

1807 rok – Jean Baptiste Joseph Fourier (Narodziny idei)

Francuski matematyk Joseph Fourier przedstawia przed Instytutem Francuskim swoją przełomową pracę dotyczącą rozchodzenia się ciepła w ciałach stałych.

• Co zrobił: Sformułował twierdzenie, że każdą funkcję okresową można przedstawić jako sumę nieskończonego szeregu sinusów i cosinusów.
• Odrzucenie przez elity: Najwięksi matematycy epoki (Lagrange, Laplace, Legendre) odrzucają jego pracę. Lagrange kategorycznie twierdził, że zaokrąglonymi sinusami nie da się opisać funkcji o ostrych, "kwadratowych" załamaniach (tzw. efekt Gibbsa).

1822 rok – Publikacja „Analitycznej teorii ciepła”

Po 15 latach walki z recenzentami i dopracowywania teorii, Fourier publikuje swoje opus magnum: Théorie analytique de la chaleur.

• Co to zmienia: W tej książce formalnie definiuje to, co dziś nazywamy Szeregiem Fouriera (dla sygnałów okresowych) oraz Transformacją Fouriera (dla sygnałów nieokresowych, rozciągniętych od \(-\infty \) do \(\infty \)). Świat matematyki zaczyna rozumieć potęgę tego narzędzia.

1829 rok – Peter Gustav Lejeune Dirichlet (Rygor matematyczny)

Fourier działał intuicyjnie, co dawało pole do krytyki. Niemiecki matematyk Dirichlet jako pierwszy podaje ścisły dowód matematyczny zbieżności szeregów Fouriera.

• Co zrobił: Sformułował tzw. warunki Dirichleta. Określił dokładnie, jakie kryteria musi spełniać funkcja (np. skończona liczba minimów, maksimów i punktów nieciągłości), aby transformacja Fouriera była dla niej w 100% poprawna matematycznie.

1965 rok – James Cooley i John Tukey (Ponowne odkrycie FFT i cyfrowa rewolucja)

Wraz z pojawieniem się pierwszych komputerów, klasyczna Dyskretna Transformacja Fouriera (DFT) okazała się zbyt ciężka obliczeniowo – dla \(N\) próbek wymagała \(N^{2}\) operacji.

• Co zrobili: Cooley i Tukey publikują niezależnie algorytm FFT (ang. Fast Fourier Transformation – Szybka Transformacja Fouriera), redukujący liczbę obliczeń do \(N \log_2 N\).
• Efekt: To najważniejszy moment dla współczesnej technologii. Bez tego algorytmu komputery nie byłyby w stanie przetwarzać dźwięku i obrazu w czasie rzeczywistym.

Dopiero po publikacji odkryto, że powtórzyli oni schemat Gaussa z 1805 roku.

Przełom XX/XXI wieku – Powstanie wersji okienkowej (STFT) i falkowej

Klasyczna FT ma jedną wadę: pokazuje jakie częstotliwości są w sygnale, ale nie pokazuje, kiedy wystąpiły. W odpowiedzi na to ograniczenie matematycy rozwinęli transformację:

• STFT (ang. Short-Time Fourier Transform – Krótkookresowa Transformacja Fouriera): Dzieli sygnał na małe "okienka" czasowe i analizuje je po kolei.
• Ciągła transformacja falkowa (ang. Continuous wavelet transform): Zastępuje sinusy tzw. falkami, co pozwala na jednoczesną, precyzyjną analizę czasu i częstotliwości. Stała się podstawą formatów takich jak JPEG 2000.

Intuicja geometryczna: Szukanie przepisu

Wyobraźmy sobie, że mamy gotowe ciasto. Transformacja Fouriera to algorytm, który działa jak wykrywacz składników – mówi nam dokładnie, ile gramów mąki, cukru i jajek użyto do wypieku tego ciasta.

W świecie inżynierii i kodu dostajemy poszarpany, trudny do zanalizowania wykres funkcji w dziedzinie czasu \(f(t)\). Transformacja Fouriera przekształca go na wykres częstotliwości \(F(\omega)\), dając jasną informację: „Ten sygnał składa się w 70% z fali o częstotliwości 50 Hz i w 30% z szumu o częstotliwości 120 Hz”.

Jak zmierzyć ukrytą częstotliwość? (Mnożenie i sumowanie)

Aby sprawdzić, czy w nieznanym sygnale ukrywa się fala o konkretnej częstotliwości (np. 100 Hz), matematyka wykonuje prosty test dopasowania:

1. Bierzemy nasz sygnał i mnożymy go punkt po punkcie przez idealną, czystą sinusoidę o częstotliwości 100 Hz.
2. Sumujemy (całkujemy) wynik tego mnożenia z całego przedziału czasu.

Co z tego wynika?

• Brak dopasowania: Jeśli w sygnale nie ma fali 100 Hz, wykresy będą się mijać. Wartości dodatnie pomnożone przez ujemne dadzą minus, ujemne przez dodatnie dadzą minus. Po zsumowaniu wszystkiego wartości dodatnie i ujemne wzajemnie się zniosą. Wynik końcowy wyniesie zero.
• Idealne dopasowanie: Jeśli w sygnale jest ukryta fala 100 Hz, jej grzbiety spotkają się z grzbietami naszej fali testowej. Plus razy plus da plus, a minus razy minus da również plus. Wszystkie punkty zaczną współgrać w tym samym kierunku. Po zsumowaniu otrzymamy bardzo dużą wartość – to geometryczny dowód, że ta częstotliwość tam istnieje.

Skrót inżynierski: Wzór Eulera

Zamiast mnożyć sygnał osobno przez sinusy (odpowiedzialne za przesunięcie w fazie) i kosinusy (odpowiedzialne za amplitudę), matematycy używają liczb urojonych. Słynny wzór Eulera wiąże funkcje trygonometryczne z funkcją wykładniczą:

Zapis \(e^{-i\omega t}\) w transformacji Fouriera nie powinien przerażać. W interpretacji geometrycznej to po prostu "wskazówka zegara", która kręci się po okręgu płaszczyzny zespolonej w lewą stronę z określoną prędkością (częstotliwością kołową) \(\omega \).

Klasyczny wzór krok po kroku

Oto pełne, ciągłe równanie transformacji Fouriera:

Rozbijmy je na czynniki pierwsze:

• \(\omega \) – Częstotliwość kołowa (prędkość kątowa). Dla uproszczenia myślowego możemy traktować ją jak badaną częstotliwość, pamiętając jednak o inżynierskiej zależności \(\omega = 2\pi f\), gdzie \(f\) to klasyczna częstotliwość wyrażona w hercach.
• \(F(\omega)\) – Wynik końcowy dla konkretnej badanej częstotliwości \(\omega \). Jest to liczba zespolona. Jej długość (moduł) mówi nam, jak głośna jest ta częstotliwość (amplituda), a kąt mówi o jej przesunięciu (fazie).
• \(\int\limits_{-\infty}^{\infty} \dots dt\) – Całka, czyli nasze sumowanie. Zbieramy sygnał z całego dostępnego czasu, aby wyciągnąć średnią.
• \(f(t)\) – Wejściowy, skomplikowany, surowy sygnał (funkcja czasu).
• \(e^{-i\omega t}\) – Nasza kręcąca się „probówka” (wskazówka zegara) o częstotliwości \(\omega \), na którą nawijamy sygnał, żeby sprawdzić stopień ich dopasowania.

W świecie cyfrowym (np. przy obróbce dźwięku w formacie MP3 lub filtrowaniu obrazu) całkę zastępuje się zwykłą pętlą for, która sumuje dyskretne próbki – jest to tzw. Dyskretna Transformacja Fouriera (DFT), najczęściej optymalizowana algorytmem FFT.

Wartość \(\omega \) (czyli badaną częstotliwość) bierzemy bezpośrednio z sufitu – to my, jako programiści lub badacze, decydujemy, jaką częstotliwość chcemy w danym momencie sprawdzić.

Transformacja Fouriera działa jak "radiowy odbiornik analogowy" z pokrętłem strojenia. Sygnał \(f(t)\) to całe pasmo elektromagnetyczne docierające do anteny. My kręcimy pokrętłem i ustawiamy je na konkretną wartość \(\omega \) (np. 100 MHz), aby sprawdzić, czy na tej fali nadaje jakaś stacja.

Z punktu widzenia kodu i algorytmu sprawa wygląda następująco:

W wersji ciągłej (Teoria matematyczna)

Matematyczna całka sprawdza po kolei wszystkie możliwe częstotliwości od minus do plus nieskończoności. Tworzy z tego nową, ciągłą funkcję \(F(\omega)\), gdzie osią poziomą wykresu staje się częstotliwość \(\omega\).

W wersji cyfrowej / kodzie (Praktyka programistyczna)

Gdy mamy ponumerowane próbki sygnału (np. tablicę wartości kolejnych próbek), nie możemy sprawdzić nieskończonej liczby częstotliwości. Wtedy \(\omega \) generuje się automatycznie w pętli for.

Jeśli mamy \(N\) próbek sygnału, to sprawdzamy dokładnie \(N\) różnych częstotliwości, krok po kroku:

• \(\omega _0\) = 0 Hz (składowa stała, czyli średni poziom sygnału)
• \(\omega _1\) = podstawowa częstotliwość (taka, której jeden pełny okres fali mieści się w całym okienku pomiarowym)
• \(\omega _2\) = dwukrotność podstawowej częstotliwości (dwa okresy fali w okienku)
• \(\omega _3\) = trzykrotność podstawowej częstotliwości (trzy okresy fali w okienku)
• ... i tak dalej, aż do tzw. częstotliwości Nyquista (połowy częstotliwości próbkowania).

Mówiąc krótko: \(\omega\) to parametr wejściowy pętli. Pętla zewnętrzna wybiera kolejną częstotliwość testową \(\omega \), a pętla wewnętrzna „nawija” na nią cały sygnał czasu, żeby sprawdzić, co z tego wyjdzie.

Inżynierski problem z wydajnością: Koszt \(O(N^2)\)

Ponieważ dla każdej z \(N\) częstotliwości badamy \(N\) próbek, procesor komputera musi wykonać \({N \times N = N^2}\) operacji mnożenia i dodawania.

• Jeśli sygnał ma \(1000\) próbek, komputer wykona \({1{'}000{'}000}\) operacji (milion – dla współczesnych procesorów to ułamek sekundy).
• Jeśli jednak nagramy tylko 25 sekund dźwięku w jakości CD (gdzie mamy 44 100 próbek na sekundę, czyli łącznie \({1{'}102{'}500}\) próbek), to \(N^{2}\) da kosmiczną liczbę \({1{'}215{'}506{'}250{'}000}\) operacji. Tego żaden procesor nie przeliczy w locie.

Jak to oszukano? Szybka Transformacja Fouriera (FFT)

Aby uniknąć tego potężnego \(N^{2}\) sprawdzania, w latach 60. XX wieku matematycy Cooley i Tukey wymyślili genialny algorytm FFT (ang. Fast Fourier Transform).
Wykorzystuje on fakt, że "wskazówka zegara" w ułamkach zespolonych kręci się w kółko i jej wartości cyklicznie się powtarzają. Algorytm FFT dzieli całą tablicę \(N\) próbek na pół (na próbki parzyste i nieparzyste), przelicza je osobno, a potem sprytnie składa wyniki.

Dzięki temu zamiast sprawdzać wszystko \({N\times N}\) razy, liczba operacji spada do \({N \log_2 N}\).

Dla wspomnianego miliona próbek zamiast biliona operacji, FFT wykonuje ich tylko około 20 milionów – to dlatego np. telefon potrafi w ułamku sekundy narysować skaczące słupki korektora graficznego podczas odtwarzania muzyki.

xxx

xxx

xxx

xxx

xxx