Speedtest on DFT and FFT's

For finished programs
Post Reply
User avatar
Dutchman
Posts: 92
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

Speedtest on DFT and FFT's

Post by Dutchman » Thu Oct 10, 2019 2:19 pm

The following program does a speed test with the three versions of the Discrete Fourier Transform as published in viewtopic.php?f=4&t=41

Code: Select all

'Speedtest FFT-versions and DFT
'by Ton Nillesen, October 2019
'FFT according to Wikipedia Cooley–Tukey FFT algorithm
'  and Henko's translation of Fortran algorithm
'
'---- determine size ----
pwr=12 ' size is 2^pwr
'
'--- fill array with test-data
OPTION BASE 0
size=2^pwr
DIM A(size)
RANDOMIZE
FOR i=1 TO size
  A(i)=2*RND-1 +2i*RND-1i'complex
NEXT i
FOR mode=1 TO 3
  IF mode=3 THEN
    PRINT "Start DFT on ";size;" samples."
    TIMER SET : CALL DFT(A,size)
  ELSE IF mode=2 THEN
    PRINT "Start wiki_FFT on ";size;" samples."
    TIMER SET : CALL Wiki_FFT(A,size)
  ELSE IF mode=1 THEN
    PRINT "Start Henko_FFT on ";size;" samples."
    TIMER SET : CALL henko_FFT(A,size)
  ENDIF
  t=TIMER GET
  PRINT "Finished in "+STRTEXT (t, "0.00")+" seconds"
  PRINT
NEXT mode
END

DATA 1,1,1,1,0,0,0,0

'
DEF DFT(A,n) 'straight DFT
IF n<2 THEN RETURN
DIM B(n)
FOR k=0 TO n-1
  q=-1i*2*PI*k/n : T=0 ' improve speed
  FOR i=0 TO n-1
    T+=A(i)*EXP(q*i)
  NEXT i
  B(k)=T
NEXT k
RETURN B
END DEF
'
DEF Wiki_FFT(A,n) 'Wikipedia: Cooley–Tukey FFT algorithm
' Input: Array A with complex or real values A(0)-A(n-1)
' size n should be power of 2
' Return value: array with the DFT
' Algorithm first sorts bit-reverse, then transforms
'--- check parameters
pwr=INT(log(n)/log(2))
IF 2^pwr<>n
  PRINT "ERROR in FFT array-size"
  STOP
ENDIF
IF OPTIONBASE()<>0
  PRINT "ERROR: option base should be zero for FFT"
  STOP
ENDIF
'--- sort bit-reverse. Translated from FORTRAN-program
nd2=N/2 : j=nd2
FOR i=1 TO N-2
  IF i<j : tr=A(j) : A(j)=A(i) : A(i)=tr : ENDIF
  k=nd2
  DO IF k<=j : j-=k : k/=2 : REDO
  j+=k
NEXT i
'--- transform
FOR s=1 TO pwr
  m=2^s
  Pm=exp(-2*1i*pi/m)
  FOR k=0 TO n-1 STEP m
    P=1
    FOR j=0 TO m/2-1
      t=P*A(k+j+m/2)
      u=A(k+j)
      A(k+j)=u+t
      A(k+j+m/2)=u-t
      P=P*Pm
    NEXT j
  NEXT k
NEXT s   
RETURN A
END DEF
'
DEF Henko_FFT(A n)
' Based on Henko's translation of FORTRAN-program
' Input: Array A() with complex values A(1) - A(n)
' n should be power of 2
' Return value: array with the DFT
' speed is two times faster then Wikipedia's Cooley–Tukey FFT algorithm
' Algorithm first transforms, then sorts bit-reverse
'--- check parameters
m=INT(log(n)/log(2))
IF 2^m<>n
  PRINT "ERROR in FFT array-size"
  STOP
ENDIF
'--- Transform
FOR k=1 TO M              'loop for each stage
  ke=2^(M-k+1) : ke2=ke/2
  u=1+0i : angl=pi/ke2
  's=COS(angl)-SIN(angl)*1i 
  s=EXP(-1i*angl)
  FOR j=1 TO ke2          'loop for each sub-DFT
    FOR i=j TO N STEP ke  'loop for each butterfly
      ip=i+ke2 : t=A(i)+A(ip)
      A(ip)=u*(A(i)-A(ip)) : A(i)=t
    NEXT i
    u*=s
  NEXT j
NEXT k
'--- Sort bit-reverse
nd2=N/2 : nm1=N-1 : j=1
FOR i=1 TO nm1
  IF i<j : t=A(j) : A(j)=A(i) : A(i)=t : ENDIF
  k=nd2
  DO IF k<j : j-=k : k/=2 : REDO
  j+=k
NEXT i
RETURN A
END DEF
.
For a long-term test I set the variable 'pwr' to 18, which gives a size of 262144 samples.
.
Speedtest screenshot.JPG
Speedtest screenshot.JPG (28.95 KiB) Viewed 45 times
.
This test shows clearly
• "Henko" 's version of the FFT is twice as fast as the Wikipedia algorithm.
• The straight DFT is unusable for such large data series due to the hours required for a single transformation.
It is still a long way to go

Post Reply