/*_______________________________________________________________________
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]= "noise64.dat";
.VAR real[N]= "zeros64.dat";
.VAR imag[N]= "zeros64.dat";


.SECTION/PM     pm_data;        /* Declare variables in program memory */
.VAR sine[N]= "sin64.dat";      /* Cosine is derived using a shifted */ 
.VAR modul[N]= "zeros64.dat";

.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;                 /*           ''            */
		nop;
		call modl;
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;

/*_________________________modul computing Subroutine___________________________*/
modl:
			I0=real;
			L0=0;
			I1=imag;
			L1=0;
			I8=modul;
			L8=0;
			F2=DM(I0,M1);
			F6=F2;
			F3=DM(I1,M1);
			F7=F3;
			F8=F2*F6;
			F13=F3*F7;
			
		LCNTR=N, DO mod_cal UNTIL LCE;
			F8=F2*F6,F13=F8+F13,F3=DM(I1,M1);
			F4=RSQRTS F13, F7=F3;
			F4=F4*F13,F2=DM(I0,M1);
			F13=F3*F7,F6=F2;
mod_cal:		PM(I8,M9)=F4;
			nop;
			nop;
			RTS;