DFT - przekształcenie Fouriera w wersji cyfrowej

Ta druga to filtracja cyfrowa - o której też sobie opowiemy w jednym z kolejnych artykułów serii - zaś ta pierwsza, która będzie nas interesować w tym numerze, to

DYSKRETNE PRZEKSZTAŁCENIE FOURIERA

czyli DFT. Oczywiście, skrót DFT pochodzi z angielskiej "wersji językowej", czyli od słów Discrete Fourier Transform. Ktoś powie: owszem, transformacja Fouriera - to obiło mi się o uszy, ale pierwszy raz słyszę o jakimś DFT! Co to takiego jest? No tak, może bardziej "przemawiający" będzie więc skrót FFT? Czy teraz już lepiej?

No więc FFT i DFT to prawie to samo. Mówiąc najoględniej, FFT to szczególny przypadek algorytmu DFT, który dzięki zastosowaniu pewnego "tricku" jest znacznie szybszy od standardowego algorytmu DFT, a więc efektywniejszy. Może jednak czas zacząć operować konkretami.

Niestety, nie da się ukryć, że DFT (podobnie jak i przekształcenie Fouriera w ogólnym rozumieniu) jest procedurą matematyczną - nie uda nam się więc wymigać się od "odrobiny" matematyki. Postaram się jednak, aby była ona przystępna nawet dla tych, którzy na widok całek i innych "ślimaczków" dostają drgawek.

Przede wszystkim należy powiedzieć, że DFT służy do wyznaczania zawartości harmonicznej lub częstotliwościowej sygnału dyskretnego (te dwa ostatnie wyrazy są bardzo istotne w tym stwierdzeniu). Pochodzenie DFT bierze oczywiście swój początek od ciągłego przekształcenia Fouriera, zdefiniowanego jako:

gdzie x(t) jest pewnym sygnałem ciągłym w dziedzinie czasu (czyli np. napięciem mierzonym na zaciskach wyjściowych jakiegoś układu analogowego).

Z innych występujących we wzorze "ślimaczków" mogą nas zaniepokoić literki:j =  √-1, czyli coś, co tak naprawdę nie istnieje (o ile pamiętamy ze szkoły, pierwiastek kwadratowy z liczby ujemnej nie istnieje), czyli jest to kolejny wymysł chorych umysłów matematyków. A tak na poważnie, jest to operator używany w rachunku liczb zespolonych - nie musimy bynajmniej zawracać sobie nim zbytnio głowy.

e – czyli podstawa logarytmu naturalnego.

Dzięki temu przekształceniu możemy przetransformować sygnał x(t) z dziedziny czasu w dziedzinę częstotliwości, czyli, mówiąc "po ludzku", wyznaczyć jego widmo w postaci ciągłej. To z kolei pozwala nam określić zawartość częstotliwościową dowolnego sygnału nas interesującego. Dowolnego?

Teoretycznie tak, ale spróbujcie obliczyć sobie takie wyrażenie:

no, i jak idzie? Prawda, że to nie takie proste?

A przecież w tym przypadku nasz sygnał x(t) to zwykły sygnał sinusoidalny o częstotliwości 1 kHz, czyli ton prosty. Co dopiero, jeśli nasz x(t) jest złożeniem kilku, kilkudziesięciu, a nawet kilkuset (nie mówiąc już o kilku tysiącach) różnych sygnałów  sinusoidalnych o różnych częstotliwościach i różnych fazach początkowych? No, bo tak po prawdzie rzeczywiste sygnały audio, nawet te wydawałoby się proste, jak głos ludzki czy pojedyncze instrumenty, to właśnie takie sygnały, będące sumą bardzo wielu składowych sinusoid.

Dzięki komputeryzacji i cyfryzacji zyskaliśmy potężnego sojusznika, który nie tylko pozwala uzyskać dźwięk czysty, klarowny, pozbawiony szumów, brumów i zakłóceń, ale i też daje nam możliwość w o wiele łatwiejszy sposób dokonywać różnych operacji na sygnale spróbkowanym, między innymi transformacji Fouriera - szybciej, łatwiej i przyjemniej.

Dlaczego? Ano dlatego, że w sygnale cyfrowym nie musimy już dokonywać całkowania, lecz zwykłego sumowania, i to też nie od -∞ do +∞, ale w ograniczonym jego zakresie. Czas więc ujawnić

RÓWNANIE OKREŚLAJĄCE DFT

a wygląda ono następująco:

... oznacza sumowanie (czyli dodawanie) kolejnych wyrażeń pod znakiem sumy, zmieniając parametr n (i odpowiadające mu wartości próbek x(n)) od wartości 0 do N-1, wartość N zaś zakładamy z góry.

Owszem, i w tym równaniu mamy zarówno "magiczne j", jak i "bolesne e". Jednak, jak już wspomniałem, nie musimy sobie zawracać głowy owym "j" - jest to tylko wygodna (jak dla kogo!!) abstrakcja, która pomaga porównywać zależności fazowe pomiędzy różnymi składowymi sinusoidalnymi sygnału. A z "bolesnym e" można sobie poradzić (czyt. wyeliminować je) w pewien sprytny sposób.

Z tej "wzorologi" wynikają dwa, bardzo ważne dla nas wnioski. Po pierwsze: jeśli potraktujemy DFT jako pewną "czarną skrzynkę", w której dokonywane są pewne operacje, dzięki którym sygnał w dziedzinie czasu na wejściu jest zamieniany na sygnał w dziedzinie częstotliwości na wyjściu, to zarówno sygnał, który wchodzi do tej "skrzynki", jak i wychodzący z niej są sygnałami dyskretnymi, czyli spróbkowanymi (mówiąc inaczej, są to sygnały cyfrowe) - Rysunek 1.

Po drugie: liczba próbek wychodzących z naszej "czarnej skrzynki" jest zawsze taka sama, jak liczba wchodzących. Wydawałoby się na pierwszy rzut oka, że to oczywiste. Ale weźmy pod uwagę, że próbki wejściowe to spróbkowane wartości pewnego ciągłego sygnału rzeczywistego (jakim może być np. napięcie), natomiast na wyjściu otrzymujemy "prążki" oznaczające pewne konkretne wartości na skali częstotliwościowej (a jak się okaże, odpowiadające wartościom spróbkowanego widma, które otrzymalibyśmy, obliczając ciągłe przekształcenie Fouriera z naszego nie spróbkowanego jeszcze, ciągłego sygnału wejściowego).

Może trochę zamieszałem, ale najważniejsze jest to, że mając N próbek wejściowych w dziedzinie czasu, DFT wyznacza zawartość widmową sygnału wejściowego w N równomiernie rozłożonych punktach osi częstotliwości. Wartość N określa nam więc ile wymaganych jest próbek wejściowych, aby osiągnąć odpowiednią rozdzielczość wyników w dziedzinie częstotliwości. Dzięki niej będziemy też mogli określić czas potrzebny do obliczenia N-punktowej DFT.

Jak wszystko na tym świecie - wszak nie ma ideałów - tak i DFT nie jest wolne od pewnych problemów. Jednym z nich, i bodajże najważniejszym, jest tzw.

PRZECIEK DFT

Choć brzmi to bardzo "hydraulicznie", nie ma to nic wspólnego z wodą ani inną cieczą. DFT daje prawidłowe wyniki jedynie wtedy, kiedy ciąg danych wejściowych zawiera energię rozłożoną dokładnie przy częstotliwościach będących całkowitymi wielokrotnościami częstotliwości podstawowej, określonej wzorem fS/N. Jeśli zaś sygnał wejściowy zawiera składową o pewnej częstotliwości pośredniej, to ujawni się ona w pewnym stopniu przy wszystkich N prążkach (czyli wszystkich częstotliwościach wyjściowych), dla których przeprowadzamy analizę tego sygnału.

Oczywiście, najbardziej będzie on wpływał na prążki będące w pobliżu "jego" częstotliwości (np. jeśli mamy prążki co 1 kHz, a składowa ma wartość 3,6 kHz, to wpłynie ona najbardziej na prążki 2, 3, 4 i 5 kHz), jednakże wpływa to też na wartość nawet bardzo oddalonych prążków (proporcjonalnie do ich "oddalenia").

MAŁY PRZYKŁAD

Załóżmy, że mamy sygnał wejściowy, będący spróbkowanym sygnałem sinusoidalnym o 3 okresach zawartych w N = 64 próbkach. Tyle próbek będzie nam potrzebnych, gdyż chcemy wyznaczyć 64 punktową DFT tego przebiegu. Na rysunku 2a pokazano wygląd takiego spróbkowanego sygnału sinusoidalnego, a na rysunku 2b wynik przeprowadzonej analizy DFT. Jak widać, sygnał nasz nie zawiera składowej stałej (a więc i jego wartość średnia jest również zerowa - nie ma prążka w punkcie 0) ani składowej częstotliwościowej innej niż częstotliwość odpowiadająca m = 3. A więc wszystko się zgadza.

Załóżmy teraz, że chcemy dokonać DFT ciągu, który w zakresie N = 64 próbek zawiera 3,4 okresów (rysunek 3a). Ponieważ sygnał wejściowy nie ma całkowitej liczby okresów w przedziale 64 próbek, energia wejściowa "przecieka" do wszystkich innych prążków, jak to pokazuje Rysunek 3b. Niezerowy jest nie tylko prążek odpowiadający m = 3, ale też i m = 1, 2, 4, 5, 6, a nawet prążek 0. Również do prążków o wartości m > 6 "przecieka" energia, choć oczywiście w mniejszym już stopniu.

To jest właśnie przeciek - powoduje on, że dowolny sygnał wejściowy, którego częstotliwości nie jest dokładnie równa częstotliwość, dla której jest wyznaczany dany prążek DFT, przecieka do wszystkich innych wyznaczanych prążków DFT. Co więcej, przeciek jest nie do uniknięcia, kiedy wyznaczamy DFT rzeczywistego ciągu czasowego o skończonej długości.

Nie wnikajmy szczegółowo, dlaczego akurat tak, a nie inaczej wyglądają prążki sygnału "przeciekającego", i dlaczego tak się dzieje, że sygnał, który w przedziale analizy nie ma całkowitej liczby okresów, nie będzie jednoznacznie określony, ale będzie się "rozpływał" po całym przedziale częstotliwości. Wspomnę tylko, że jest to związane z funkcją sinc = (sin x)/x, którą już znamy z poprzednich artykułów z tej serii (dla przypomnienia, wygląda ona tak, jak na rysunku 4a). Z tego powodu tworzą się tak zwane listki boczne, pod którymi "prześlizgują" się prążki dla sygnałów "przeciekających", co uwidacznia rysunek 4b.

Zobaczmy lepiej, co zrobić, aby zminimalizować (gdyż całkiem wyeliminować niestety nie możemy) przeciek DFT. Otwórzmy więc

OKNA

i wpuśćmy trochę wiedzy na ten temat.

Okienkowanie to proces mający na celu zmniejszenie przecieku poprzez zminimalizowanie wspomnianych nieco wyżej listków bocznych funkcji sinc = (sin x)/x. Dokonujemy tego przez wymuszenie na wartościach wejściowego ciągu czasowego, zarówno na początku, jak i na końcu przedziału próbkowania, aby gładko zmierzały do pojedynczej wspólnej wartości.

Po pierwsze jednak okno służy nam do "wykrojenia" z długiego ciągu próbek tej liczby N próbek, które są nam potrzebne do obliczenia N-punktowej DFT. Wystarczy proste obliczenie: jeśli mamy fragment nagrania próbkowanego z częstotliwością 44,1 kHz, to oznacza, że w 1 s tego nagrania dysponujemy 44.100 próbkami. Jeśli będziemy mieć fragment 10 sekundowy, to automatycznie tych próbek będzie 10 razy więcej.

Wyznaczając więc 1.024-punktową DFT będziemy potrzebować tylko 1.024 próbki w jednej "porcji", będzie nam więc potrzebny fragment nagrania mający w przybliżeniu długość 0,023 s. To tylko tak dla orientacji, gdyż program dokonujący analizy sam sobie "wytnie", i to nie "około 0,023 s", ale dokładnie 1.024 próbki. To wycięcie to nic innego, jak przemnożenie tego długiego ciągu kilku tysięcy próbek przez okno o długości N próbek.

Spójrzmy na rysunek 5a, który przedstawia przykładowy, teoretycznie nieskończony (a praktycznie bardzo, bardzo długi) ciąg próbek wejściowych. Poniżej, na rysunku 5b, jest funkcja okna prostokątnego (rectangle), która ma wartość 1 w przedziale próbkowania i 0 poza tym przedziałem. To oznacza, że jeśli teraz przemnożymy nasz sygnał wejściowy przez to okno, to w przedziale próbkowania (gdzie wartość okna, przypomnę, ma wartość 1) sygnał wejściowy pozostanie bez zmian, a poza tym przedziałem wyzerowany (no, bo jakiegoż innego wyniku mogliśmy się spodziewać, mnożąc cokolwiek przez 0).

Reasumując, wycięliśmy z sygnału wejściowego tyle, ile nam trzeba. Wychodzi więc to, co pokazuje rysunek 5c. Za każdym więc razem, gdy wyznaczamy DFT sygnału wejściowego o skończonym czasie trwania, to tak jakbyśmy przemnażali ten sygnał przez okno samych jedynek, a wartości wejściowe poza tym oknem mnożyli przez zero. Jak się jednak okazuje, właśnie z tego powodu robi się cały ten galimatias.

Wspominana wcześniej funkcja sinc, przez którą przeciek ma miejsce, jest właśnie transformatą Fouriera funkcji prostokątnej, czyli naszego okna prostokątnego. I właśnie z powodu nieciągłych zmian okna prostokątnego pomiędzy wartością jeden a wartością zero biorą się te nieszczęsne listki boczne funkcji sinc. Aby więc zminimalizować przeciek widma, spowodowany przez listki boczne, musimy ich wartość zmniejszyć, używając okna innego niż prostokątne.

Jeżeli np. przemnożymy nasz sygnał wejściowy przez okno trójkątne (triangular), takie jak pokazuje rysunek 5d, to otrzymamy okienkowany sygnał, który będziemy poddawać DFT w formie pokazanej na rysunku 5e. Zauważyć możemy od razu, że wartości tego wynikowego sygnału wyjściowego stają się takie same na początku i na końcu przedziału próbkowania.

Zredukowana nieciągłość zmniejsza poziomy prążków listków bocznych. To oczywiście nie jedyne okno, które możemy zastosować. Oprócz okna trójkątnego powszechnie stosuje się okno Hanninga (rysunek 5f oraz wynikowy sygnał powstały po przemnożeniu okna Hanninga z ciągiem wejściowym - rysunek 5g), a także podobne i z nazwy, i z "wyglądu" okno Hamminga (rysunek 5h), zwane też podniesionym oknem Hanninga, co widać jeśli spojrzymy na podstawę okna.

Przyjrzyjmy się nieco bliżej nie tyle samym oknom, co efektom ich działania. Pomoże nam w tym rysunek 6. Charakterystyka okna prostokątnego stanowi miarę porównawczą dla oszacowania poprawy działania danego okna w stosunku do okna prostokątnego. Na rysunku 6a charakterystyka amplitudowa okna prostokątnego narysowana jest linią kropkową. Na rysunku tym ujęto też charakterystyki okna trójkątnego, Hanninga oraz Hamminga.

Widać od razu "gołym okiem", że te trzy okna dają zmniejszony poziom listków bocznych względem okna prostokątnego. Ponieważ okna Hamminga, Hanninga i okno trójkątne zmniejszają poziomy sygnałów poddawanych DFT, wartości maksymalne ich listków głównych są zmniejszone względem okna prostokątnego. Poza tym na skali liniowej trudno nam zaobserwować szybkość zanikania listków bocznych.

Dlatego lepiej przedstawiać tego typu charakterystyki w unormowanej (względem wartości maksymalnej okna prostokątnego - wartość maksymalna każdego okna będzie wynosić 0 dB) skali logarytmicznej, co też przedstawia nam rysunek 6b. Teraz dokładnie widzimy, że listek główny okna prostokątnego jest najwęższy, a więc i rozdzielczość DFT przy okienkowaniu za pomocą okna prostokątnego będzie lepsza od pozostałych.

Ale to koniec zalet tego okna. Zauważmy bowiem, że pierwszy listek boczny leży jedynie o -13 dB poniżej szczytu listka głównego, a ponadto spadek poziomu kolejnych listków jest bardzo powolny, tak że np. listek 7 ma poziom tylko o -25 dB mniejszy od poziomu listka głównego. Co to oznacza?? Ano to, że pomimo nienajgorszej rozdzielczości pod listkami bocznymi będzie nam przeciekać tyle sygnału, że pośród nich może nam zginąć prążek sygnału, który będzie miał o wiele mniejszą amplitudę, leżący w pobliżu sygnału o większej amplitudzie, tak jak to widać na rysunku 7.

Okno trójkątne ma zmniejszone poziomy listków bocznych, ale ceną jaką za to musimy zapłacić jest to, że szerokość listka głównego okna trójkątnego jest dwa razy większa, niż szerokość listka głównego okna prostokątnego. Również szerokość listków bocznych jest dwa razy większa, niż listków pozostałych okien, a więc, jak się okazuje, okno trójkątne jest lepsze pod względem zmniejszania przecieku, niż okno prostokątne, ale są jeszcze lepsze.

Do takich należeć będzie okno Hanninga, w przypadku którego pierwszy listek boczny jest na poziomie -31 dB w stosunku do listka głównego, a kolejne listki maleją szybko, tak że wygasają już niemalże po szóstym, siódmym listku. Natomiast okno Hamminga ma jeszcze niższy poziom pierwszego listka bocznego (na poziomie -40 dB), jednak wygasanie kolejnych listków jest o wiele, wiele wolniejsze, niż w przypadku poprzednim.

Oznacza to, że przeciek w odległości trzech lub czterech prążków od prążka środkowego jest mniejszy dla okna Hamminga, niż dla okna Hanninga, przeciek zaś dla mniej więcej tuzina prążków od prążka środkowego jest mniejszy dla okna Hanninga, niż dla okna Hamminga.

Jeśli więc zastosujemy, dla przykładu, okno Hanninga dla sygnału złożonego z dwóch składowych (jednej o dużym poziomie, a drugiej o małym poziomie w pobliżu tej pierwszej) to, jak się można spodziewać, zauważymy nie tylko mniejszy poziom prążków i zawężenie przecieku w stosunku do okna prostokątnego, ale także uda nam się z tego przecieku "wyłuskać" ów sygnał o małym poziomie, który w przypadku okna prostokątnego praktycznie całkowicie ginął nam w przecieku (rysunek 7).

To nie wszystkie okna. Powstało ich już tak wiele, że zostały one nazwane praktycznie nazwiskiem prawie każdego autora działającego w obszarze cyfrowego przetwarzania sygnałów. Nie każde z nich jest jednak tak powszechnie stosowane, jak podane powyżej, oraz kilka innych (np. Blackmana, Blackmana-Harrisa czy Welscha - zwanego też oknem gaussowskim). Jeśli nieco poeksperymentujemy w wyborze okna, zauważymy, że różne okna mają swoje zalety i wady.

Choć różnice pomiędzy nimi nie są jakieś kolosalne, to jednak dzięki odpowiedniemu doborowi okna będziemy mogli uzyskać odpowiedni kompromis pomiędzy szerokością listka głównego, poziomami pierwszego listka bocznego oraz tego, jak szybko maleją listki boczne wraz ze wzrostem częstotliwości.

Praktycznie rzecz ujmując trzeba stwierdzić, iż całościowa rozdzielczość częstotliwościowa i czułość na wykrywanie małych składowych sygnału są uzależnione nie tylko jedynie od długości DFT, ale i w dużej mierze od długości i kształtu funkcji okna. Stosowanie okien przy DFT to nie jedyne pole, gdzie okienkowanie jest bardzo przydatne. W dziedzinie przetwarzania cyfrowego audio okna wykorzystuje się również przy projektowaniu filtrów cyfrowych (o czym będziemy jeszcze mówić w jednym z najbliższych odcinków cyklu).

To tyle o oknach, ale nie koniec o DFT. W kolejnym artykule, za miesiąc, wrócimy do tematu. Powiemy też sobie co to jest FFT i czym różni się (jeśli w ogóle się różni) od DFT.

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.