Page 1 of 1

### Speedtest on DFT and FFT's

Posted: 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 (28.95 KiB) Viewed 393 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.