{_______________________________________________________________________
DFT.ASM		ADSP-21020 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)
_________________________________________________________________________}
{ Constant for number of points in input }
#define N 64

.SEGMENT/DM	dm_data;	{ Declare variables in data memory }
.VAR input[N]= "test64.dat";
.VAR real[N];
.VAR imag[N];
.ENDSEG;

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

.SEGMENT/PM	pm_rsti;	{ The reset vector resides in this space }
		DMWAIT=0x21;	{ Set data memory waitstates to zero }
		PMWAIT=0x21;	{ Set program memory waitstates to zero }
		JUMP start;
.ENDSEG;

.SEGMENT/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;
.ENDSEG;