DFT spektrumanalizátor – nullverzió

Ne legyenek illúziók, eddig se voltak kész projektek, és ez se lesz az…

DFT spektrum analizátor
DFT spektrum analizátor

De azért ez a projekt is elindult és van elkövetett működő első példány, nullverzió, prototípus, ahogy tetszik. Igazából amiért érdemesnek látom ebben a kezdetleges fázisban megosztani, az az egyszerűségének tudható be. Ugyan lehetne sok ponton fejleszteni-bővíteni, de akkor átláthatatlanabb is lesz, most pedig még letisztult, könnyen átlátható és így könnyebb tanulmányozni is.


Rövid videó a működésről

Lényegében három funkciót lát el, és eléggé lazán van ez a három funkció időzítve, egymás után kerülnek végrehajtásra, nincs átfedés, párhuzamosítás. Ebből fakadóan megszakításkezelésre se lenne szükség, de a mintavételező mégis azzal működik, egyszerűen mert a korábbi munkákból így hoztam át és miért variáljam meg azt, ami már egyszer jól dolgozik. A három egymást követő funkció tehát a következő:

  1. Mintavételezés (N db analóg minta bemintavételezése egy RAM pufferbe)
  2. Fourier-transzformáció (csúcsdetektálással egybekötve kimeneti pufferbe dolgozva)
  3. Kijelzés (oszlopok felrajzolása a DFT kimeneti pufferéből és az oszlopok visszahullásának vezérlése)

A mintavételező

N db mintát vesz az A/D konverterrel Fs mintavételi frekvencián, ami jelen esetben N=32 elemet jelent és Fs=12kHz mintavételi frekvenciával. A minták 8 bitesek, és egy 32 bájtos pufferbe íródnak be a RAM-ba. Amikor összegyűlt a 32 minta, akkor leáll a mintavételezés (kikapcsoljuk a megszakításokat) és a Fourier-transzformációra fut a végrehajtás. A mintavételezés alatt a megszakításkezelés engedélyezve van, miközben a főprogram áll egy önmagára ugró brts utasításon. A mintavételezés befejezését a T jelzőbittel jelezzük a főprogramnak. (Tehát a T a mintavételezés folyamata alatt 1, és 0-ra váltással lépünk ki belőle, pont fordítva, ahogy a korábbi projektekben csináltam)

Azért az érezhető, hogy amíg megy a mintavételezés, addig közben mást is csinálhatna a főprogram, minthogy egy önmagára ugró utasításon álljon, azaz itt el van vesztegetve nem kevés kihasználatlan számítási kapacitás. Mivel azonban itt most az egyszerű alapműködés elérése volt a cél, így a párhuzamosításokat inkább mellőztem.

A Fourier-transzformáció (DFT)

Csak nagyon röviden programozás szemszögből nézzük a dolgot, nem megyünk bele matematikai részletekbe. A mintaszám tehát N=32, ebből 16 spektrumvonal számítható ki, pont ennyit tudunk az LCD kijelzőn megjeleníteni. Az így létrejött spektrum lineáris felbontású, ezért túl szép nem lesz, mert az audio spektrumanalizátornál logaritmikus frekvencialéptékezés kellene. (ez az ami a korábbi IIR szűrős spektrumanalizátornál könnyen megoldható volt) Sajnos ez mindig is hátrány lesz az FFT (DFT) alapú spektrumanalizátorok esetében, bár sokan sokféleképpen próbálnak ezen javítani. Pl. kellően nagy N elemszám választása után átszámolhatóak a spektrumvonalak logaritmus léptékezésre, vagy lehet 2-3 lineáris szakaszban működő több különböző mintavételezés és transzformáció, ahol pl van egy 100-1000Hz közötti szakasz 100,200,300…800Hz spektrumvonalakkal, majd jön egy 1k-8k szakasz kHz-es egységekkel. Természetesen ehhez elő kell állítani egy kisebb mintavételezésű lekvantált időtartományt is, amit pl. egyik neten keringő változat esetében külső analóg szűrővel oldottak meg, de ezt szoftveresen digit jelfeldolgozással is meg lehetne csinálni. Ezek a lehetőségek egyelőre félre vannak téve a későbbi továbbfejlesztésekre. Mintavételi frekvenciának 12kHz-et választva (pontos értéke 12019.23Hz) 6kHz a felső határ, a létrejött spektrumvonalak pedig: 376, 751, 1127, 1502, 1878, 2254, 2629, 3005, 3380, 3756, 4132, 4507, 4883, 5258, 5634, 6010 Hz. A spektrum közepének belőtt 1kHz a 16 oszlopból a harmadik, vagyis jól látjuk, hogy nem a megszokott log karakterisztikát követi, és 376Hz-es első spektrumvonalból azt is látjuk, hogy az igazán mély hangok már kívül esnek a mérésen, mint ahogy a 6kHz feletti igazán magas hangok is.

A diszkrét Fourier-transzformáció lényege, hogy veszünk egy adott mintaszámú diszkrét időtartományú bemeneti mintát (tömböt) és ezt végigszorozgatjuk a pontosan illeszkedő szinuszos és koszinuszos harmonikusokkal elemről elemre haladva, miközben ezeket a szorzateredményeket összegezzük (akkumuláljuk). A művelet végén minden harmonikusra kapunk egy szinuszos és egy koszinuszos összetevőnk, tehát két N elemű tömbünk keletkezik. (Amennyiben komplex algebrában gondolkodunk, úgy csak egy N elemű komplex tömbünk, de itt most külön van a komplex érték valós és képzetes része) N elemű időtartományból N elemű frekvenciatartományt kapunk, mely N/2-re szimmetrikus, N/2 felett a spektrum tükörképét kapjuk, amire nincs szükség, így a spektrumot csak N/2-ig számoljuk ki. N=32 elem esetén csak az alsó 16 spektrumvonal kell. Basic nyelvjárásban valahogy így néz ki egy DFT algoritmus:

LET N=32   ; elemszám
DIM X[N]   ; bemeneti tömb (időtartomány)
DIM YS[N]  ; kimeneti tömb a szinusz összetevőkkel
DIM YC[N]  ; kimeneti tömb a koszinusz összetevőkkel
FOR H=1 TO N/2 ; harmonikusok ciklus
  FOR T=1 TO N ; időtartomány ciklus
    YS[H] += X[T] * -SIN(2*PI*H*(T-1)/N) ; szinusz harmonikussal szorzás
    YC[H] += X[T] *  COS(2*PI*H*(T-1)/N) ; koszinusz harmonikussal szorzás
  NEXT T
NEXT H

A fenti kódban a tömb indexelése 1-től indul, a szinusz és koszinusz függvényekben azonban 0-tól kell indulnia, ezért látunk ott (T-1) kifejezést a T helyett. Nyilván más nyelvjárásban más lehet a szitu. Lényegében a (T-1)/N alak egy 0-1 közötti értéket ad, ez az időrész, H a harmonikus, 2*PI a teljes kör radián. Ha egyben nézzük a 2*PI*H*(T-1)/N részt, akkor sorra előállítódik az adott H-adik harmonikus, azaz az első külső ciklusban a belső ciklus leír egy teljes szinuszgörbét, a második külső ciklusban a belső ciklus pontosan két teljes szinuszgörbét és így tovább. YS[] tömbbe a szinusz, YC[] tömbbe a koszinusz harmonikusokat összegezzük be, azonban a szinuszos összetevőt negatív előjellel kell venni. Amikor egyszer végigfut a belső ciklus, mindig létrejön egy harmonikus, így a külső ciklusban sorra létrejön az első, második stb harmonikusok szinusz és koszinusz összetevői. Természetesen mi nem ezeket az összetevőket akarjuk kigrafikonozni, hanem az abszolút nagyságot. A két összetevő abszolútértékét pedig a Pitagorasz tétellel tudjuk számítani:

√(YS2 + YC2)

Azaz a magnitúdó (abszolút nagyság) számításához négyzetre kell emelni a szinuszos és koszinuszos összetevőt, össze kell adni őket, majd gyököt kell vonni. Ezután ez az érték kijelezhet akár lineáris, akár logaritmikus (dB-es) skálán. Mivel azonban a kijelzés komparátorozással történik, így a gyökvonást el is lehet hagyni, a nyers négyzetes értékre is készíthetünk 16 komparálási műveletet. A gyökvonás amúgy is mumus, nincs rá hardveres végrehajtás, vagy le kell programozni az algoritmusát, ami több-10 órajel, vagy táblázatot kell csinálni és onnan kiolvasni. (A sebesség céljából ilyenkor elég gyakran használunk macerásabb függvények esetén értéktáblázatot, bár ezek néha kényelmetlenül nagy adathalmazok is lehetnek, amik esetleg már nem férnek be a RAM-ba sem) Ha szükséges, akkor természetesen a fázis is számítható árkusztangenssel, amit szintén táblázatból kiolvashatóra lehet csinálni, de ez most nem játszik. (esetleg egy későbbi grafikus LCD-re dolgozó komolyabb, mérésre szánt DFT esetében – pl. van elképzelésem digit szkópra dsPIC-el…)

Igazából ha jól megnézzük, semmi extra nincs, gyerekjáték még assemblyben is a lekódolása. Egy dologra kell csak nagyon (mondhatni rohadtul) figyelni, hogy hol mi mennyi biten jön létre. Piszokul el lehet vele szállni, engem is szanaszét szívatott eleinte. Pontosan ki kell számítani, hogy az összegzőregiszterek max határa hol lesz. Én pl először 24 bites regisztereket hoztam létre az összegzőnek, amibe amúgy tuti befér az összegzett eredmény. (8 bites bemeneti mintát szorzunk 8 bites együtthatóval, ami 16 bit, és ezt 32-szer összeadjuk, ami +5 bit, azaz 21 bitnél soha semmilyen körülmények között nem lehet az összegzés eredménye. Valójában azonban még ennyi se, mert az együtthatók szinusz függvények, melyeknek van középértékük, és inkább azzal kellene számolni. Továbbá az ablakfüggvénynek is van középértéke, azt is bele kellene számolni, és ehhez jön, hogy ez is csak morbid, pl. négyszögjelnél az első harmonikusra ad baromi nagy összegértéket. Szinusz sweepre már kevesebbet, fehérzajra meg pláne még kevesebbet. Végül olyan megoldás született, hogy az összegző regiszter 16 bites lett, de az előtte lévő szorzás művelet 16 bitről csak a felső 8-ra lett csonkolva. Mivel ez csak egy spektrumanalizátor, aminek csak kb 20dB dinamikatartomány kell, így ezzel is jól működik, és a számításigény sokkal kevesebb. Mellesleg megjegyzendő még, hogy előjeles számokkal van dolgunk, és kicsit nehezebb pl túlcsordulást lekezelni, mintha csak egészek lennének.

Mértem amúgy a DFT futásidőt, kb 14 ezer cpu órajel (egy korábbi 24 bites akkumuláló regiszterekkel), ami kb 0.7ms. Ezzel szemben a 32 minta begyűjtése kb 57 ezer cpu órajel (kb 2.6ms), így látható, hogy a transzformáció párhuzamosítható a mintavételezéssel, akár röptében trafózással, akár úgy, hogy dupla puffereléssel mindig a korábbi ciklus trafója folyik miközben a következő trafózás adatait mintavételezzük, ezzel folyamatos mérést lehetne megvalósítani a jelenlegi szakaszos helyett. Ez is a fejlesztési célok egyike.

Ablakfüggvény

Fourier-trafó esetében megkerülhetetlen az ablakolás, azaz hogy a time window szélein leszintezzük a jelet. Ezt egy haranggörbe alakú erősítésvezérléssel oldják meg, ami annyit tesz, hogy az időtartománybeli jelet végig kell szorozni egy ilyen függvénnyel. Ezt a szorzást én beépítettem az együtthatókba, így kevesebbet kell szorozni a Fourier-trafó alatt. Az ablakolás amúgy azért kell, mert a trafó valójában egy periodikus jel spektrumát számolja ki, mely a mintavételezett jelünk végtelenített egymásutánban játszásából alakulna, így a jel része lenne az ablak végének és a következő periódusnak számítva az elejének a találkozási pontja, ami hamis főleg magasfrekvenciás spektrumvonalakat eredményezne. Az ablakfüggvény megtisztítja ettől a spektrumot, azonban ilyenkor a spektrum első 1-2 spektrumvonala sérül egy kicsit. Sokféle ablakfüggvény létezik, én amolyan Hamming féleséget választottam, a széleken nem zérusra csökkentve a jelet, hanem csak a 20 százalékára. Az ablakfüggvény formáját egy koszinusz görbe adja. Az így létrehozott harmonikusok együtthatóit egy C programmal generáltam le PC-n, a mikróvezérlőben pedig  a programmemóriába kerül, és onnan lesz inkrementált címzéssel végigolvasva a transzformáció alatt.

Kijelzés

Tehát 2.6ms-ig tart a mérés ami amúgy pont a legelső spektrumvonal periódusideje. Azt is lehet mondani, hogy 376Hz-en ismétlődne csak a mérés. Persze a transzformáció ideje is hozzájön, meg igazából a kijelzőkezelés ideje is, utóbbinál a sok kényszerű delay függvényt is nézni kell. Kb 200-250Hz környékén pörögne a teljes folyamat és ez lenne a kijelzőre kiíratás frissítése is. Ez lelassítható 30-60Hz környékére, de ha nem teszünk be lassítás akkor is működik rendesen, abból nincs gond, hogy gyorsabban irkáljuk az LCD vezérlő RAM-ját, mint ahogy az hardveresen frissít. Én is több variációt próbáltam végig, mint pl. 16 bites hardveres számláló átfordulásának figyelését, vagy egyszerű számlálót (a megosztott kódban ez működik jelenleg). Ha semmilyen kijelzéslassítás nem volt, akkor is működik, csak kicsit gyorsan esnek vissza az oszlopok, és a visszahullás kezelésnél nehéz lenne megoldani a vezérlést. Maga a kijelzőkezelés első felindulásból 16 egymás követő komparálási művelet volt, mely oszlopról-oszlopra ment, és rajzolta meg az oszlop alsó és felső részét (a két soros LCD-n ugyebár) Na a sok delay miatt ezzel alaposan lelassult a main, amit eleinte nem értettem, csak menet közben esett le, hogy a kijelzés a belassulás okozója. Most soronként írja ki a kijelzőt, és inkább kétszer olvassa ki az oszlopok értékeit a RAM-ból, de így legalább nem kell 32-szer kurzorpozicionálás parancsot küldeni. Másfelől a komparálási sorozat helyett gyorsabb, táblázatos módszert találtam ki. 16 vagyis 0-val együtt 17 értéket tudunk megjeleníteni egy oszlopban. Karakterhelyenként 9 lehetőség van. Kódoljuk be egy-egy bájtba a kiíratás képét, azaz hogy az alsó és a felső helyre milyen érték megy 0-8 között, és készítsünk egy olyan 256 elemű tömböt, amit ha megindexelünk a spektrumvonal értékével, akkor a kiolvasott bájtból azonnal kivehető, hogy milyen ASCII kódot kell a kijelzőre küldeni. Ezzel a megoldással óriásit gyorsult a kijelzés szakasz. A kiolvasott bájt alsó 4 bitje 0-8 között az alsó karakterhely, a felső 4 bit a felső karakterhely értékét tartalmazza. 0 esetén szóközt kell írni, minden más esetben a kiolvasott érték mínusz egyet. Végül fontos még megjegyezni, hogy a DFT algoritmusban a spektrum pufferbe úgy íródik be az adott spektrumvonal nagysága, hogy ha nagyobb, mint amit a puffer tartalmaz, akkor felülírja, egyébként eldobja, azaz csúcsdetektálással építi a spektrumot. Azonban ezzel beragad a kijelzőre a legnagyobb érték, tehát gondoskodni kell róla, hogy a pufferben tárolt spektrumvonalak értékei folyamatosan csökkenjenek, ennek vezérlése pedig a kijelzési szekcióban van kódolva; minden kiírás után törtszorzással csökkentjük spektrumvonalakat, ezzel lágy visszahullást kapunk.

Végül következzen a teljes asm kód:

.include "m88adef.inc"


; makrók

.macro delay40us
	ldi	temp3,3
	rcall	delay
.endmacro
.macro delay100us                 ; fix idejû késleltetõk
	ldi	temp3,8
	rcall	delay
.endmacro

.macro delay200us
	ldi	temp3,16
	rcall	delay
.endmacro

.macro delay10ms
	ldi	temp3,255
	rcall	delay
	rcall	delay
	rcall	delay
.endmacro

.macro lcd0                      ; LCD port 0-ba vezérlése
	cbi	lcdPort@0,lcd@0
.endmacro

.macro lcd1                      ; LCD port 1-be vezérlése
	sbi	lcdPort@0,lcd@0
.endmacro

.macro lcdCmd                    ; karakteres módba váltás
	cbi	lcdPortRS,lcdRS
.endmacro

.macro lcdChr                    ; parancsmódba váltás
	sbi	lcdPortRS,lcdRS
.endmacro

.macro write
	ldi	temp, @0
	rcall	lcdWrite
.endmacro

.macro lcdSendE
;	delay40us
	lcd1	E
	delay40us
	lcd0	E
.endmacro

; LCD vezérlés portok

.equ	lcdPort4    = portD
.equ	lcdPort5    = portD
.equ	lcdPort6    = portD
.equ	lcdPort7    = portB
.equ	lcdPortE    = portB
.equ	lcdPortRS   = portB

.equ	lcdDdr4     = ddrD
.equ	lcdDdr5     = ddrD
.equ	lcdDdr6     = ddrD
.equ	lcdDdr7     = ddrB
.equ	lcdDdrE     = ddrB
.equ	lcdDdrRS    = ddrB

.equ	lcd4        = portD5 ; 11 láb - LCD egység csatlakozásai a uC-hez
.equ	lcd5        = portD6 ; 12 láb   LCD modul D0-D3 adatlábak és az
.equ	lcd6        = portD7 ; 13 láb   R/W láb GND-hez van kötve!
.equ	lcd7        = portB0 ; 14 láb
.equ	lcdE        = portB1 ; 15 láb
.equ	lcdRS       = portB2 ; 16 láb

; konstansok

.equ	N = 32               ; idõablak elemszáma

; változók regiszterekben

.def	sin0           = R20
.def	sin1           = R2
.def	cos0           = R21
.def	cos1           = R3
.def	dispTimer      = R4
.def	temp           = R16
.def	temp2          = R17
.def	temp3          = R18
.def	temp4          = R19
.def	aksi2          = R22
.def	aksi3          = R23
.def	freqCount      = R24 ; frek mintaszámláló
.def	timeCount      = R25 ; idõ mintaszámláló


; vátozók a memóriában

.dseg
.org SRAM_START

time_buf: .byte N   ; idõtart buffer 
freq_buf: .byte N/2 ; FFT freq buffer (peak)


; kódmemória

.cseg
.org 0
	rjmp	init


; *** sampler IRQ ***

.org	ADCCaddr
	push	temp
	in	temp,SREG
	push	temp
	lds	temp,ADCH
	subi	temp,128
	st	Y+,temp
	dec	timeCount
	brne	skip_samp_1
		pop	temp
		out	SREG,temp
		pop	temp
		clt
		reti
skip_samp_1:
		pop	temp
		out	SREG,temp
		pop	temp
		reti

; *** init ***

init:

; verem
	ldi	temp, low(RAMEND)
	out	SPL, temp
	ldi	temp, high(RAMEND)
	out	SPH, temp

;rjmp main

; delay idõzítõ elindítása
	ldi	temp3,4; osztás: clk/1024 (51.2us) MOD: clk/256 (12,8us)
	out	TCCR0B,temp3

; kijelzõ idõzítõ elindítása
/*	ldi	temp3,2; clk/8 -> 16bites számláló átford.frek:38.15Hz
	sts	TCCR1B,temp3
*/
; LCD init
	ldi	temp3,255
	rcall	delay
	ldi	temp3,255
	rcall	delay

	; LCD portok kimenetre állítása
	sbi	lcdDdr7,lcd7
	sbi	lcdDdr6,lcd6
	sbi	lcdDdr5,lcd5
	sbi	lcdDdr4,lcd4
	sbi	lcdDdrE,lcdE
	sbi	lcdDdrRS,lcdRS
	lcd0	E
	lcd0	RS

	; Func reset - itt még 8 bites módban
	lcd0	7   ; function reset parancs bebitelése
	lcd0	6
	lcd1	5
	lcd1	4
	lcdSendE           ; elso reset kiküldés
	delay10ms
	lcdSendE           ; második reset kiküldés
	delay200us
	lcdSendE           ; harmadik reset kiküldés
	delay200us
	lcd0	4          ; D4=0 -> 4 bites módra váltunk
	lcdSendE           ; negyedik reset kiküldés 4 bit móddal
	delay200us

	; Init folyt 4 bites módban
	write	0b00101000 ; func set
	write	0b00001000 ; Off
	write	0b00000001 ; Clear
	delay10ms
	write	0b00000110 ; Entry Mode
	write	0b00001100 ; On
	write	0b10000000 ; pos $0


	; egyedi karakterek létrehozása
	write	$40 ; váltás egyedi karakter bitképének megadására
	lcdChr
	
	write	0b00000 ; 0
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b11111

	write	0b00000 ; 1
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b11111
	write	0b11111

	write	0b00000 ; 2
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b11111
	write	0b11111
	write	0b11111

	write	0b00000 ; 3
	write	0b00000
	write	0b00000
	write	0b00000
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111

	write	0b00000 ; 4
	write	0b00000
	write	0b00000
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111

	write	0b00000 ; 5
	write	0b00000
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111

	write	0b00000 ; 6
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111

	write	0b11111 ; 7 (ilyen amúgy van, de így egyszerûbb)
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111
	write	0b11111


; mintavételezés indítása

	sei
	ldi	temp,(1<<REFS0)|(1<<ADLAR)
	sts	ADMUX,temp
	ldi	temp,(1<<ADPS0)|(1<<ADPS1)|(1<<ADPS2)|(1<<ADIE)|(1<<ADATE)|(1<<ADEN)
	sts	ADCSRA,temp
	ori	temp,(1<<ADSC)
	sts	ADCSRA,temp

	clr	dispTimer

; *** main ***

main:

; egy ablak minta mintavételezése
	ldi	YL,low(time_buf)
	ldi	YH,high(time_buf)
	ldi	timeCount,N
	sei
	set
	brts	PC
	cli

; DFT és peak keresés

	ldi	ZL,low(data_dft*2)           ; együtthatók pointer (ez egymenetben végigmegy)
	ldi	ZH,high(data_dft*2)

	ldi	XL,low(freq_buf)             ; frekvencia pointer (abszolútértékes és csúcsértékes!)
	ldi	XH,high(freq_buf)

	ldi	freqCount,N/2                ; freq ciklusváltozó (lefelé számoló)

	loop_dft_freq:                       ; freq ciklus kezdete (DFT külsõ ciklus)
		ldi	YL,low(time_buf)     ; idõ
		ldi	YH,high(time_buf)

		ldi	timeCount,N          ; time ciklusváltozó (lefelé számoló)

		clr	sin0                 ; szinusz összetevõ összegzõregiszterek törlése
		clr	sin1
		clr	cos0                 ; koszinusz összetevõ összegzõregiszterek törlése
		clr	cos1

		loop_dft_time:               ; time ciklus kezdete (DFT belsõ ciklus)
			ld	aksi2, Y+    ; minta beolvasása a time pufferbõl
			lpm	aksi3, Z+    ; együttható (szinusz) beolvasása a programmemóriából (konstans tömb)
			fmuls	aksi2, aksi3 ; szorzás az együtthatóval
			clr	R0
			sbrc	R1,7
			com	R0
			add	sin0, R1     ; elõjeles három bájtos összegzés végrehajtása
			adc	sin1, R0
			lpm	aksi3, Z+    ; együttható (koszinusz) beolvasása
			fmuls	aksi2, aksi3 ; szorzás az együtthatóval
			clr	R0
			sbrc	R1,7
			com	R0
			add	cos0, R1     ; elõjeles három bájtos összegzés végrehajtása
			adc	cos1, R0

			dec	timeCount    ; ciklus vezérlés
			brne	loop_dft_time; time ciklus vége

		; sin^2 + cos^2 végrehajtása
		; a pitagorasz mûveletbõl a gyökvonás elmarad, enélkül is lehet komparátorszinteket definiálni

		asr	sin1 ; értékigazítás (kimenet a sin0 és cos0, túlcsordulásvédelem nem lesz, de annyira nem is okoz bajt!)
		ror	sin0
		asr	cos1
		ror	cos0

		fmuls	sin0,sin0 ; szinusz összetevő négyzetre emelése

		push	R1

		fmuls	cos0,cos0 ; koszinusz összetevő négyzetre emelése

		pop	temp
		add	R1,temp ; négyzetértékek összege

		ld	R0,X ; korábbi értékkel összehasonlítás és a spektrumvonal frissítése (csúcsdetektálás)
		cp	R0,R1
		brcc	PC+2
		mov	R0,R1 
		st	X+,R0

		dec	freqCount ; ciklus vezérlés
		breq	PC+2
		rjmp	loop_dft_freq ; freq ciklus vége


; kijelzés idõzítés
/*
	sbis	TIFR1, TOV1 ; időzítő segítségével 38Hz-re korlátozza a kijelzőre írást
	rjmp	main
	sbi	TIFR1, TOV1
*/
	inc	dispTimer
	sbrs	dispTimer,2 ; egyszerűbb kijlzéslassító
	rjmp	main
	clr	dispTimer

; kijelzés

	; felsõ kijelzõ sor

	lcdCmd
	write	$80
	lcdChr

	ldi	freqCount,16

	ldi	XL,low(freq_buf)
	ldi	XH,high(freq_buf)

disp_loop_1:
	ldi	ZL,low(data_scale*2)
	ldi	ZH,high(data_scale*2)
	ld	temp,X+
	add	ZL,temp
	clr	temp
	adc	ZH,temp
	lpm	temp,Z
	swap	temp
	andi	temp,15
	dec	temp
	brpl	PC + 2
	ldi	temp, ' '
	rcall	lcdWrite
	dec	freqCount
	brne	disp_loop_1

	; alsó kijelzõ sor

	lcdCmd
	write	$80+$40
	lcdChr

	ldi	freqCount,16

	ldi	XL,low(freq_buf)
	ldi	XH,high(freq_buf)

disp_loop_2:
	ldi	ZL,low(data_scale*2)
	ldi	ZH,high(data_scale*2)
	ld	temp,X
	add	ZL,temp
	clr	temp
	adc	ZH,temp
	lpm	temp,Z
	andi	temp,15
	dec	temp
	brpl	PC + 2
	ldi	temp, ' '
	rcall	lcdWrite
	ld	temp,X     ; visszahullás vezérlés
	ldi	temp2,235
	mul	temp,temp2
	st	X+,R1
	dec	freqCount
	brne	disp_loop_2

rjmp main


; *** szubrutinok ***

lcdWrite:
	; felso 4 bit
	lcd0	7
	lcd0	6
	lcd0	5
	lcd0	4
	sbrc	temp,7
	lcd1	7
	sbrc	temp,6
	lcd1	6
	sbrc	temp,5
	lcd1	5
	sbrc	temp,4
	lcd1	4
	rcall	writeE

	; alsó 4 bit
	lcd0	7
	lcd0	6
	lcd0	5
	lcd0	4
	sbrc	temp,3
	lcd1	7
	sbrc	temp,2
	lcd1	6
	sbrc	temp,1
	lcd1	5
	sbrc	temp,0
	lcd1	4

writeE:
;	delay100us
	lcd1	E
	delay40us
	lcd0	E
	ret

delay:
	com	temp3
	out	TCNT0, temp3
	sbi	TIFR0, TOV0
	sbis	TIFR0, TOV0
	rjmp	PC - 1
	ret


; DFT8 együtthatók

; minden sor egy harmonikus
; minden sorban szinusz és koszinusz párok vannak, összesen tehát 64 érték/sor és a 16 sor (32 elembõl ennyi spektrumvonal)
; az együtthatók elõre vannak szorozva a HAMMING ablakfüggvénnyel
; a táblázatot a teljes DFT algoritmus alatt egyszer olvassuk végig, az adatok sorrendben követik egymást, mindig az az adat jön, amivel éppen dolgozni kell

data_dft:

.db 0,25,-5,25,-11,27,-18,28,-28,28,-39,26,-52,21,-65,12,-76,0,-84,-16,-88,-36,-86,-58,-79,-79,-65,-98,-47,-113,-24,-123,0,-127,24,-123,47,-113,65,-98,79,-79,86,-58,88,-36,84,-16,76,0,65,12,52,21,39,26,28,28,18,28,11,27,5,25
.db 0,25,-10,24,-20,20,-31,12,-40,0,-44,-18,-40,-40,-25,-61,0,-76,32,-79,67,-67,96,-39,112,0,109,45,87,87,48,116,0,127,-48,116,-87,87,-109,45,-112,0,-96,-39,-67,-67,-32,-79,0,-76,25,-61,40,-40,44,-18,40,0,31,12,20,20,10,24
.db 0,25,-14,21,-27,11,-33,-6,-28,-28,-9,-47,21,-52,55,-36,76,0,71,47,36,88,-20,102,-79,79,-116,23,-113,-47,-70,-104,0,-127,70,-104,113,-47,116,23,79,79,20,102,-36,88,-71,47,-76,0,-55,-36,-21,-52,9,-47,28,-28,33,-6,27,11,14,21
.db 0,25,-18,18,-29,0,-24,-24,0,-40,33,-33,56,0,46,46,0,76,-60,60,-95,0,-73,-73,0,-112,83,-83,123,0,89,89,0,127,-89,89,-123,0,-83,-83,0,-112,73,-73,95,0,60,60,0,76,-46,46,-56,0,-33,-33,0,-40,24,-24,29,0,18,18
.db 0,25,-21,14,-27,-11,-6,-33,28,-28,47,9,21,52,-36,55,-76,0,-47,-71,36,-88,102,-20,79,79,-23,116,-113,47,-104,-70,0,-127,104,-70,113,47,23,116,-79,79,-102,-20,-36,-88,47,-71,76,0,36,55,-21,52,-47,9,-28,-28,6,-33,27,-11,21,14
.db 0,25,-24,10,-20,-20,12,-31,40,0,18,44,-40,40,-61,-25,0,-76,79,-32,67,67,-39,96,-112,0,-45,-109,87,-87,116,48,0,127,-116,48,-87,-87,45,-109,112,0,39,96,-67,67,-79,-32,0,-76,61,-25,40,40,-18,44,-40,0,-12,-31,20,-20,24,10
.db 0,25,-25,5,-11,-27,28,-18,28,28,-26,39,-52,-21,12,-65,76,0,16,84,-88,36,-58,-86,79,-79,98,65,-47,113,-123,-24,0,-127,123,-24,47,113,-98,65,-79,-79,58,-86,88,36,-16,84,-76,0,-12,-65,52,-21,26,39,-28,28,-28,-18,11,-27,25,5
.db 0,25,-26,0,0,-29,33,0,0,40,-47,0,0,-56,66,0,0,76,-86,0,0,-95,104,0,0,112,-118,0,0,-123,126,0,0,127,-126,0,0,-123,118,0,0,112,-104,0,0,-95,86,0,0,76,-66,0,0,-56,47,0,0,40,-33,0,0,-29,26,0
.db 0,25,-25,-5,11,-27,28,18,-28,28,-26,-39,52,-21,12,65,-76,0,16,-84,88,36,-58,86,-79,-79,98,-65,47,113,-123,24,0,-127,123,24,-47,113,-98,-65,79,-79,58,86,-88,36,-16,-84,76,0,-12,65,-52,-21,26,-39,28,28,-28,18,-11,-27,25,-5
.db 0,25,-24,-10,20,-20,12,31,-40,0,18,-44,40,40,-61,25,0,-76,79,32,-67,67,-39,-96,112,0,-45,109,-87,-87,116,-48,0,127,-116,-48,87,-87,45,109,-112,0,39,-96,67,67,-79,32,0,-76,61,25,-40,40,-18,-44,40,0,-12,31,-20,-20,24,-10
.db 0,25,-21,-14,27,-11,-6,33,-28,-28,47,-9,-21,52,-36,-55,76,0,-47,71,-36,-88,102,20,-79,79,-23,-116,113,47,-104,70,0,-127,104,70,-113,47,23,-116,79,79,-102,20,36,-88,47,71,-76,0,36,-55,21,52,-47,-9,28,-28,6,33,-27,-11,21,-14
.db 0,25,-18,-18,29,0,-24,24,0,-40,33,33,-56,0,46,-46,0,76,-60,-60,95,0,-73,73,0,-112,83,83,-123,0,89,-89,0,127,-89,-89,123,0,-83,83,0,-112,73,73,-95,0,60,-60,0,76,-46,-46,56,0,-33,33,0,-40,24,24,-29,0,18,-18
.db 0,25,-14,-21,27,11,-33,6,28,-28,-9,47,-21,-52,55,36,-76,0,71,-47,-36,88,-20,-102,79,79,-116,-23,113,-47,-70,104,0,-127,70,104,-113,-47,116,-23,-79,79,20,-102,36,88,-71,-47,76,0,-55,36,21,-52,9,47,-28,-28,33,6,-27,11,14,-21
.db 0,25,-10,-24,20,20,-31,-12,40,0,-44,18,40,-40,-25,61,0,-76,32,79,-67,-67,96,39,-112,0,109,-45,-87,87,48,-116,0,127,-48,-116,87,87,-109,-45,112,0,-96,39,67,-67,-32,79,0,-76,25,61,-40,-40,44,18,-40,0,31,-12,-20,20,10,-24
.db 0,25,-5,-25,11,27,-18,-28,28,28,-39,-26,52,21,-65,-12,76,0,-84,16,88,-36,-86,58,79,-79,-65,98,47,-113,-24,123,0,-127,24,123,-47,-113,65,98,-79,-79,86,58,-88,-36,84,16,-76,0,65,-12,-52,21,39,-26,-28,28,18,-28,-11,27,5,-25
.db 0,25,0,-26,0,29,0,-33,0,40,0,-47,0,56,0,-66,0,76,0,-86,0,95,0,-104,0,112,0,-118,0,123,0,-126,0,127,0,-126,0,123,0,-118,0,112,0,-104,0,95,0,-86,0,76,0,-66,0,56,0,-47,0,40,0,-33,0,29,0,-26

; oszlopkijelzés kódolt adatok
; az adott spektrumvonal egy bájtos értékével megindexelve az adatokat kapunk egy kódot, mely felső- és alsó 4 bitje kódolja a kijelezni szükséges érték felső- és alsó karakterét

data_scale:

.db 0,0,0,0,1,2,2,3,4,4,5,5,6,6,6,7,7,7,7,8,8,8,8,24,24,24,24,24,24,40,40,40,40,40,40,40,40,56,56,56,56,56,56,56,56,72,72,72,72,72,72,72,72,72,72,72,72,88,88,88,88,88,88,88,88,88,88,88,88,88,88,88,104,104,104,104,104,104,104,104,104,104,104,104,104,104,104,104,104,104,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136,136

A kód egyelőre nagyon alfa verzió, akár égbekiáltó hibát is tartalmazhat (sőt tartalmaz is, bár inkább csak hiányosság*1), de működik.

Frissítés 2017.04.13.: Frissült a kód, egy-két javítás vagy inkább csak módosítással. A jelenlegi verzióban 13284 órajel a DFT futásideje, ez 664.2μs.

Frissítés 2017.04.18.: Frissült az együtthatós adatrész, most már nem Hamming szerű, hanem tiszta Hamming ablakfüggvény dolgozik. Nagy változást amúgy nem hozott, elvben kicsit javít a spektrumvonalak szeparációjában. Az új együtthatókat érdemes lehet külön fájlba szervezni és a .include direktívával behúzni a kód végére, de akár felül is írhatod vele a kódot.
Az együtthatók fájlja: datas.asm
A létrehozó C forráskód: datas.c