Principe for Discrete Fouriertransform using cordic

Reinhard

I am an absolute P2 beginner and learn by myself with concrete examples the p2asm language.
some of them I like to share.

{{
	simple demo for a 8point DFT 
	show only the princip to do it with cordic
	remember DFT
	
           N-1
	X[k] = sum (x[n] * exp(-j * 360° * n*k/N))  0 <= k <= N-1
	       n=0	
	       
	testcase: x[n(0...7)]= 100  0   0   0   0   0   0   0
	expected: X[k(0...7)]= 100 100 100 100 100 100 100 100 
	proof: ok
	metric: test with spinsim and analyse logfile       
}}

dat
	orgh	200
	samples		res 8*2			
	
	org 0
	mov		ptra,##$200			'this is a testsignal (impuls)
	wrword	#100,ptra++			'8 samples a 16 bit
	wrword	#0,ptra++			'normaly this input is from an adc
	wrword	#0,ptra++			'impuls, because this has an 
	wrword	#0,ptra++			'known dft response (FT)
	wrword	#0,ptra++			'make it easier to test w/o hardware	
	wrword	#0,ptra++			'if pulse	| . . . . . . .
	wrword	#0,ptra++			'the FT is  | | | | | | | | 
	wrword	#0,ptra++			'the DC component is not calculated
								'in this example
								'because it is the simple sum
								'of the input samples
	
'---------------------------------

	mov	s,#0					'at begin the shift operator is zero
	mov	k,#7					'calculate 7 frequency components
								'btw. 3 was enough for real app
								'4,5,6,7 are mirror or negative freq.
	nop							'the DC component is not calculated
	
LP	mov	ptra,##$200				'point to first sample
	mov	arg,#0					'the twiddle factor starts at 0 degree
	mov	arg2,#0					'needed for temporary calculations		
	mov	n,#8					'8 samples
	mov	sum_real,#0
	mov	sum_imag,#0

l1	rdword	xn,ptra++			'get next sample
	qrotate	xn,arg				'calculate x[n] * exp(-j * arg)		
	getqx	x					'get real part
	getqy	y					'get imaginary part
	add		sum_real,x			'sum real
	add		sum_imag,y			'sum imaginary	
	add		arg2,##$2000_0000	'increment the phase 45°  360/8
	mov		arg,arg2			'(8 point dft)
	shl		arg,s				'adjust the speed of rotation
	djnz	n,#l1				'all samples done ?
'---------------------------------
								'(sum_real,sum_imag) can here go for
								'further processing
								'i.e. out to dac or serial or ... 
								'... build abs/arg with qvector
	add		s,#1				'increment the speed of rotation
	djnz	k,#LP				'all frequency components done?

	nop
	nop
'---------------------------------

	sum_real	res	1
	sum_imag	res	1
	arg		res	1
	arg2	res	1
	xn		res	1
	x		res	1
	y		res	1
	n		res	1
	s		res	1
	k		res	1

cgracey

Yes, CORDIC lets you do two things at once: generate a pre-scaled sine and cosine, all ready for addition into the accumulators. No need to post-scale sine and cosine values.

Eventually, that can be recoded to feed the CORDIC every 8 clocks with scaled sine/cosine commands.

Reinhard

cgracey wrote: »

Eventually, that can be recoded to feed the CORDIC every 8 clocks with scaled sine/cosine commands.

Yes, next I 'll try to interleave the cordic feed for performance.
But first I fixed a bug in my code and expand it to 256point DFT.

dat

ORG        0           					'cog code
'///////////////////////////////////////////////////////////////////////
'    			FILL SAMPLE BUFFER WITH TEST DATA
'///////////////////////////////////////////////////////////////////////
	mov			ptra,##$02000			'fill sample buffer
	wrword		#100,ptra++
	rep			#1,#256					
	'wrword		#0,ptra++				'puls  ok
	wrword		#100,ptra++				'step  ok
'///////////////////////////////////////////////////////////////////////
'    			SIMPLE CALCULATION OF DC COMPONENT
'///////////////////////////////////////////////////////////////////////
DFT_DC	
	mov			ptra,##$02000
	rep			#2,#256
	rdword		xn,ptra++
	add			sum_real,xn
	wrlong		sum_real,##$3000
'///////////////////////////////////////////////////////////////////////
'    			DISCRET FOURIER TRANSFORM
'///////////////////////////////////////////////////////////////////////
DFT_AC	
	mov		k,#255
	mov		s,#1
L1	
	mov		n,#256
	mov		arg,#0
	mov		arg2,#0
	mov		sum_real,#0
	mov		sum_imag,#0
	mov		ptra,##$02000
	
L0	

	rdword		xn,ptra++
	qrotate		xn,arg			'calculate x[n] * exp(-j * arg)	
	getqx	x					'get real part
	getqy	y					'get imaginary part
	add		arg2,dphi
	mov		arg,arg2	
	mul 	arg,s
	shl		arg,#16
	sumc	sum_real,x			'sum real
	sumc	sum_imag,y			'sum imaginary	
	djnz 	n,#L0
	
	wrlong	sum_real,##$3000	'write result to somewhere
	wrlong	sum_imag,##$3004
	
	add		s,#1
	djnz 	k,#L1
		
'///////////////////////////////////////////////////////////////////////
'    		  END DFT
'///////////////////////////////////////////////////////////////////////
	jmp	$



dphi	long	$100		'($100 << 16) * 256 = $1_0000_0000 = 360°
k		res		1
arg		res		1
arg2	res		1 
s		res		1
n		res		1
sum_real	res	1	
sum_imag	res	1
xn		res	1
x		res	1
y		res	1

ORGH       $02000      						'hub code
samples

Reinhard

the bug was

	mov		arg,arg2			'(8 point dft)
	shl		arg,s		

right is

         mov		arg,arg2	
	mul 	               arg,s
	shl		     arg,#16