FFT - szybkie Przekształcenie Fouriera

Co tak naprawdę dokonuje się w ich wnętrzach i - najważniejsze - czy wykorzystujemy ich możliwości w 100-procentach? Wydaliśmy kupę pieniędzy na zakup niezbędnego dla nas urządzenia lub oprogramowania, ale często nie potrafimy skorzystać z wielu jego opcji czy funkcji. A bierze się to z braku wiedzy na jego temat lub na temat realizowanych przez niego "zadań".

Temu pierwszemu zaradzić nie możemy - każdy potencjalny użytkownik nowozakupionego RTA czy programu powinien przede wszystkim zapoznać się z "Manualem", czyli Instrukcją Użytkowania. Jednak nawet znając od podszewki naszą nową "zabawkę", jeśli nie rozumiemy, o co w niej tak naprawdę chodzi, nie będziemy w stanie efektywnie z niej korzystać.

W ostatnim numerze zagłębialiśmy się w tajniki procedury zwanej Dyskretnym Przekształceniem Fouriera, znanym również jako DFT. Wprawdzie DFT jest najbardziej bezpośrednią procedurą matematyczną do określania częstotliwościowej zawartości sygnału, ale jest ona bardzo nieefektywna. Ponieważ liczba punktów DFT obecnie przekracza często dziesiątki tysięcy, liczba wymaganych do przetworzenia danych bardzo szybko staje się gigantyczna. Dlatego DFT w "czystej" formie nie jest używane w praktycznych zastosowaniach - służy tylko do wyjaśnienia zasady przetwarzania za pomocą tego rodzaju procedury.

W roku 1965, a więc w czasach gdy wielu z nas (z piszącym te słowa włącznie) nie było jeszcze na świecie, został opublikowany przez dwóch panów o nazwiskach rodem z niemych filmów komediowych z lat 30-tych - Cooleya i Tuckeya - artykuł przedstawiający bardzo wydajny algorytm implementujący DFT. Algorytm ten jest obecnie znany jako Szybkie Przekształcenie Fouriera (ang. Fast Fourier Transform).

Przed wprowadzeniem w życie FFT obliczanie tysiącpunktowych DFT wymagało wielu godzin pracy specjalistycznych komputerów, w wielkich centrach komputerowych lub akademickich. Teraz obliczenie 1.024-punktowej FFT zajmuje naszym domowym komputerom ułamki sekund, tak że można to robić "w locie", czyli w czasie rzeczywistym. Fakt, że technika komputerowa od tamtych czasów posunęła się drastycznie naprzód, ale faktem jest również to, że algorytm FFT jest znacznie, znacznie bardziej "szybki" w porównaniu do DFT.

ZWIĄZEK POMIĘDZY FFT I DFT

Wprawdzie zaproponowano wiele różnych algorytmów FFT, jednak najpopularniejszym i powszechnie stosowanym jest algorytm FFT o podstawie 2. Algorytm ten jest bardzo efektywną procedurą wyznaczania DFT pod warunkiem, że rozmiar DFT będzie całkowitą potęgą liczby dwa - czyli mówiąc matematycznie, liczba punktów transformaty wynosi N = 2k, gdzie k jest pewną liczbą naturalną.

Spróbuję teraz pokazać tę efektywność FFT w stosunku do tradycyjnej DFT. Jeśli przypomnimy sobie wzór na obliczenia DFT, to zauważymy, że aby obliczyć N-punktową transformatę musimy przeprowadzić N2 mnożeń (i to, żeby jeszcze było bardziej skomplikowane, mnożeń zespolonych, gdyż to one głównie pochłaniają moc obliczeniowa naszego procesora - dodawanie jest o wiele szybszym i mniej "pracochłonnym" procesem, dlatego jest tutaj pomijane).

Na przykład dla 8-punktowej DFT musieliśmy dokonać 82 czyli 64 mnożenia zespolone. W przypadku FFT okazuje się, że liczba mnożeń zespolonych dla N-punktowej FFT wynosi około N/2 * log2N. Napisałem w przybliżeniu, ponieważ niektóre mnożenia okazują się mnożeniami przez +1 lub -1, co sprowadza się jedynie do zmiany znaku.

Oznacza to znaczącą poprawę wydajności poprzez zmniejszenie liczby mnożeń w stosunku do ilości mnożeń dla klasycznej DFT. Skąd bierze się takie zmniejszenie mnożeń? Ano np. z tego, że wiele operacji, które wykonujemy przy obliczaniu transformaty Fouriera "na piechotę", ze wzoru DFT jest nadmiarowych, często np. wykonujemy wielokrotnie mnożenie tych samych czynników, podczas gdy wystarczyłoby zrobić to tylko raz, a potem przy ponownym mnożeniu tych samych liczb odczytywać już gotowy, wcześniej obliczony wynik.

Algorytm FFT o podstawie 2 eliminuje tę nadmiarowość i przez to tak znacząco zmniejsza liczbę operacji arytmetycznych. Aby lepiej to sobie wyobrazić, niech przemówią liczby, gdyż takie suche wzory przeważnie nie są zbyt przejrzyste.

MAŁY PRZYKŁADZIK

Jeśli, na przykład, chcielibyśmy obliczyć 512 punktową (pamiętamy, że przy FFT liczba punktów transformaty musi być potęgą liczby 2) transformatę Fouriera, musimy wykonać 200 razy więcej mnożeń zespolonych przy obliczaniu DFT w stosunku do wymaganych w algorytmie FFT.

Nie brzmi powalająco? Ano nie, ale jeśli będziemy mieli już nie 512, a 8.192 punkty, to w DFT musimy obliczyć 1.000 mnożeń zespolonych na każde(!) mnożenie zespolone w FFT! Teraz już lepiej, prawda? Żeby sobie to jeszcze lepiej uzmysłowić wystarczy spojrzeć na rysunek na poprzedniej stronie. Widać na nim jak na dłoni zysk, jaki daje nam FFT w stosunku do DFT.

W tym miejscu należy jasno stwierdzić, że FFT nie jest przybliżeniem DFT, tak jak to ma miejsce np. przy przetwarzaniu sygnałów analogowych na cyfrowe. Jest ona równoważna DFT i to JEST DFT w stu procentach. Ponadto wszystkie właściwości, o których pisałem w poprzednim odcinku, zarówno te pozytywne (symetria, liniowość), jak i negatywne (przeciek), dotyczą również zachowania się FFT.

Nie ma sensu zagłębiać się w strukturę algorytmu FFT. Do niczego nam to potrzebne nie będzie, a zainteresowani zawsze mogą poczytać o tym w fachowej literaturze, której jest na kilogramy, jeśli nie na tony wręcz. Zajmijmy się może bardziej praktycznymi wskazówkami dotyczącymi FFT.

PRÓBKOWANIE PRZY FFT

Jak wiadomo, przy przetwarzaniu A/C, aby ustrzec się przed aliasingiem w dziedzinie częstotliwości, stosuje się kryterium Nyquista, według którego powinniśmy próbkować sygnał analogowy z częstotliwością co najmniej dwa razy większą od szerokości pasma tego sygnału.

W zależności od zastosowania próbkuje się z szybkością 2,5 do 4 razy większą, niż szerokość pasma (pomijając oversampling). Jeśli znamy szerokość pasma, sprawa jest prosta. Gorzej, jeśli szerokość pasma wejściowego sygnału analogowego jest nam nieznana. Jak wtedy możemy stwierdzić czy występują problemy z aliasingiem?

Przede wszystkim należy podchodzić z dużą rezerwą do wyników FFT, które mają znaczące składowe widmowe przy częstotliwościach bliskich połowie częstotliwości próbkowania. W przypadku idealnym chcielibyśmy mieć do czynienia z sygnałami, których wartości widmowe maleją ze wzrostem częstotliwości. Takie w zasadzie są sygnały audio, których wartości widmowe maleją niemalże do zera przy częstotliwościach wyższych od 16 kHz.

Jeśli więc będziemy stosować częstotliwość próbkowania przynajmniej równą 44,1 kHz, możemy być stosunkowo spokojni. Jeśli jednak istnieją jakiekolwiek składowe widmowe, których częstotliwości okazują się być zależnymi od częstotliwości próbkowania, wtedy powinniśmy podejrzewać powstanie aliasingu. Aby temu zapobiec musimy przed próbkowaniem stosować analogowy, dolnoprzepustowy filtr antyaliasingowy, o czym już kiedyś była mowa.

ROZDZIELCZOŚĆ

Wprawdzie wiemy, że N-punktowa FFT o podstawie 2 wymaga N = 2k próbek wejściowych, ale należy zadać sobie pytanie, ile próbek musimy zebrać zanim wyznaczymy FFT? Przedział czasu, w którym będziemy pobierali dane musi być wystarczająco długi, aby spełnione były nasze wymagania odnośnie rozdzielczości częstotliwościowej przy zadanej szybkości próbkowania. Przedział czasu, w którym pobieramy dane jest odwrotnością wymaganej rozdzielczości FFT i im dłużej próbkujemy z ustaloną szybkością próbkowania, tym lepsza będzie rozdzielczość częstotliwościowa.

Jeśli przy ustalonej częstotliwości próbkowania fs potrzebujemy konkretnej rozdzielczości r, to ilość potrzebnych próbek obliczymy ze wzoru:

N = fs/r

W przypadku gdy częstotliwość próbkowania wynosi 44,1 kHz, a chcemy uzyskać rozdzielczość 5 Hz, wtedy

N = fs/5 = 0,2 fs = 0,2 × 44.100 = 8.820

Skoro N musi być potęgą 2, to najbliższa wartość wynosi 8.192. I tutaj mamy następny problem - czy wybrać tę wartość, która jest najbliższa, ale będzie wymagała obcięcia kilkuset próbek, czy "podciągnąć" w górę? Zasada jest taka, że zawsze lepiej "w górę". Dzięki temu unikniemy zdegradowania rozdzielczości. No więc dobrze, w górę. Ale skąd wziąć te brakujące? Zwłaszcza że następna dostępna liczba N to 16.384? A więc brakuje nam, bagatela, 7.564 próbek!

Fakt, że przykład niezbyt udany, bo akurat w tym przypadku może i lepiej byłoby odciąć te 628 próbek, niż dodawać zerowe, bo tak właśnie możemy postąpić. Ale jeśli np. dysponujemy 12.000 próbek, to w takim przypadku stanowczo lepiej dodać do wymaganych w następnym progu 16.384 próbek (doklejając na końcu próbki zerowe), niż ujmować.

OKIENKOWANIE

FFT cierpi na takie same negatywne skutki przecieku widmowego, jakie opisywałem w ubiegłym miesiącu przy DFT. Można więc, a nawet trzeba, przemnażać dane czasowe przez funkcje okien, aby złagodzić ten problem, choć powinniśmy się liczyć i zawsze pamiętać o nieuchronnej degradacji rozdzielczości częstotliwościowej, zawsze występującej, jeśli używane są okna.

Przy okazji dołączania próbek zerowych, jeśli ma to miejsce, musimy być pewni, że dołączamy zera po przemnożeniu oryginalnego ciągu danych czasowych przez funkcję okna. Zastosowanie funkcji okna do ciągu próbek uzupełnionych próbkami o zerowych wartościach zniekształci wynikowe okno i powiększy przeciek FFT.

Wprawdzie okienkowanie zmniejsza przeciek, jednak nie eliminuje go całkowicie. Nawet gdy wykorzystuje się okienkowanie, składowe widmowe o wysokim poziomie mogą maskować pobliskie składowe widmowe o poziomie niskim. Jest to szczególnie wyraźnie widoczne, kiedy oryginalne dane czasowe mają niezerową wartość średnią, tj. są przesunięte o składową stałą. Jeśli w takim przypadku wyznaczamy FFT, to składowa stała przy częstotliwości 0 Hz o dużej amplitudzie przesłoni jej widmowych sąsiadów.

Możemy wyeliminować ten efekt przez obliczenie wartości średniej okienkowanego sygnału i odjęcie tej wartości średniej od każdej próbki oryginalnego sygnału okienkowanego. Operację te nazywa się często centrowaniem sygnału. Sprawia ona, że wartość średnia nowego ciągu czasowego staje się równa zeru i eliminuje w wynikach FFT każdą składową wysokiego poziomu przy częstotliwości 0 Hz.

UŚREDNIANIE

Jeśli używamy FFT, aby wykryć sygnał w obecności szumu, i dostępna jest wystarczająca liczba danych w dziedzinie czasu, możemy poprawić czułość przetwarzania przez uśrednianie wielokrotnych procedur FFT. Ta technika może być zaimplementowana bardzo efektywnie, aby wykryć sygnał, który w rzeczywistości znajduje się poniżej uśrednionego poziomu szumu. W ten sposób, mając wystarczająco dużo danych w dziedzinie czasu, możemy wykryć składowe sygnału, które charakteryzują się ujemną wartością stosunku sygnał/szum.

Na koniec małe

PODSUMOWANIE

Transformacja Fouriera w wydaniu cyfrowym, czy FFT (bo tego algorytmu używamy na co dzień w różnego rodzaju urządzeniach czy programach), jest doskonałym narzędziem, służącym do analizy częstotliwościowej badanego sygnału. W zależności od tego, ile punktów (prążków) wyjściowych chcemy uzyskać poprzez analizę, tyle też próbek wejściowych musimy posiadać. A więc aby mieć N-punktową FFT, musimy dysponować N próbkami wejściowymi. Ilość punktów transformaty oczywiście wpływa na wynikową rozdzielczość analizy, a więc im ich więcej, tym lepiej.

Ale oczywiście nie możemy zwiększać ich w nieskończoność, gdyż wraz ze wzrostem liczby punktów transformaty rośnie nam ilość wykonywanych operacji, a przez to wydłuża się czas analizy (dlatego przy zbyt "gęstej" FFT nie zawsze będziemy mogli robić ją w czasie rzeczywistym - to oczywiście zależy w dużej mierze od możliwości sprzętowych, głównie chodzi o możliwości procesora).

W przypadku "klasycznej" DFT ilość mnożeń potrzebnych do obliczenia N-punktowej transformaty wynosi N2, a w przypadku algorytmu FFT o podstawie 2 ilość tych operacji określona jest wzorem N/2 * log2N. Zarówno w klasycznej DFT, jak i FFT jesteśmy narażeni na zjawisko przecieku, które to zjawisko fałszuje nam wynik końcowy analizy. Zamiast więc otrzymywać cienkie prążki odpowiadające konkretnym częstotliwościom, o wysokościach zależnych od wielkości składowej tej konkretnej częstotliwości, otrzymujemy mniej lub bardziej "rozlane" słupki.

Tak więc końcowy wykres wygląda często jak rozlany budyń. Możemy to zjawisko minimalizować, stosując okienkowanie, jednak nigdy nie uda nam się go pozbyć całkowicie. Okna zmniejszają co prawda rozdzielczość, ale poprzez zmniejszanie skutków przecieku często jesteśmy w stanie wykryć małe składowe (zwłaszcza te, które znajdują się w bliskim sąsiedztwie dużych składowych), które normalnie "giną' nam w przecieku FFT czy DFT.

Dobór rodzaju okna oraz jego długości jest więc kompromisem, który każdy sobie musi wypracować w zależności od potrzeby analizy. A żeby to wypracować, trzeba próbować, próbować i jeszcze raz próbować. A więc do dzieła!

Jan Erhard

---

Jan Erhard z wykształcenia jest informatykiem i specjalistą od sieci cyfrowych, zaś z zamiłowania muzykiem. Zajmuje się też realizacją dźwięku, stąd jego zainteresowania i duża wiedza na temat cyfrowego przetwarzania sygnałów.