{_______________________________________________________________________
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)

Author: 27-APR-91, Steven Cox, Analog Devices DSP Div. (617) 461-3672
_________________________________________________________________________}
{ Constant for number of points in input }
#define N 64

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

.SEGMENT/PM	pm_data;	{ Sine array is N points of one cycle }
.VAR sine[N]= "sine.dat";	{ Cosine is derived using a shifted } 
.ENDSEG;			{ pointer to this circular buffer.}

.SEGMENT/PM	rst_svc;	{ 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;