DFT code
The DFT is defined by
(1.1)
Relations (1.2) and (1.3)
can be used for values of i outside the N/2
to N/2
1
range. The result is that
(1.4)
The FFT uses d or D to return all N independent values of the other. It is fast owing to the fact that it does this for only a small 2-100 factor times the time required to find D or d for even a single value using (1.2) or (1.3). The DFT is mostly of pedagogical interest, but for a few values it is actually competitive.
Fortran 1
Rewrite (1.2) slightly to eliminate the time
The di’s are assumed to be in order from i = -N/2 to N/2-1
TDFT.FOR DFT.FOR
those for this are in for\dft1.zip
Test code
The test code produces a cosine that should integrate to N/2 or 0.
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (N=1024)
COMPLEX*16 D(1024),FSUM(1024)
DATA PI/3.141592653589793D0/
DATA TIME/0.37D0/
MD=21
FD=MD/TIME
IT=-N/2
DO I=1,N
T=IT*TIME/N
D(I)=COS(T*FD*2*PI)
IT=IT+1
ENDDO
M1=-512
M2=511
CALL DFTF(D,FSUM,N,M1,M2)
OPEN(1,FILE='TDFTR.OUT')
OPEN(2,FILE='TDFTI.OUT')
IM=1
DO M=M1,M2
WRITE(1,'(I7,G15.6)')M,REAL(FSUM(IM))
WRITE(2,'(I7,G15.6)')M,IMAG(FSUM(IM))
IM=IM+1
ENDDO
END
C$INCLUDE
DFT set up for Watfor
The dft using complex values
SUBROUTINE DFTF(DAT,FSUM,N,M1,M2)
IMPLICIT REAL*8 (A-H,O-Z)
COMPLEX*16 DAT(*),FSUM(*),TEMP
DATA PI/3.141592653589793D0/
IM=1
IB=-N/2
IE=IB+N-1
DO M=M1,M2
FSUM(IM)=0
IDAT=1
DO I=IB,IE
ARG=(2*PI*I*M)/N
Note that the π makes the argument in parenthesis double precision so that the division by N is not integer math.
TEMP=EXP((0,1)*ARG) complex
FSUM(IM)=FSUM(IM)+TEMP*DAT(IDAT)
IDAT=IDAT+1
ENDDO
IM=IM+1
ENDDO
RETURN
END
Figure 1 Forward transform of cos(21*2 pi * t). x axis is
m as in fm = m/T. Full range
of 512
to 511 is in the data.
Fortran 2
This code uses the cosine and sine rather than the complex exponential
SUBROUTINE DFTF(DAT,FSUM,N,M1,M2)
IMPLICIT REAL*8 (A-H,O-Z)
COMPLEX*16 DAT(*),FSUM(*),TEMP
DATA PI/3.141592653589793D0/
IM=1
IB=-N/2
IE=IB+N-1
DO M=M1,M2
FSUM(IM)=0
IDAT=1
AREAL=0
AIMAG=0
DO I=IB,IE
IMT=MOD(I*M,N)
The ability to use an integer mod is due to the argument
becoming i×m/m rather
than tifm. The sin
cannot utilize the integer mod in its operation owing to the fact that its
argument is double precision. Note that
the integers are limited to ±21474836472×109, thus for N/2 > 46,340, the
multiplication will wrap the integers.
ARG=(2*PI*IMT)/N
CT=COS(ARG)
ST=SIN(ARG)
FSUM(IM)=FSUM(IM)+(
IDAT=IDAT+1
ENDDO
IM=IM+1
ENDDO
RETURN
END
This had T=.37 and N=1024
Slight modifications in test code eliminating the complex exponentials for\dft2.zip
Assignment 
The integer multiplication of i*m is exact for i and m < 46,000. It is followed by a mod function limiting the product to N in value. Write
FUNCTION MOD(I,M,N)
This function should work for I and M up to 2×109.
C
The C function modf is somewhat more complicated than the Fortran mod. Since an integer mod is extremely simple, it is included in cpp\Cmod.zip
Tmod
#include <stdio.h>
int mod(int itot,int idiv);
void main(void)
{int itot,idiv,irem;
redbeg:
printf(" enter itot,idiv >");
fflush(stdin);
scanf("%d %d",&itot,&idiv);
if(idiv == 0)return;
irem = mod(itot,idiv);
printf(" irem = %d ",irem);
goto redbeg;}
The mod routine
int mod(int itot,int idiv)
{int it,iret;
it=itot/idiv;
iret = itot - idiv*it;
return iret;}
Assignment C_mod
The integer multiplication of i*m is exact for i and m < 46,000. It is followed by a mod function limiting the product to N in value. Write
int mod(int i,int m, int n)
This function should work for |i| and |m| < 2×109.
Cdft
The C version of the dft is in cpp\cdft.zip.
The test code
#include <stdio.h>
#include <math.h>
int mod(int itot,int idiv);
void dft(double datai[],double datao[],int nn,int m1,int m2);
void main(void)
{double d[2048],fsum[2048];
/* real values are 0,2,4,...,2046. imag values are 1,3,5,...,2047 */
double pi=3.141592653589793;
double arg;
int md=21,m1=-30,m2=30;
int n=1024;
int i,it,itm,m,im;
FILE *fpr,*fpi;
it=-n/2;
for(i=0;i<2*n;i+=2)
{itm=mod(it*md,n);
arg=itm*2*pi/n;
d[i]=sin(arg);
d[i+1]=0;
it+=1;}
dft(d,fsum,n,m1,m2);
fpr=fopen("tdftr.out","w");
fpi=fopen("tdfti.out","w");
im=0;
for(m=m1;m<=m2;m+=1)
{fprintf(fpr," %d %lf ",m,fsum[im]);
fprintf(fpi," %d %lf ",m,fsum[im+1]);
im+=2;}
return;}
The dft.c code
#include <math.h>
#include <stdio.h>
int mod(int itot,int idiv);
void dft(double datai[],double datao[],int n,int m1,int m2)
{double pi=3.141592653589793,arg,ct,st;
int i,it,m,mt,imt,itest;
mt=0;
for (m=m1;m<=m2;m+=1)
{datao[mt]=0;
datao[mt+1]=0;
it=0;
for (i=-n/2;i<n/2;i+=1)
{imt=mod(i*m,n);
arg=2*pi*imt/n;
ct=cos(arg);
st=sin(arg);
datao[mt]+=datai[it]*ct-datai[it+1]*st;
datao[mt+1]+=datai[it+1]*ct+datai[it]*st;
it+=2;}
mt+=2;}
return;}
Re-invention of radio
Figure 2 This is data with a Gaussian times the relevant frequency plus an irrelevant frequency.
MD=128
MD2=140
FD=MD/TIME
IT=-N/2
DO I=1,N
T=TIME*IT/N
ITM=MOD(IT*MD,N)
ITM2=MOD(IT*MD2,N)
ARG=ITM*2*PI/N
ARG2=ITM2*2*PI/N
DT(I)=COS(ARG)*EXP(-T*T*200)+
Figure 3 Data above after frequency shifting and filtering.
c Notch filtering
IM=1
DO M=-N/2,N/2-1
IF(ABS(M).GE.1)THEN
DF(IM)=0
ENDIF
IM=IM+1
ENDDO
c Error Function Conjugate Filtering
F=1D0*M/T
H=AIGAU((F+FH)/W)-AIGAU((F-FH)/W)
DF(IM)=DF(IM)*H