$mod51
$nolist
$nosymbols
$include(c:\cet\inc\controls.inc)
$list
; monitor routines:
funbase	equ	100h		; monitor jump table
chkbrk	equ	funbase+42	; keyboard break routine
sndchr	equ	funbase+12	; output character to screen
ilprt	equ	funbase+18	; print string pointed by stack to screen
prtrad	equ	funbase+21	; print hex of value in A 
dvide8	equ	funbase+39	; 16 by 8 bit division
AD0	equ	0e000h		; address of A/D chip

codeseg		equ	3000h
tableseg	equ	codeseg + (8*256)
arrayseg	equ	tableseg + (4*256)
rawdat		equ	4000h
rate		equ	8000		; sample rate
srate		equ	not 125	; 1/8000
offset		equ	48	; data input offset

temp		equ	r2		; temp storage
sine		equ	r3		; current sine value
cosine		equ	r4		; current cosine value
treal		equ	r5		; temp real comp
timag		equ	r6		; temp imag comp
lctr1		equ	r7		; loop counter 1

flags		data	21h
pwramp		bit	21h.0
normalized	bit	21h.1
ready		bit	21h.2

celct		data	24h		; cell counter
celnum		data	25h		; number of cells
pairnm		data	26h		; pairs per cell
celdis		data	27h		; span between pairs
delta		data	28h		; angle increment
sclfct		data	29h		; multiply output amplitude by 2^sclfct


counter		data	45h		; plotter loop
rdh		data	46h		; raw data hi addr
rdl		data	47h		; raw data lo addr

; since all arrays are on page boundaries
; only the low byte is used for indexing
head	data	4Ah		; plot pointer
rlpt1	data	4Bh		; real ptr 1
rlpt2	data	4Ch		; real ptr 2
impt1	data	4Dh		; imag ptr 1
impt2	data	4Eh		; imag ptr 2
sinpt	data	4Fh		; sine table ptr


	cseg at codeseg
start:
	clr	normalized	; 0 = perform 2's comp of data move
	setb	pwramp		; do amplitude, 0=do power
	mov	psw,#0		; reg base 0
	mov	dptr,#rawdat
	mov	rdl,dpl		; raw data pointer
	mov	rdh,dph
;	call	wave		; create data
	call	init_sampler	; setup data collector
	jmp	main
;
window:	mov	a,rdl
	add	a,#16
	mov	rdl,a
	mov	a,rdh		; raw data high byte address
	addc	a,#0
	cjne	a,#80h,notbot
	mov	a,#41h
notbot:	mov	rdh,a
	ret
;
;$include(wave.asm)		; create test data
;	OR
$include(sampler.asm)		; collect test data
;
$include(power.asm)		; power spectrum
;
$include(scale.asm)		; keep data from overflowing
;
$include(sysmpy.asm)		; signed multiply
;		
; Use one of the following main loops
;============================================================================
; this will collect data from A/D, perform FFT, display and repeat
;============================================================================
main:	call	sampler		; get data from A/D
	call	fft
	call	power		; spectrum
	call	ppas		; plot power/amplitude spectrum
	call	chkbrk		; quitter ? control-c to exit
	jmp	main
;============================================================================
; This will move existing data from raw data buff, perform fft, display,
; move to next blob of raw data and repeat 
;============================================================================
;main:	call	fft
;	call	power		; spectrum
;	call	ppas		; plot power/amplitude spectrum
;	call	chkbrk		; quitter ? control-c to exit
;	call	window		; next blob
;	jmp	main
;============================================================================
;============================================================================
	
fft:	mov	sclfct,#0		; no scale factor yet

; clear imag array
clrimg:	mov	dptr,#imag
	mov	r0,#0
	clr	a
clri10:	movx	@dptr,a
	inc	dptr
	djnz	r0,clri10

; move input data into real array
	mov	dptr,#real
	mov	p2,rdh
	mov	r0,rdl
	mov	temp,#0
movdi:	movx	a,@r0			; get input data
	clr	c
	jb	normalized,no2cmp
	subb	a,#80h			; convert to 2's compliment
no2cmp:	movx	@dptr,a			; put to real
	inc	dptr
	inc	r0
	djnz	temp,movdi
;
$include(perm2.asm)		; decimation in time
;
; FFT first pass
;
Pass1:	call	scale
	mov	dptr,#real
Pa1:	movx	a,@dptr
	mov	r1,a		; get real[i]
	inc	dpl
	movx	a,@dptr
	mov	r0,a		; get real[i+1]
	dec	dpl
	add	a,r1		; real[i+1] + real[i]
	movx	@dptr,a
	mov	a,r1		; get real[i]
	clr	c
	subb	a,r0		; real[i] - real[i+1]
	inc	dpl
	movx	@dptr,a
	inc	dpl
	mov	a,dpl
	jnz	Pa1
;
; computation of fft pass 2 thru n.
;
Fpass:	mov	a,#64			; 90 degrees
	mov	celnum,a
	mov	delta,a
	mov	a,#2
	mov	pairnm,a
	mov	celdis,a
;
Npass:	call	scale
	mov	celct,celnum

; initialize data pointers
	clr	a
	mov	rlpt1,a
	mov	rlpt2,a
	mov	impt1,a
	mov	impt2,a
;
Ncell:	mov	sinpt,#0
	mov	lctr1,pairnm		; loop counter 1
;
Nc1:	mov	a,rlpt1
	add	a,celdis		; 2,4,8,16,32,64
	mov	rlpt2,a
	mov	impt2,a
;
	mov	dph,#high stadr
	mov	dpl,sinpt
	movx	a,@dptr
	mov	sine,a
	mov	a,dpl			; bugfix - add current offset
	add	a,#64			; 90 degrees
	mov	dpl,a
	movx	a,@dptr
	mov	cosine,a
;
	mov	dph,#high real
	mov	dpl,rlpt2
	movx	a,@dptr
	mov	b,a		
	mov	temp,a
	mov	a,cosine
	call	mpy
	mov	treal,a			; tr = rn * cos(z)
	mov	b,temp
	mov	a,sine
	call	mpy
	mov	timag,a			; ti = rn * sin(z)
;
	inc	dph			; high imag
	mov	dpl,impt2
	movx	a,@dptr	
	mov	b,a
	mov	temp,a
	mov	a,sine
	call	mpy			; a = in * sin(z)
	add	a,treal
	mov	treal,a			; tr = rn*cos + in*sin
	mov	b,temp
	mov	a,cosine
	call	mpy
	clr	c
	subb	a,timag
	mov	timag,a			; ti = in*cos - rn*sin	
;
	dec	dph			; high real
	mov	dpl,rlpt1
	movx	a,@dptr
	mov	temp,a
	add	a,treal
	movx	@dptr,a			; rl1 = rl0 + treal
	mov	a,temp
	clr	c
	subb	a,treal
	mov	dpl,rlpt2
	movx	@dptr,a			; rl2 = rl0 - treal
;
	inc	dph			; high imag
	mov	dpl,impt1
	movx	a,@dptr
	mov	temp,a
	add	a,timag
	movx	@dptr,a			; im1 = im0 + timag
	mov	a,temp
	clr	c
	subb	a,timag
	mov	dpl,impt2
	movx	@dptr,a			; im2 = im0 - timag
;
	mov	a,sinpt
	add	a,delta
	mov	sinpt,a
	inc	rlpt1
	inc	impt1
	djnz	lctr1,Nc1a
	sjmp	Nc2
Nc1a:	jmp	Nc1	
;
Nc2:	mov	a,rlpt1
	add	a,celdis
	mov	rlpt1,a
	mov	impt1,a
	djnz	celct,Np1a
	sjmp	Np1
Np1a:	jmp	Ncell
;
; change parameters for next pass
;
Np1:	mov	a,celnum		; 64,32,16,8,4,2
	clr	c
	rrc	a			; halve no. cells
	mov	celnum,a
	jz	fftdun			; done out of cells
	mov	a,pairnm
	clr	c
	rlc	a
	mov	pairnm,a		; twice as many pairs
	mov	a,celdis
	clr	c
	rlc	a
	mov	celdis,a		; twice as far apart
	mov	a,delta
	clr	c
	rrc	a
	mov	delta,a			; halve the angle
	jmp	Npass	

fftdun:	ret
;
$include(plotpas.asm)

	cseg at tableseg
$include(sqrt.dat)		; used in power spectrum
$include(sint.dat)		; for FFT computation
$include(sqrut.dat)		; used for amplitude spectrum
$include(pert.dat)		; permute lookup 

sqadr	equ	sqrt		; square root table address
stadr	equ	sint		; sine table address

	dseg at arrayseg
; these two arrays MUST stay together else oops.
real:	ds	256		; result = cycles per unit
imag:	ds	256
	ds	127		; filler so loop can inc to 0
pabuf:	ds	129		; power/amplitude spectrum


	end