Eksperymenty z FPGA (12). Dyskretna transformata Fouriera

Eksperymenty z FPGA (12). Dyskretna transformata Fouriera
Pobierz PDF Download icon

Zanim połączymy bloki naszego układu FPGA w układ analizatora widma dźwięku, rozpoczniemy od teoretycznych podstaw dyskretnej transformaty Fouriera. Jest to konieczne dla dobrego zrozumienia tego trudnego zagadnienia.

Przed przystąpieniem do wykonywania eksperymentów, zachęcam do aktualizacji repozytorium z przykładami (poprzez wywołanie polecenia git pull) [1]. Rysunek 1 pokazuje przykładowy sygnał analogowy.

Rysunek 1. Ciągły sygnał wejściowy z zaznaczonymi miejscami próbkowania

Został on spróbkowany w ośmiu, równo od siebie oddalonych punktach – są one zaznaczone na wykresie kropkami. Natomiast to, co zostaje z sygnału za przetwornikiem ADC, ilustruje rysunek 2.

Rysunek 2. Spróbkowany sygnał z rysunku 1

Na pierwszy rzut oka ciężko dopatrzeć się w pojedynczych punktach pierwotnego sygnału, ale jeżeli spełnione jest twierdzenie o próbkowaniu, to zawarta jest w nich pełna informacja o przebiegu ciągłym. Twierdzenie mówi o tym, że częstotliwość próbkowania musi być dwa razy większa niż szerokość pasma badanego sygnału.

Wrócimy jednak do naszych próbek, które możemy przedstawić jako wektor o długości N=8:

W tym zapisie współczynniki odnoszą się do najprostszej bazy, składającej się z wektorów:

Krócej możemy zapisać ją za pomocą delty Kronecker:

(1)

Gdzie vk to kolejny wektor bazy, natomiast n jest indeksem elementu wewnątrz niego. Nie jest to jednak jedyna możliwa baza. Możemy na przykład zaproponować taką, która składa się z ciągów odpowiadających kolejnym (zespolonym) sygnałom sinusoidalnym:

 (2)

Kolejne wektory bazy pokazuje rysunek 3. Na wszystkich wykresach kolorem niebieskim oznaczono część rzeczywistą, a zielonym – urojoną.

Rysunek 3. Baza złożona z zespolonych funkcji sinusoidalnych

Jeżeli więc zapiszemy nasz wektor w nowej bazie, kolejne współczynniki (teraz już zespolone) będą oznaczały amplitudę i fazę kolejnych składowych, występujących w sygnale. Zwracam uwagę, że częstotliwość k jest przy próbkowaniu równoważna częstotliwości k–N:

 (3)

Ostatnia równość jest prawdziwa, ponieważ n jest zawsze liczbą całkowitą, a dodanie 2π do argumentu funkcji trygonometrycznej, nie powoduje zmiany wyniku. Dlatego częstotliwości 5, 6 i 7 możemy interpretować także jako –3, –2 i –1 – takie podejście pokazuje rysunek 4.

Rysunek 4. Baza z rysunku 3, ale przy interpretacji części częstotliwości jako ujemnych

Warto zwrócić uwagę, że częstotliwości ujemne różnią się od dodatnich jedynie znakiem części urojonej. Jest to możliwe, ponieważ tak naprawdę interesują nas jedynie wartości w punktach, a nie całe krzywe. Zostały one narysowane na wykresie, ponieważ bez oznaczeń ciężko dojść, skąd się wzięły oraz co reprezentują poszczególne punkty – o czym można przekonać się, spoglądając na rysunek 5.

Rysunek 5. Wektory bazy bez zaznaczonych pomocniczych przebiegów

Aby przejść z jednej bazy do drugiej, potrzebujemy jeszcze macierzy przejścia. Składa się ona z kolejnych wektorów nowej bazy, wyrażonych we współczynnikach starej. W naszym przypadku jest ona dana wzorem:

 (4)

Oznacza to, że jej kolejne elementy to:

 (5)

Aby uzyskać wartości współczynników naszego wektora x, w nowej bazie musimy pomnożyć go przez macierz:

 (6)

Rozpisując, otrzymujemy:

 (7)

Po podstawieniu elementów macierzy z (5) otrzymujemy dobrze znany wzór na dyskretną transformatę Fouriera:

 (8)

Gdzie zarówno n, jak i k, zmieniają się od 0 do N–1. Dla uproszczenia zapisów czasami stosuje się też oznaczenie:

 (9)

Prowadzi to do krótszego zapisu:

 (10)

Warto zauważyć, że dla sygnału rzeczywistego części urojone muszą się zredukować. Oznacza to, że współczynniki przy częstotliwościach dodatnich będą sprzężeniem odpowiadających im współczynników przy częstotliwościach ujemnych. Nie dotyczy to częstotliwości 0 oraz , które nie mają odpowiadających im składowych ujemnych. Nie jest to jednak problemem, ponieważ oba przebiegi są rzeczywiste.

Pozostaje jeszcze wykonanie obliczeń dla danych z rysunku 1. Otrzymana transformata to:

Kształt poszczególnych składowych pokazuje rysunek 6.

Rysunek 6. Składowe sygnału z rysunku 1

FFT, czyli sprytne obliczenia

Korzystając bezpośrednio z definicji, dojdziemy do tego, że aby obliczyć pojedynczy współczynnik, musimy wykonać N mnożeń zespolonych.

Ponieważ samych współczynników też jest N, łatwo obliczyć, że przy naiwnym podejściu niezbędnych jest N2 mnożeń zespolonych. Proponuję zamiast tego prosty trik, który przyspieszy obliczenia:

(11)

Stosując inne podstawienie, można uzyskać inne wersje szybkiej transformaty Fouriera [2]. My jednak prześledzimy tylko podstawowy przypadek, w którym uzyskamy obliczenie transformaty o rozmiarze N za pomocą dwóch transformat o połowę krótszych. Podstawiając do (10), otrzymamy:

 (12)

Na początku zajmijmy się samą potęgą. Po wymnożeniu nawiasów otrzymamy:

 (13)

Sumę w wykładniku możemy rozpisać jako mnożenie:

 (14)

Spróbujmy uprościć poszczególne składniki:

 (15)

Dzieje się tak, ponieważ iloczyn k1n1 może przyjąć jedynie wartość 0 albo 1. Jeszcze bardziej możemy uprościć trzeci iloczyn:

 (16)

Ponieważ zarówno k2, jak i n1 są liczbami całkowitymi, mamy zawsze do czynienia z obrotem o pełny kąt 2π. Dlatego cały czynnik upraszcza się i jest zawsze równy 1. Korzystając z tych zależności, możemy teraz rozpisać (4) jako:

(17)

Rozpiszmy jeszcze wewnętrzną sumę. Ma ona w końcu jedynie dwa składniki:

 (18)

Po usunięciu mnożeń razy 1 otrzymamy:

 (19)

Gdy zasłonimy sobie nawias kwadratowy i przyjrzymy się pozostałej części wzoru, zauważymy coś ciekawego. Otóż cała reszta przypomina także transformatę Fouriera, ale tym razem o długości o połowę mniejszej. Rozpiszmy teraz jeszcze indeks k, osobno dla parzystych i nieparzystych współczynników, a otrzymamy:

 (20)

Rozpisaliśmy transformatę Fouriera o długości N na dwie transformaty o mniejszej długości. Teraz musimy wykonać pewne działania na współczynnikach, a następnie policzyć jedną transformatę dla współczynników parzystych, a drugą dla nieparzystych. Jednak sądzę, że to dobry moment, aby zakończyć ten odcinek i dać czytelnikom czas na spokojne przeanalizowanie opisanych obliczeń.

Rafał Kozik
rafkozik@gmail.com

Przypisy:
[1] Repozytorium z przykładami: http://bit.ly/33uYPxs
[2] Lyons r.G., Wprowadzenie do cyfrowego przetwarzania sygnałów, Wydawnictwo Komunikacji i Łączności, Warszawa 2010
[3] Film prezentujący działanie FFT: https://bit.ly/35H0yTp

Artykuł ukazał się w
Listopad 2020
DO POBRANIA
Pobierz PDF Download icon

Elektronika Praktyczna Plus lipiec - grudzień 2012

Elektronika Praktyczna Plus

Monograficzne wydania specjalne

Elektronik kwiecień 2021

Elektronik

Magazyn elektroniki profesjonalnej

Raspberry Pi 2015

Raspberry Pi

Wykorzystaj wszystkie możliwości wyjątkowego minikomputera

Świat Radio kwiecień - maj 2021

Świat Radio

Magazyn krótkofalowców i amatorów CB

Automatyka Podzespoły Aplikacje kwiecień 2021

Automatyka Podzespoły Aplikacje

Technika i rynek systemów automatyki

Elektronika Praktyczna kwiecień 2021

Elektronika Praktyczna

Międzynarodowy magazyn elektroników konstruktorów

Elektronika dla Wszystkich maj 2021

Elektronika dla Wszystkich

Interesująca elektronika dla pasjonatów