example 3-11.asm

《基于TI DSP的通用算法实现》程序代码

;Example 3 - 11. Float-Point Complex FFT Implementation Based on BFP Method for TMS320C54x

;*********************************************************************
; Function: cfft64
; Version : 1.00
; Description: complex FFT
;
; Copyright Texas instruments Inc, 1998
;–––––––––––––––––––––––––––––––––
; Revision History:
;
; 0.00 M. Christ. Original code
;
; 0.01 M. Chishtie.12/96.
;– Improved radix–2 bfly code form 9 cycles to 8.
;– Combined bit–reversal in to COMBO5XX macro to save cycles.
;– Improved STAGE3 macro to 31 cycles
;
; 1.00Beta R. Piedra, 8/31/98.
;– C–callable version.
;– Removed bit–reversing and made it a separate optional function
; that also support in–place bit–reversing. In this way the FFT can
; be computed 100% in–place (memory savings)
;– Modifed STAGE3 macro to correct functional problem
;– Modified order of xmem, ymem operands in butterfly code
; to reduce number of cycles
;
; 1.00 A. Aboagye 10/15/98
; – added scale option as a parameter
;
; 1.00BFP A. Chhabra 11/09/99
;– incorporated Block Floating Point concept
;– adjustable scaling factor dictated by “cmprval_2”
;
;*********************************************************************



	N .set 64 ; NUMBER OF POINTS FOR FFT


	.include “macros.asm”
	.include “sintab.q15”
	.mmregs
		; Far–mode adjustment
		; –––––––––––––––––––
	.if __far_mode
		offset .set 1 ; far mode uses one extra location for ret addr ll
	.else
		offset .set 0
	.endif

	.asg (0), DATA
	.asg (1), SIN45
	.asg (2), save_ar7 	; stack description
	.asg (3), save_ar6 	; stack description
	.asg (4), save_ar1
	.asg (5), ret_addr
	.asg (6+offset), scale
				; x in A
;*********************************************************************
;Setting the bit growth test value.
;Depending on the bit growth quantity desired for implementation,
;include an appropriate “cmprval_2” value.
;*********************************************************************
;cmprval_2: .equ 32768*5/10 ; 0.5 decimal
cmprval_2: .equ 32768*8284/10000 ; 0.8284 decimal
;*********************************************************************

	.ref InvYeTable
	.def _cfft64

	.text
_cfft64
		; Preserve registers
		; ––––––––––––––––––
	pshm ar1
	pshm ar6
	pshm ar7
		; Preserve local variables
		; ––––––––––––––––––––––––
	frame –2
	nop

		; Get Arguments
		; –––––––––––––
	stl a,*sp(DATA) ; DATA = *SP(DATA)

   .if N>4 ; ??? no need
	st #5a82h,*sp(SIN45)
   .endif

		; Set modes
		; –––––––––
	stm #0100010101011110b,ST1 	; ASM=–2 , FRACT=1, sbx=1;CPL=1(compiler)
	stm #0, *ar5 			; initialize the contents of AR5

		; Execute
		; –––––––

********* Modifications by AC 03/25/99 **********
	mvdk *sp(DATA), ar1 		; Transfer first value of input buffer – which
					; currently
					; contain the inputs to the first stage butterflies
					; – into AR1.
					; Further manipulation of input set can be done
					; by addressing AR1
	call max_abs
	exp a 				; determine the scale up shift quantity
	stm #127, brc 			; scale for whole input array
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	rptb end_upscale–1 		; begin loop
	ld *ar1, a
	norm a 				; Scale up all input values
	sth a, *ar1+ 			; Put rescaled value back into memory;
					; increment counter to shift next value
end_upscale:
	nop
	combo5xx 			; FFT CODE for STAGES 1 and 2
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	call max_abs
	ld #cmprval_2, b 		; load threshold 0.4167 into Acc B	
	max a 				; Acc A will contain the larger of the earlier
					; MaxAbs value and the current threshold value
					; of 0.4167
	sub #cmprval_2,a,b 		; If diff > 0, thenAcc A = maxabs value
					; GOTO scaling_2
					; If diff = 0, then Acc A = threshold value
					; all values in input array are less than this...
					; GOTO performing regular next stage of bfly
   .if cmprval_2 = (32768*8284/10000)
	cc scaling_2, bgt 		; perform scaling if MaxAbs > cmprval_2
   .else

	stm #127, brc 			; scale for whole input array
	stm #0800h, ar1 		; Reset pointer to beginning of input array

	bc loop1, beq 			; execute next 1 instruction if diff>0.
	rptb loop1–1 			; begin loop
	ld *ar1, a
	sfta a, –1 			; shift down by factor 2
	stl a, *ar1+ 			; restore new value into ar1
loop1: .endif

stage3 					; MACRO WITH CODE FOR STAGE 3
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	call max_abs
	ld #cmprval_2, b 		; load threshold cmprval_2 into Acc B
	max a 				; Acc A will contain the larger of the earlier
					; MaxAbs value and the current threshold value
					; of cmprval_2
	sub #cmprval_2,a,b 		; If diff > 0, then Acc A = maxabs value
					; GOTO scaling_2
					; If diff = 0, then Acc A = threshold value
					; GOTO performing regular next stage of bfly
    .if cmprval_2 = (32768*8284/10000)
	cc scaling_2, bgt 		; perform scaling if MaxAbs > cmprval_2
   .else

	stm #127, brc 			; scale for whole input array
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	bc loop2, beq 			; execute next 1 instruction if diff>0.
	rptb loop2–1 			; begin loop
	ld *ar1, a
	sfta a, –1 			; shift down by factor 2
	stl a, *ar1+ 			; restore new value into ar1
loop2: .endif

	stdmacro 4,4,8,16,sin4,cos4 	; stage,outloopcnter,loopcnter,index
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	call max_abs
	ld #cmprval_2, b 		; load threshold cmprval_2 into Acc B
	max a 				; Acc A will contain the larger of the earlier
					; MaxAbs value and the current threshold value
					; of cmprval_2
	sub #cmprval_2,a,b 		; If diff > 0, then Acc A = maxabs value
					; GOTO scaling_2
					; If diff = 0, then Acc A = threshold value
					; GOTO performing regular next stage of bfly
   .if cmprval_2 = (32768*8284/10000)
	cc scaling_2, bgt ; perform scaling if MaxAbs > cmprval_2
   .else

	stm #127, brc 			; scale for whole input array
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	bc loop3, beq 			; execute next 1 instruction if diff>0.
	rptb loop3–1 			; begin loop
	ld *ar1, a
	sfta a, –1 			; shift down by factor 2
	stl a, *ar1+ 			; restore new value into ar1
loop3: .endif

	stdmacro 5,2,16,32,sin5,cos5 	; stage,outloopcnter,loopcnter,index
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	call max_abs
	ld #cmprval_2, b 		; load threshold cmprval_2 into Acc B
	max a 				; Acc A will contain the larger of the earlier
					; MaxAbs value and the current threshold value
					; of cmprval_2
	sub #cmprval_2,a,b 		; If diff > 0, then Acc A = maxabs value
					; GOTO scaling_2
					; If diff = 0, then Acc A = threshold value
					; GOTO performing regular next stage of bfly
   .if cmprval_2 = (32768*8284/10000)
	cc scaling_2, bgt ; perform scaling if MaxAbs > cmprval_2
   .else

	stm #127, brc 			; scale for whole input array
	stm #0800h, ar1 		; Reset pointer to beginning of input array
	bc loop4, beq 			; execute next 1 instruction if diff>0.
	rptb loop4–1 			; begin loop
	ld *ar1, a
	sfta a, –1 			; shift down by factor 2
	stl a, *ar1+ 			; restore new value into ar1
loop4: .endif

	laststag 6,sin6,cos6 		; MACRO WITH CODE FOR STAGE 7
	bd end_lab
	nop
	nop


*********************************************************************
;MAX_ABS
;=======
;
;Perform comparison of consecutive values in order to obtain maximum
;absolute value in the array of inputs. Steps to do this:
;
;(i) Place consecutive values in acc A and B respectively
;(ii) Compute their absolute values
;(iii) Find the MAX of these two accumulators values
;(iv) Monitor the Carry bit and determine which acc contains MAX value.
;(v) Store the max value in acc A
;(vi) Take in the next value in the input array and load into acc B.
; Compute its absolute.
;(vii) Go back to step (ii)
;
;*** Steps (iv) and (v) above are performed in combination
; as a result of the C54x “MAX” instruction.
**********************************************************************

max_abs:
					; Set breakpoint to verify that AR1 does point to correct address
	ld *ar1+, a
	ld *ar1+, b
	absa 				; setup absolute value for max comparison
	absb
	stm #126, brc 			; 126 values remain to be read of the 128 value
					; input array. The loop executes 127 times.
	rptb find_max–1
	max a
	ld *ar1+, b 			; enter the next value in the input array into acc B
	abs b 				; setup next value for absolute max comparison
find_max
	ret 				; returns the maximum absolute value in the
; Acc A
************************ END of “MAX_ABS” routine ********************



**********************************************************************
;
;SCALING
;=======
;
;This routine performs the following in order:
;
;(i) Now, since reaching this routine implies that bit growth is
; likely at the output of this stage, scale the mantissa values.
; Scaling factor is determined by the reciprocal of the expected
; bit growth.
; Expected bit growth = 2.4
; Reciprocal = 0.4167
; If maximum of input block to a stage is > 0.4167, then scale
; down once by shifting right.
; If maximum of input block is twice > 0.4167, then scale down
; by shifting right twice.
; Else ignore scaling and proceed as normal.
;
;Arriving at this routine implies that Acc A already contains
;the maximum absolute value of the input array.
;
**********************************************************************

scaling_2: 				; Acc A contains the MaxAbs value
	pshm ar1
	pshm ar2
	ssbx SXM
	rsbx FRCT
					;rsbx TC ; this is ID for div–by–1.207 selection
					; in the reciprocal routine
	call reciprocal 		; Acc A contains MaxAbs
					; COMPUTE SCALING FACTOR = (1/(MaxAbs*2.41421356))
					; now, remember that the original DSPLIB FFT code contains
					; automatic scaling down within each radix–2 stage
					; by a factor of 2, i.e. one bit place.
					; As a result, what we manually need to tweak with for
					; scaling is a value = (2.41421356/2) = 1.207
					; => Scaling factor needed here = (1/(MaxAbs * 1.207))
	ld a, b 			; temporarily store Acc A into B
	add *ar5, a 			; add scaling factor to Acc A
	stl a, *ar5 			; save the scaling factor into AR5
	stm #0800h, ar1 		; re–align to the beginning of the array
	stl a, *ar2
	ld *ar2, t
	stm #127, brc
	rptb end_scale2–1
	mpy *ar1, a
	sth a, *ar1+
end_scale2:

	popm ar2
	popm ar1
	ret
*************************End of “SCALING” routine ********************


**********************************************************************
;
;RECIPROCAL
;==========
;
;This routine will perform the following:
;
; 1) Take the value in Acc A and compute its reciprocal
;2) The result is available in two portions, r and rexp.
;3) MPY r and rexp. however this may lead to a value greater
; than 32768.
;4) Now MPY (3) by reciprocal of 4, i.e. by 0.25. This is also
; equivalent to right shifting by 2 bit locations.
;5) Since (3) will likely lead to a value greater than 32768, it
; is probably better to perform (4) on value “r” from (2).
; This way the value is decreased by a factor of 4.
; Now, we can mutiply by rexp without exceeding 32768.
;
;NOTE: On entering this routine, the Acc A contains the MaxAbs value
; on which reciprocal computation has to be performed.
;
**********************************************************************

reciprocal:
		;–––––––––––––––––––––––––––––––––––––––––––––––––––––
		; Set offsets to local function variables defined on stack
		;––––––––––––––––––––––––––––––––––––––––––––––––––––
	.asg 0, SP_INVYETABLE
	.asg 1, SP_XNORM
	.asg 2, SP_TEMP
	.asg 3, FRAME_SZ
;––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
; Assign registers to local variables
;––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
	.asg ar0, AR_X
	.asg ar1, AR_Z
		;.asg brc, AR_N
	.asg ar3, AR_ZEXP
	.asg ar4, AR_TABLE
;––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
; Process command–line arguments
;–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
		; stl a,*(AR_X) ; Acc A contains MaxAbs value
	pshm ar0
	pshm ar1
	pshm ar3
	pshm ar4
	pshm ar5
		;––––––––––––––––––––––––––––––––––––––––––––––––––––
		; Initialize constants
		;––––––––––––––––––––––––––––––––––––––––––––––––––––
	st #InvYeTable,*sp(SP_INVYETABLE)