unpacki.asm

CCS3.3自带的TI 5400系列DSP的dsplib文件。文档说明可以在TI公司网站上下载。

;***********************************************************
; Version 2.20.01                                           
;***********************************************************
;****************************************************************************
; Function name :   unpacki
; Description   :   unpacks an NREAL-point RIFFT input data before
;                   an N/2-point CIFFT
;
; Copyright Texas instruments Inc, 1998
;----------------------------------------------------------------------------
; History	:
;
; 1.00, A. Aboagye, 10/31/98. inverse FFT changes from unpack.asm
; 1.10, A. Aboagye, 08/23/99. fixed bugs in code. changed from macro to function
;                             use new twiddles for computation
; 1.20, A. Aboagye, 06/30/00. negated sign of im(Y(N/2))
;--------------------------------------------------------------------------

                                    ; AR0: not used
	.asg	ar1,ar_brcsave          ; preserve brc
	.asg	ar2,ar_x1               ; pointer to k
	.asg	ar3,ar_x2               ; pointer to (N-k)
	.asg	ar4,ar_x11              ; pointer to &x[1]Re
	.asg	ar5,ar_sin              ; pointer to w
	.asg	ar6,ar_x2save           ; preserve x2
	.asg	ar2,ar_cos

invpack .macro  sin,cos

;        .global    UNPACKI

UNPACKI:
    mvdk   *sp(DATA),*(ar_x1)		; pointer to x[0] Re

    stm    #(N/2 -2), brc		; brc = (n/2) -2
					; Initialize ar_x2
    stm    #2,ar0			; ar0 = 2

    ldm    ar_x1, a			; a = &x(0)
    add    #(N*2-2),a			; a = ar_x1 + (2*n-2)
    stlm   a, ar_x2			; ar_x2 = ar_x1 + (2*n-2)
					; point to x(N-1) Re

    stlm   a, ar_x2save			; preserve ar_x2

					; repeats (n/2 -1 times)
  
; ------------------------------------------------------
; Step 1: precompute y(0)Re (DC) and y(N) Re (Nyquist)
; y(0) Re		= x(0)Re + x(0)Im  = store in y(0) Re
; y(N) Re (Nyquist)	= x(0)Re - x(0)Im  = store in y(0) Im
; -------------------------------------------------------

    ld     *ar_x1+, a		        ; al = x(0)Re
    add    *ar_x1, a, b		        ; bl = x(0)Re + x(0)Im
    sub    *ar_x1, a
    stl    a, -2, *ar_x1- 	        ; y(0)Im = al	(nyquist)
    stl	   b, -2, *ar_x1+0	        ; y(0)Re = bl	(DC)
					                ; ar_x1 points to x(1)


; -----------------------------------------------
; Step 2: Compute rp, ip, rm, im  for 1 <= i <N/2
; rp = 0.5 * (x1Re + x2Re)
; ip = 0.5 * (x1Im - x2Im)
; rm = -0.5 * (x1Im + x2Im)
; im = 0.5 * (x1Re - x2Re)
; Final order in memory: rp-ip .... rm-im
; -----------------------------------------------
    stm    #3, ar0                  ; ar0 = 3
    rptbd  end_factor-1
    mvmm   ar_x1, ar_x11		    ; preserve &x(1)
    nop                             ; real

    add    *ar_x1, *ar_x2, a	    ; A = rp = x1Re + x2Re 
    sub    *ar_x1 , *ar_x2+, b	    ; B = im = x1Re - x2Re
    sth    a, -1, *ar_x1+ 	        ; x1Re = rp = A
    add    *ar_x1, *ar_x2-, a	    ; A = x1Im + x2Im
    neg    a                        ; rm = -(x1Im + x2Im)
    sth    a, -1, *ar_x2+ 	        ; x2Re = rm = A
    sub    *ar_x1 , *ar_x2 , a	    ; A = ip = x1Im - x2Im
    sth    a, -1, *ar_x1+ 	        ; x1Im = ip = A; point to next x1Re
    sth	   b, -1, *ar_x2-0	        ; x2Im = im = B; point to next x2Re
end_factor

; ---------------------------------------------
; Step 3: pre-compute y(N/2)
; rm = 0 im = 0 cos(pi/2) = 0 sin(-pi/2) = -1   aka
; rp = x(N/2)Re     ip = x(N/2)Im
; y(N/2)Re = rp +cos*ip - sin* rm = rp = x(N/2)Re
; y(N/2)Im = im -cos*rm - sin*ip =  -ip = x(N/2)Im
; ---------------------------------------------
					                ; ar_x1 = ar_x2 = &x(N/2)Re
    ld     *ar_x1, a
    stl	   a, -1, *ar_x1+ 	        ; scale x(N/2)Re and point to Im
    ld     *ar_x1, a
    neg    a                            ; negate x(N/2)Im
    stl	   a, -1, *ar_x1-


; Re-initialize registers
; --------------------------
    ; ar_x11 preserves &x[1]Re
    stm    #(N/2-2), brc                 ; brc = N/2 - 2
    stm    #cos, ar_cos                  ; start of cosine in stage 'stg'
    stm    #sin, ar_sin			 ; start of sine in stage 'stg'

; -----------------------------------------------

; Step 4:  compute y final results
; y1Re = rp + cos * rm - sin * im
; y1Im = ip + cos * im + sin * rm
; y2Re = rp - (cos * rm - sin * im)	; symmetries with y1Re
; y2Im = -ip + cos * im + sin * rm

; -----------------------------------------------

    rptbd  end_loop-1
    mvmm   ar_x2save, ar_x2		    ; ar_x2 = &x[N/2]Re
    nop
    
    mpy	   *ar_x2+, *ar_cos, b          ; b = rm * cos
    mas    *ar_x2, *ar_sin, b		; b = rm * cos - im * sin
    add    *ar_x11, 16, b, a		; a = b + rp = y1Re
    neg    b
    add	   *ar_x11+, 16, b, b		; b = b + rp = y2Re
    frame  #-2 
    nop
    sth    a, -1, *sp(0) 		    ; store y1Re on stack
    sth    b, -1, *sp(1) 		    ; store y2Re on stack
    mpy	   *ar_x2-, *ar_cos+,a	    ; a = im * cos
    mac    *ar_x2+, *ar_sin+, a    ; a = im * cos + rm * sin
    add    *ar_x11, 16, a, b        ; b =  ip + im * cos - rm * sin = y1Im
    sub    *ar_x11, 16, a           ; a =  -ip + im * cos - rm * sin = y2Im
    sth    b, -1, *ar_x11-          ; ip = y1Im
    sth    a, -1, *ar_x2-           ; im = y2Im
    popd   *ar_x11+                 ; store y1Re in rp position
    popd   *ar_x2                   ; store y2Re in rm position
    addm   #-2, *(ar_x2)
    nop
    mar    *ar_x11+
end_loop:

    .endm

;end of file. please do not remove. it is left here to ensure that no lines of code are removed by any editor