unpack.asm

seed傅立叶变换程序

;***********************************************************
; 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

unpack  .macro  sin,cos

;        .global    UNPACK

UNPACK:
    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, -1, *ar_x1- 	             ; y(0)Im = al	(nyquist)
    stl	   b, -1, *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 = x1Re - x2Re
    neg    b                            ; im = -(x1Re - x2Re)
    sth    a, -1, *ar_x1+               ; x1Re = rp = A
    add    *ar_x1, *ar_x2-, a           ; 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
    mac    *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
    mas    *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