*********************************************************************************
*       (C) COPYRIGHT TEXAS INSTRUMENTS, INC. 1996                              *
*********************************************************************************
*                                                                               *
* MODULE NAME:  unpack.asm                                                      *
*                                                                               *
* AUTHORS:      Simon Lau and Nathan Baltz                                      *
*                                                                               *
* DESCRIPTION:  PHASE THREE & FOUR      Unpacking to 2N Outputs                 *
*               This function is called from the main module of the 'C54x Real  *
*               FFT code.  It first computes four intermediate sequences (RP,   *
*               RM, IP, IM) from the resulting complex sequence at the end of   *
*               the previous phase.  Next, it uses the four intermediate        *
*               sequences to form the FFT of the original 2N-point real input.  *
*               Again, the outputs are divided by 2 to prevent overflow.        *
*                                                                               *
* REGISTER USAGE:   AR0 index of twiddle tables                                 *
* (PHASE THREE)     AR2 pointer to R[k], I[k], RP[k], IP[k]                     *
*                   AR3 pointer to R[N-k], I[N-k], RP[N-k], IP[N-k]             *
*                   AR6 pointer to RM[k], IM[k]                                 *
*                   AR7 pointer to RM[N-k], IM[N-k]                             *
*                                                                               *
* REGISTER USAGE:   AR0 index of twiddle tables                                 *
* (PHASE FOUR)      AR2 pointer to RP[k], IP[k], AR[k], AI[k], AR[0]            *
*                   AR3 pointer to RM[k], IM[k], AR[2N-k], AI[2N-k]             *
*                   AR4 pointer to cos(k*pi/N), AI[0]                           *
*                   AR5 pointer to sin(k*pi/N), AR[N], AI[N]                    *
*                                                                               *
* DATE:         7-16-1996                                                       *
*                                                                               *
*********************************************************************************

        .mmregs
        .include    "fft_size.inc"
        .def        unpack
        .ref        fft_data, sine, cosine


        .text
                                                                                                                   
unpack:

; Compute intermediate values RP, RM, IP, IM

        .asg    AR2,XP_k
        .asg    AR3,XP_Nminusk
        .asg    AR6,XM_k
        .asg    AR7,XM_Nminusk

        STM     #fft_data+2,XP_k                    ; AR2 -> R[k] (temp RP[k])
        STM     #fft_data+2*K_FFT_SIZE-2,XP_Nminusk ; AR3 -> R[N-k] (temp RP[N-k])
        STM     #fft_data+2*K_FFT_SIZE+3,XM_Nminusk ; AR7 -> temp RM[N-k]
        STM     #fft_data+4*K_FFT_SIZE-1,XM_k       ; AR6 -> temp RM[k]
        STM     #-2+K_FFT_SIZE/2,BRC
        RPTBD   phase3end-1
        STM     #3,AR0

        ADD     *XP_k,*XP_Nminusk,A             ; A := R[k]+R[N-k] = 2*RP[k]
        SUB     *XP_k,*XP_Nminusk,B             ; B := R[k]-R[N-k] = 2*RM[k] 
        STH     A,ASM,*XP_k+                    ; store RP[k] at AR[k]
        STH     A,ASM,*XP_Nminusk+              ; store RP[N-k]=RP[k] at AR[N-k]          
        STH     B,ASM,*XM_k-                    ; store RM[k] at AI[2N-k]
        NEG     B                               ; B := R[N-k]-R[k] = 2*RM[N-k]
        STH     B,ASM,*XM_Nminusk-              ; store RM[N-k] at AI[N+k]

        ADD     *XP_k,*XP_Nminusk,A             ; A := I[k]+I[N-k] = 2*IP[k]
        SUB     *XP_k,*XP_Nminusk,B             ; B := I[k]-I[N-k] = 2*IM[k]
        STH     A,ASM,*XP_k+                    ; store IP[k] at AI[k]
        STH     A,ASM,*XP_Nminusk-0             ; store IP[N-k]=IP[k] at AI[N-k]
        STH     B,ASM,*XM_k-                    ; store IM[k] at AR[2N-k]
        NEG     B                               ; B := I[N-k]-I[k] = 2*IM[N-k]
        STH     B,ASM,*XM_Nminusk+0             ; store IM[N-k] at AR[N+k]

phase3end:
        ST      #0,*XM_k-                       ; RM[N/2]=0 
        ST      #0,*XM_k                        ; IM[N/2]=0


; Compute AR[0], AI[0], AR[N], AI[N]

        .asg    AR2,AX_k
        .asg    AR4,IP_0
        .asg    AR5,AX_N

        STM     #fft_data,AX_k                  ; AR2 -> AR[0] (temp RP[0])
        STM     #fft_data+1,IP_0                ; AR4 -> AI[0] (temp IP[0])
        STM     #fft_data+2*K_FFT_SIZE+1,AX_N   ; AR5 -> AI[N]
        ADD     *AX_k,*IP_0,A                   ; A := RP[0]+IP[0]
        SUB     *AX_k,*IP_0,B                   ; B := RP[0]-IP[0]
        STH     A,ASM,*AX_k+                    ; AR[0] = (RP[0]+IP[0])/2
        ST      #0,*AX_k                        ; AI[0] = 0
        MVDD    *AX_k+,*AX_N-                   ; AI[N] = 0
        STH     B,ASM,*AX_N                     ; AR[N] = (RP[0]-IP[0])/2


; Compute final output values AR[k], AI[k]

        .asg    AR3,AX_2Nminusk
        .asg    AR4,COS
        .asg    AR5,SIN

        STM     #fft_data+4*K_FFT_SIZE-1,AX_2Nminusk    ; AR3 -> AI[2N-1] (temp RM[1])
        STM     #cosine+512/K_FFT_SIZE,COS              ; AR4 -> cos(k*pi/N)
        STM     #sine+512/K_FFT_SIZE,SIN                ; AR5 -> sin(k*pi/N)
        STM     #K_FFT_SIZE-2,BRC
        RPTBD   phase4end-1
        STM     #512/K_FFT_SIZE,AR0             ; index of twiddle tables

        LD      *AX_k+,16,A                     ; A := RP[k] || AR2->IP[k]
        MACR    *COS,*AX_k,A                    ; A := A+cos(k*pi/N)*IP[k]
        MASR    *SIN,*AX_2Nminusk-,A            ; A := A-sin(k*pi/N)*RM[k] 
                                                ; || AR3->IM[k]   
        LD      *AX_2Nminusk+,16,B              ; B := IM[k] || AR3->RM[k]
        MASR    *SIN+0%,*AX_k-,B                ; B := B-sin(k*pi/N)*IP[k]
                                                ; || AR2->RP[k]
        MASR    *COS+0%,*AX_2Nminusk,B          ; B := B-cos(k*pi/N)*RM[k]
        STH     A,ASM,*AX_k+                    ; AR[k] = A/2
        STH     B,ASM,*AX_k+                    ; AI[k] = B/2
        NEG     B                               ; B := -B
        STH     B,ASM,*AX_2Nminusk-             ; AI[2N-k] = -AI[k] = B/2
        STH     A,ASM,*AX_2Nminusk-             ; AR[2N-k] = AR[k] = A/2 

phase4end:
        RET                                     ; return to Real FFT main module

        .end