/*_______________________________________________________________________
DFT.ASM         ADSP-2106x Discrete Fourier Transform
This routine performs an N point real DFT according to the following equation:

		    N-1
real(k)+j*imag(k) = SUM input(n)[C - j*S]; k=0 to N-1
		    n=0

where: C=cos(2*pi*k*n/N), S=sin(2*pi*k*n/N), j=sqrt(-1)
_________________________________________________________________________*/

#include "def21060.h"           /* Memory Mapped IOP register definitions */
#define N 64                    /* Constant for number of points in input */

.SECTION/DM     dm_data;        /* Declare variables in data memory */
.VAR input[N]= "test64.dat";
.VAR real[N];
.VAR imag[N];


.SECTION/PM     pm_data;        /* Declare variables in program memory */
.VAR sine[N]= "sin64.dat";      /* Cosine is derived using a shifted */ 
                        /* pointer to this circular buffer.*/

.SECTION/PM     pm_rsti;           /* The reset vector resides in this space */
		NOP;
		USTAT2= 0x108421;  /* 1st instr. to be executed after reset */
		DM(WAIT)=USTAT2;   /* Set external memory waitstates to 0 */
		JUMP start;
	     


.SECTION/PM     pm_code;        /* Example setup for DFT routine */
start:          M1=1;
		M9=1;
		B0=input;
		L0=@input;              /* Input buffer is circular */
		I1=imag;
		L1=0;
		CALL dft (DB);          /* Example delayed call instruction */
		  I2=real;              /* In delay field of call */
		  L2=0;                 /*           ''            */
end:            IDLE;

/*___________________________DFT Subroutine___________________________*/
dft:            B8=sine;                /* Sine pointer */
		L8=@sine;
		B9=sine;                /* Derive cosine from sine by */
		I9=sine+N/4;            /* shifting pointer over 2pi/4 */
		L9=@sine;               /* and using a circular buffer.*/
		I10=0;                  /* I10 is used to increment the */
		L10=0;                  /* frequency of sine lookup.*/
		F15=0;                  /* Zero to clear accumulators */
		LCNTR=N, DO outer UNTIL LCE;
		  F8=PASS F15, M8=I10;          /* Update frequency */
		  F9=PASS F15, F0=DM(I0,M1), F5=PM(I9,M8);
		  F12=F0*F5, F4=PM(I8,M8);
		  LCNTR=N-1, DO inner UNTIL LCE;
		    F13=F0*F4, F9=F9+F12, F0=DM(I0,M1), F5=PM(I9,M8);
inner:              F12=F0*F5, F8=F8-F13, F4=PM(I8,M8);
		  F13=F0*F4, F9=F9+F12;
		  F8=F8-F13, DM(I2,M1)=F9;      /* Write real result */
		  MODIFY(I10,M9);               /* Increment frequency */
outer:            DM(I1,M1)=F8;                 /* Write imaginary result */
		RTS;