Rychlá Fourierova transformace na AVR

Rychlá Fourierova transformace (FFT) slouží k převodu signálů z časové do frekvenční oblasti, inverzní (zpětná) Fourierova fransformace (IFFT) slouží k opačnému převodu. Jedná se o rychlou verzi diskrétní Fourierovy transformace (DFT). FFT je převážně doménou PCček a signálových procesorů (DSP), které jsou pro tento druh výpočtu patřičně výkonově vybavené. Možná proto si málokdo dokáže představit real-time výpočet FFT na 8-bitovém jednočipovém mikropočítači. Dnešní jednočipy ale dosahují výpočetních výkonů i přes 20 MIPS a proto tato představa není nereálná. Já jsem se ze zvědavosti pokusila tuto představu realizovat a výsledek mě pozitivně překvapil. Napsaný software zvládá real-time výpočet FFT a IFFT řádu 6 (64 bodů) na procesoru ATmega32 16MHz při vzorkovací frekvenci 3200Hz, data získaná FFT zobrazuje na LCD a ještě mu zbývá přes 40% výpočetního výkonu na případné další operace. Řád FFT je omezen velikostí vnitřní SRAM (2KB), na procesorech ATmega8 a ATmega16 lze se stejnými výsledky realizovat výpočet FFT řádu 5 (mají jen 1KB SRAM), na procesorech ATmega64 a 128 (4KB SRAM) by bylo možno počítat FFT až do řádu 7 (128 bodů).

Výpočet je prováděn v pevné řádové čárce, čísla jsou ve formátu 16.16 (16 bitů celá část, 16 bitů desetinná část). Jedno komplexní číslo tedy v paměti zabírá 8 bytů. Program je optimalizován na rychlost a přesnost, je při něm využito všech 32 vnitřních registrů procesoru. Délka celého programu včetně lookup tabulek s twiddle faktory je 2560 bytů. Přesnost výpočtu je velmi dobrá (pro některé aplikace možná až zbytečně), ve srovnání s MATLABem se vypočtené hodnoty liší až na čtvrtém desetinném místě. Pro vzorkování a výpočet je použita dvojice bufferů, obsah jednoho bufferu se vysílá ven přes PWM a zároveň plní novými hodnotami z A/D převodníku, paralelně s tím se obsah druhého bufferu zpracovává. Po provedení všech výpočtů a naplnění bufferu, do kterého se vzorkuje, se buffery prohodí. Vstupní signál se tedy nepřetržitě vzorkuje, fransformuje, zobrazuje, zpětně transformuje a převádí zpět na analogovou hodnotu. Parametry výpočtu jsou následující:

Parametry transformace na ATmega32, 16MHz:

Algoritmus výpočtu: Cooley-Tukey s decimací ve frekvenci (DIF), Radix-2

Počet bodů FFT (IFFT): 64

Vzorkovací frekvence: 3200 Hz

Rozlišení ve frekvenci: 50 Hz

Zobrazení: bargraf 16 sloupců na 2*16 znakovém LCD (zobrazuje se 16 dolních spektrálních čar)

Doba vzorkování: 20 ms (doba vzorkování 64 vzorků)

Doba výpočtu FFT: 4.04 ms

Doba výpočtu IFFT: 4.16 ms

Doba řazení dat (bitinverze): 222 us (provádí se dvakrát)

Doba zobrazování na LCD: 2.31 ms (včetně výpočtu absolutních hodnot)

Zbývá: cca 9ms

Možná oblast použití je široká, od efektního zobrazování spektra na domácím audio zesilovači, přes harmonickou analýzu sítě (do 32. harmonické) až po různé typy modulací, filtrací a kódování.

Následujících několik fotografií zachycuje spektra různých signálů v podobě sloupcového grafu, tak jak jsou vypočtena realizovaným algoritmem FFT. Zobrazuje se celkem 16 frekvencí od 0Hz do 750Hz s krokem 50Hz. Ač to není z obrázků patrné, výpočet probíhá v reálném čase a údaj na displeji se obnovuje každých 20 ms. Protože by byl při analýze neharmonických signálů (např. hudby) pohyb spektrálních čar pro lidské oko příliš rychlý a nepříjemný, je pohyb sloupců směrem dolů zpomalen, někdy bývá tento režim označován jako "slow decay".

Spektrum sinusového signálu 50 Hz
Spektrum sinusového signálu 80 Hz
Spektrum sinusového signálu 100 Hz
Spektrum sinusového signálu 230 Hz
Spektrum obdélníkového signálu 50 Hz
Takto vypadá spektrum písničky Rain od Madonny někde na začátku druhé minuty :)

Na následujících oscilogramech jsou zachyceny průběhy signálů po A/D převodu, provedení dopředné a zpětné Fourierovy transformace a následném převedení zpět do analogové podoby. K D/A převodu je použita 8-bitová PWM 31,25 kHz s aktivním filtrem typu Butterworth 2. řádu se zlomovou frekvencí 1600 Hz. Z obrázků není zcela zřejmá doba zpoždění mezi vstupním a výstupním signálem, jedná se od dvě délky vzorkovacích bufferů tj. o 40ms.

Stopa 1 (horní): sinusový signál 50 Hz před transformací Stopa 2 (dolní): sinusový signál 50 Hz po průchodu řetězcem ADC - FFT - IFFT - DAC (PWM) - filtr
Stopa 1 (hladká): sinusový signál 200 Hz před transformací Stopa 2 (zubatá): sinusový signál 200 Hz po průchodu řetězcem ADC - FFT - IFFT - DAC (PWM) - filtr