Intro
Zawsze fascynowała mnie możliwość konwersji sygnału czasowego na dane w dziedzinie częstotliwości. Pamiętam jak w wieku około 16 lat dopytywałem mojego ówczesnego nauczyciela matematyki, jak działa transformata Fouriera. Dowiedziałem się wtedy tylko że jest to całka, niestety nic więcej nie zapamiętałem, lub też nie potrafiłem zrozumieć. Ta zagadka nurtowała mnie do tego stopnia, że po wielu latach, nadal wracałem do tematu, natrafiając nareszcie na łopatologiczne wytłumaczenie tego „nadprzyrodzonego” zjawiska jakim jest zamiana oscyloskopu w analizator widma
W owych przedinternetowych czasach, kiedy moim głównym zainteresowaniem była radiotechnika, miałem okazję spotykać się, z różnymi mniej lub bardziej tajemniczymi skrótami, związanymi nierozłącznie z dziedziną jaką byłem ówczas zafascynowany. Niektóre były jasne i oczywiste, inne wydawały się bardzo interesujące w swej niezrozumiałości. Pamiętam to podekscytowanie, kiedy pierwszy raz zrozumiałem działanie pętli synchronizacji fazowej PLL.Nieco później napotkałem na DDS - blisko spokrewnionego kolegę. W powyższych przypadkach budowa tych fundamentów techniki radiowej nie była zbytnio skomplikowana, nie sprawiając większych trudności w dokładnym zrozumieniu zasad ich funkcjonowania. Nieco inaczej sprawa miała się z FFT którego tajniki rozwikłałem tego lata, głównie dzięki doskonałym, intuicyjnym wizualizacjom, oraz przyswojeniu dodatkowych informacji z zakresu matematyki.
Ponieważ moje zamiłowanie do elektroniki z biegiem czasu, nierozerwalnie połączyło się z zamiłowaniem do programowania, poniższy tekst będzie raczej mało merytorycznym opowiadaniem o tym, jak zaprzyjaźniłem się bliżej z działaniem algorytmów będących w dzisiejszych czasach fundamentem technik przetwarzania sygnałów. Najbardziej wartościowym moim z daniem będzie program w którym prócz sterowania MT8880, umieściłem implementację DFT, FFT, IIR oraz filtr Goertzel'a. Kod zamieszczam też, ze względu na rzadkość występowania takich implementacji w 'przyrodzie' (poza kodem elma, nie spotkałem się nigdzie więcej z dft lub fft w asemblerze na avr,).
W dalszej części spróbuję więc opisać w skrócie, jak udało mi się zrozumieć zasadę działania DFT, zapoznając się równocześnie z rozwinięciem tego konceptu, pod postacią FFT. Implementacja drugiego algorytmu, jak już wspomniałem pochodzi z udostępnionych przez Elma Chana bibliotek. Pierwotnie był to tylko asemblerowy moduł do gcc, ja jednak z racji chęci zgłębienia tematu, zaadoptowałem ten kod do bardziej surowego środowiska AS7. Program demo zawiera przełączane w pliku konfiguracyjnym implementacje obu algorytmów z możliwością użycia 'wirtualnego' generatora sinusa do testów, bez potrzeby używania kanału ADC. Platformą sprzętową została płytka 'ATTINY817 Xplained' z dodatkowo podłączonym kwarcem zegarkowym do którego strojony jest generator zegara, aczkolwiek nie jest konieczne używanie kwarcu, po wybraniu wewnętrznego generatora. W kodzie zaimplementowane jest dokładne przeliczanie wartości preskalera, oraz strojenie dla dowolnie wybranej prędkości MASTER w zakresie do 35MHz. Procesor po ustawieniu fusebitu FREQSEL na 20MHz, dał się rozpędzić do takiej wartości zachowując względną stabilność. Można wybrać także dowolną częstotliwości próbkowania adc (wyznaczaną przez timer TCA0). Oraz odczytać rzeczywistą jej wartość, jak też ilość wykonanych pętli głównych na sekundę. Wszystko w to w celu oceny prędkości i poprawności działania programu. Dodatkowo wejście ADC zostało zaopatrzone w wzmacniacz/ogranicznik, aby avr mógł nieco dogonić sprzętowy dekoder pod względem czułości sygnału wejściowego.
DFT (ang. Discrete Fourier Transform)
Jeżeli ktoś próbował by mi kiedykolwiek wytłumaczyć jak działa FFT nie tłumacząc najpierw DFT, to prawdopodobnie, próby takie spełzły by na niepowodzeniu. Warto wiedzieć, że intuicyjny opis Dyskretnej Transformaty Fouriera jest banalny, po uzmysłowieniu sobie kilku kluczowych kwestii. Najważniejszym jest zrozumienie, że to co chcemy zrobić to zwykłe odnalezienie podobieństwa dwóch sygnałów o założonej rozdzielczości i skończonej liczbie próbek (Częstotliwością próbkowania na tym poziomie nie potrzeba sobie nawet głowy zawracać). Druga kluczowa sprawa, to fakt, że w celu wyznaczenia korelacji między przebiegami wystarczy użyć zwykłego mnożenia. Jest to bardzo logiczne, jeżeli rozumiemy jak działa operacja mnożenia.
Do kompletu potrzebujemy jeszcze tylko zsumować wyniki tych jednostkowych mnożeń w jakimś akumulatorze i można powiedzieć że mamy prawie to o co nam chodziło. Dlaczego prawie, o tym zaraz się dowiemy. W mojej implementacji DFT założyłem zbudowanie generatora sinusoidalnego o zadanej częstotliwości, którego próbki, posłużą następnie jako źródło owego sygnału wzorcowego z którym przemnożymy sygnał badany, w celu określenia owego stopnia korelacji. Generator ten jest opisany na rysunku jako NCO (ang. numerically controlled oscillator), funkcjonalnie jest to moduł identyczny z tym co znajduje się w syntezerach częstotliwości DDS. W moim rozumieniu tematyki, różnica pomiędzy DDS i NCO to brak przetwornika cyfrowo-analogowego w przypadku tego drugiego. Wracając do metody pomiaru korelacji stanowiącego podstawę obliczeń. Warto wiedzieć, że dane wyjściowe w sumatorze uzyskane w taki sposób będą mocno zależne od różnicy fazy pomiędzy mnożonymi sygnałami. Z tego powodu w celu uniezależnienia się od przesunięć fazowych, do mnożenia używane są dwa wzorcowe przebiegi sinusoidalne, przesunięte jednak między sobą o 90 stopni. W tym przypadku uzyskujemy dwie dodatnie wartości wyjściowe, o poziomach zależnych od przesunięcia fazowego sygnału badanego względem wzorca. W celu ostatecznego określenia korelacji, należy jedynie dokonać spierwiastkowania sumy kwadratów tych składowych (części rzeczywistej i urojonej
Ze względu na czasochłonność pierwiastkowania, początkowo w programie eksperymentowałem z uproszczoną aproksymacją, ostatecznie korzystając z dwójkowego logarytmu, który w tym przypadku działa jak całkiem dobre przybliżenie pierwiastka 2 stopnia (na wykresach funkcji sqr i log wyglądają podobnie, ale leżą na przeciwnych osiach współrzędnych).
Program zbudowany w ten, sposób umożliwia dowolny wybór zadanych punktów „badania” - częstotliwości (bins). Przy czym ilość próbek akumulowanych dla uzyskania wyniku jest odwrotnie proporcjonalna do wymaganej szerokości widmowej takiego punktu (nieco odwrotna analogia do dobroci filtru,). Interesujące w przypadku DFT jest skromne zapotrzebowanie na pamięć, gdyż nie jest używany bufor próbek, obliczenia stopnia korelacji są dokonywane niejako „w locie”. Umożliwia to uzyskanie wyniku pomiaru z ramki, która wykracza poza zakres dostępnej pamięci ram, jak też bardzo wąskich punktów, których obliczenie na procesorze z małą ilością pamięci, w przypadku użycia typowego FFT było by niemożliwe.
Podczas omawiania FT, nie można nie wspomnieć o zjawisku wycieku spektralnego (z ang. 'spectral leakage' – jakże malownicze określenie:) ). Polegającym na psuciu wynikowej charakterystyki przez nieciągłość próbkowanego sygnału. Zjawisko powstaje min. w przypadku, gdy sygnał próbkowany, nie zostanie ujęty równo w okresowości jego występowania (np. ucięty w pół sinusoidy, a nie w miejscu przejścia sygnału przez zero). Remedium na to niepożądane zjawisko, jest użycie tak zwanego 'oknowania' (ang. windowing). Zabieg polega na wstępnym przemnożeniu ramki wejściowej z tablicą okna o kształcie zależnym od pożądanej aproksymacji. Zazwyczaj początek i koniec dyskretnych wartości tej tablicy leży w pobliżu zera, jednak są wyjątki. Stosowanych jest wiele kształtów, różniących się właściwościami filtrującymi. Wstyd się przyznawać, ale jeszcze niedawno nie wiedziałem, co dokładnie oznaczają owe tajemnicze nazwy tj. min: „hanning, hamming, blackman, flat-top” oraz inne, które wpływają na pewne, czasem bardzo subtelne zmiany w wyglądzie obrazowanego sygnału. Teraz już wiem, że chodzi o kształt owego okna zabezpieczającego wykres widma przed wyciekiem
Program liczący DFT skonstruowany jest w taki sposób, aby możliwe było wybranie rozmiaru próbkowanej ramki, oraz ilości i częstotliwości obliczanych punktów. Rzeczywista prędkość obliczeń jest niewielka i dla zegara 4MHz umożliwia pomiar zaledwie 8 punktów 256 próbkowej transformaty, przy częstotliwości próbkowania 4kHz.
W celach testowych możliwe jest użycie emulowanego sygnału wejściowego (dodatkowy generator NCO, zamiast próbek z adc). Oprócz surowych danych wyjściowych w postaci wartości uint16, rysowany jest prosty wykres składający się z kolejno wybranych punktów częstotliwości. Rozmiar w osi pionowej to 8 punktów (logarytm z 8 bitowej wartości).
FFT (ang. Fast Fourier transform)
Wszystkie powyżej opisane działania od momentu oknowania, możliwe są do wykonania w sposób o wiele szybszy, jednakże angażujący znacznie większą ilość pamięci RAM. W przypadku FFT operujemy na buforze, który jest wprost proporcjonalny do rozmiaru obliczanej jednorazowo ramki danych. W implementacji zaczerpniętej od Elm-a, mamy do czynienia z 16 bitową transformacją na liczbach całkowitych. Dla porządku wspomnę tylko, że tego typu tranformata wymaga ramki danych o rozmiarze będącym potęgą dwójki, oraz bufora próbek o rozmiarze dwa razy większym niż normalnie używany do zapamiętania przebiegu. Jest tak ze względu na fakt, iż FFT będzie używać specyficznego algorytmu, w którym połowa próbek wejściowych to część rzeczywista, druga połowa natomiast reprezentuje urojoną część obrabianego sygnału. W naszym przykładzie próbki dla części urojonej, są w momencie zapisu tego bufora zwykłymi kopiami części rzeczywistej. Proces ten jest wykonywany w module 'FFT INPUT', w momencie otrzymania próbki, zwykle z przetwornika ADC. Po zapisaniu bufora oczekiwaną liczbą próbek, na owym buforze zostanie dokonana transformacja za pomocą modułu 'FFT EXECUTE'. Po przejściu przez „algorytm motylkowy” bufor w którym znajdowały się próbki w domenie czasowej, zostaje zapełniony danymi wynikowymi, przedstawiającymi rozmiar prążków w domenie częstotliwościowej. Jednak będą one ułożone w kolejności wymagającej dodatkowego uporządkowania do postaci liniowej w funkcji częstotliwości, co skolei dokona się w ostatnim module 'FFT OUTPUT', wraz z równoczesnym wysłaniem danych wynikowych na konsole w celu ich wyświetlenia.
Na rysunku zaznaczone są zatem 3 główne części przetwarzania danych naszego algorytmu. Pierwsza o nazwie 'FFT INPUT' zawiera procedury pobierania próbek, mnożenia ich z oknem anty wyciekowym, oraz umieszczania ich w buforze danych. Drugi moduł to niejako serce szybkiej tranformaty Fouriera 'FFT EXECUTE'. Zajmuje się mnożeniem i dodawaniem danych z bufora z tablicami sinusa i cosinusa znajdującymi się w pamięci ROM. Operacje matematyczne przebiegają w ciekawy sposób, co w przypadku bardziej obrazowego przedstawienia metody, przypomina charakterystyczne motylki. Na koniec tak zapisane próbki należy odczytać w specyficznej kolejności. Odpowiada za to moduł 'FFT OUTPUT', który ze względu na ograniczenia pamięci małego procesora, podczas pobierania danych wyjściowych, jest na bieżąco angażowany w celu wyświetlenia wyników w prawidłowej kolejności. W innym razie należało by użyć dodatkowego bufora.
Podobnie jak w przypadku DFT, program liczący FFT umożliwia obsługę rozmiaru próbkowanej ramki, oraz częstotliwości samplowania. Jednakże obsługiwanymi przez procesor z 512B pamięci, są jedynie najmniejsze rozmiary ramki o wartościach 64 oraz 32. Odpowiada to 16 oraz 32 punktom wartości wynikowych. Z DFT wspólne jest też końcowe sumowanie części rzeczywistej i urojonej, co w tym przypadku jest robione w momencie pobierania próbek z bufora w kolejności podyktowanej tablicą 'tbl_bitrev'. Tym razem obliczenia dokonywane są błyskawicznie, bez problemu osiągana jest częstotliwość próbkowania równa 48kHz dla 16MHz zegara.
Podsumowując na zakończenie. Gdybym miał dzisiaj powiedzieć w kilku prostych słowach, czym czym jest owa transformat królowa (w powyżej przedstawionej realizacji, bo to nie jedyny sposób jej obliczania):
FFT jest dekompozycją/transformacją wejściowego, dyskretnie próbkowanego sygnału i przeniesienie go z pierwotnej domeny czasowej, do domeny częstotliwościowej, z użyciem tablic sinus/cosinus oraz sprytnego motylkowego algorytmu zwieńczonego zmianą kolejności próbek
Filtr Goertzel'a
Na zakończenie eksperymentów, postanowiłem napisać implementację dekodera DTMF z użyciem owego filtra. Okazał się on bardzo zbliżony do filtra z nieskończoną odpowiedzią impulsową, więc od budowy takiegoż właśnie zacząłem. Gdy już zrobiłem najzwyklejszy biquadowy filtr IIR w postaci 2, wystarczyło go nieco zmodyfikować by mógł posłużyć za wąski filtr pasmowy do badania składowych widma. W celu pomiaru wartości, program został wyposażony w miernik wartości szczytowej sygnału, oraz licznik, który resetuje pamięć próbek, oraz 'peak hold' równocześnie zapamiętując wartość wynikową. Program Gfilter operuje na 8bit próbkach wejściowych i 16bit próbkach w samym filtrze, współczynniki mnożenia są również 8 bitowe, przy czym max wartość wzmocnienia w sprzężeniach wynosi x4 (taki FMULSU x2). Niska rozdzielczość wymusza konieczność dobrania wartości, które nie spowodują zbyt szybkiego nasycenia się 16b rejestrów. W mojej implementacji usunąłem bezpośrednie połączenie sygnału wejściowego z wyjściowym, w celu lepszej obserwacji danych będących wynikiem działania algorytmu. Optymalna ilość obliczanych próbek wynosi około 200, co umożliwia podobnie jak w przypadku implementacji z DFT użycie ”filtracji większościowej” co poprawia odporność na przypadkowe interpretacje znaków, kosztem dłuższego czasu analizy sygnału. W przypadku programu Gfilter, użyta została częstotliwość próbkowania o wartości 8kHz, dla której czas jednej ramki pomiarowej wynosi 12.5ms. Tak krótki czas pomiaru widma umożliwia bezproblemową detekcje sygnałów o długości 50ms nawet po użyciu wspomnianego wcześniej uśredniania. W przypadku detektora dtmf w oparciu o filtr byłem nieco rozczarowany faktem, iż ilość mocy obliczeniowej potrzebnej do uzyskania należycie niezawodnych wyników była porównywalna do implementacji z DFT, jednak reakcja dekodera w tym przypadku była o wiele szybsza (można by dekodować nawet 25ms sygnały).
G8880 vs DFT
Po wyprodukowaniu powyższego tekstu, wraz z programem 'demo'. Uzmysłowiłem sobie praktyczną jego bezużyteczność. Co spowodowało w krótkim czasie dopisanie małego programiku zajmującego się kodowaniem i dekodowaniem pary tonów w telefonicznym standardzie DTMF. Implementacja DFT aż prosiła się, by użyć jej w celu dekodowania 16 par tonów, jak również w celu ich generowania – właściwie wszystko do tego było już gotowe. Program 'pół demo' na ten moment na porcie A6 generuje za pomocą 8 bitowego daca kolejno 16 kodów, które są w tym samym czasie dekodowane przez adc na porcie A5. Efekt działania tego można zobaczyć na załączonych filmach. Ciekawy byłem, z jak bardzo zniekształconym i zaszumionym sygnałem algorytm detekcji potrafi sobie poradzić. W Dalszej części eksperymentu postanowiłem uruchomić zabytkowy transceiver dtmf w ceramicznej obudowie i zobaczyć, jaka będzie różnica w działaniu dekodera sprzętowego w porównaniu z programową 8bitową implementacją. Uwagę zwracałem głównie na niezawodność detekcji (minimalny wymagany poziom sygnał/szum), podatność na przekłamania, czas reakcji po podaniu poprawnego tonu. G8880 był sterowany tą samą płytką z procesorem, podłączony do niej poprzez złącze goldpin.
TX
W przypadku generatora par tonów, dla G8880 widać wyraźne zaszumienie sygnału wyjściowego, ograniczające dynamikę takiego sygnału do około 40dB. Poziom zniekształceń jest również niewielki. Aby uzyskać zbliżone parametry sygnału na wyjściu daca procesora (za pojedynczym filtrem RC), należało użyć częstotliwości próbkowania rzędu 22kHz. Po zoptymalizowaniu procedur wystarczający byłby zegar rzędu 4MHz. Ze względu na to, że w programie demo pracuje cały czas 8 generatorów (równoczesne generowanie i dekodowanie). Wymagana jest dużo większa prędkość do sprostania zadaniu. W temacie nadajników to w zasadzie tyle. Noty katalogowe nie wspominają dokładnie o tym, w jaki sposób jest tworzony sygnał w scalaku. Procesor natomiast korzysta z 8bitowych tablic sinusa, które są sumowane w 8b sygnał wyjściowy dostępny na wyjściu daca (port A6).
RX
W przypadku odbiorników jest już szersze pole do eksperymentowania. W celu uzyskania optymalnych rezultatów, warto użyć odpowiednio niskiej częstotliwości próbkowania oraz odpowiedniego rozmiaru ilości próbek dla DFT. Należy pamiętać, że zwiększając częstotliwość próbkowania, zmniejszamy rozdzielczość widmową transformacji, co z kolei wymusza zwiększenie ilości zebranych próbek dla zachowania możliwości rozróżniania zadanych częstotliwości. Tak więc, również ze względu na pobór mocy, częstotliwość próbkowania warto wybrać taką która tylko z małym zapasem zapewnia poprawne samplowanie zakresu dekodowanych tonów. Umożliwi to użycie stosunkowo niskiej ilości próbek transformaty. SR=6kHz oraz 128 sampli jest wystarczające do prawidłowego działania tak skonstruowanego dekodera. Czas jednej transformacji wynosi 6kHz / 128 = 21ms. Z taką teoretyczną prędkością można by zmieniać tony dla poprawnego ich rozpoznawania. W praktyce jednak bardziej zależy nam na nie występowaniu „fałszywych alarmów”, a te są nieuniknione w przypadku braku 'uśredniania' zdekodowanych znaków. Po doprowadzeniu sygnału szumu białego o wysokiej amplitudzie, dekoder zwraca przypadkowo zinterpretowane znaki. Zmniejszenie amplitudy sygnału, przekłada się na proporcjonalne zmniejszenie występowania fałszywych interpretacji szumu jako sygnał użyteczny. Dopiero po zastosowaniu sprawdzania kilku kolejnych takich samych znaków (minimum 3) sytuacja przypomina bardziej to, z czym mamy do czynienia na wyjściu dekodera sprzętowego, czyli brak znaków wyjściowych w reakcji na sygnał szumu. W tym miejscu jednak spotykamy się z ograniczeniem możliwości filtracji wyjściowej w kontekście wymaganej minimalnej długości czasu sygnału wejściowego. Podane wyżej parametry SR oraz ilości próbek, są właściwie wartością graniczną dla 50ms czasu tonu wejściowego (dla pełnej czułości dekodera należałoby zwiększyć SR do 8KHz co daje czas 3 transformacji na poziomie 48ms). Uzyskanie braku fałszywych interpretacji szumu jest możliwe dopiero po zmniejszeniu czułości dekodera (parametry w programie), co spowodowało, że rozkład mocy szumu, nigdy nie osiągnie progów detekcji, jednak dla poprawnych sygnałów dtmf pozostanie zapas kilku dB. Byćmoże z wyżej wymienionych przyczyn do celów dekodowania dtmf w procesorach sygnałowych używane są raczej algorytmy, oparte o filtrację sygnału, np. algorytm Goertzela. Być może są też algorytmy działające podobnie do układu MT8880 (w skrócie dwa strome filtry dla grup tonów wysokich i niskich, oraz liczniki przejścia przez zero fali w celu pomiaru częstotliwości składowych). Bardzo ciekawy przykład dekodowania dtmf oparty jedynie o czasy pomiędzy przejściami przez zero znalazłem w ty miejscu:http://www.romanblack.com/DTMF/DTMF_alg.htm
Tak jak się spodziewałem implementacja dekodera za pomocą DFT miała przewagę rozpoznawania kodów w przypadku bardzo małego stosunku sygnał szum, jednak różnica nie była jakaś znacząca w praktycznym zastosowaniu. Sprzętowy dekoder wypadał zncznie lepiej w zakresie nadinterpretacji szumu. Mimo dużego zakresu dynamiki nie zauważyłem aby kiedykolwiek spreparowany w przykładach szum został zinterpretowany jako sygnał dtmf, mimo bezproblemowej detekcji 50ms kodów. Na filmach chciałem pokazać działanie wizualizacji FFT/DFT oraz zachowanie dekodera dla różnych implementacji filtrów jak i sprzętowego dekodera. Niestety przygotowane krótkie screencasty nie zostały zaakceptowane przez elektrodowy player, więc zamieszczam je w załączniku. Na stronie tylko zajawka w jakiejś marnej rozdzielczości i klatkarzu.
Fajne? Ranking DIY