In fact, imho, it should also be a basic function in every modern programming app.
The DFT (Discrete Fourier Transform) is the basic algorithm and can be used for short series of data.
For longer series, the time required increases quadratically and the algorithm becomes unsuitable for use.
In 1965 a faster algorithm was developed by Cooley and Tukey.
"The motivation for it [FFT algorithm] was provided by Dr. Richard L. Garwin at IBM Watson Research who was concerned about verifying a Nuclear arms treaty with the Soviet Union for the SALT talks. Garwin thought that if he had a very much faster Fourier Transform he could plant sensors in the ground in countries surrounding the Soviet Union. He suggested the idea of how Fourier transforms could be programmed to be much faster to both Cooley and Tukey. They did the work, the sensors were planted, and he was able to locate nuclear explosions to within 15 kilometers of where they were occurring" (Wikipedia article on James Cooley).
I have translated some functions to BestBASIC:
• The straight DFT in two versions, with start index at 0 and 1
• The FFT according to the algorithm in Wikipedia; "Cooley – Tukey FFT algorithm". Start index is 0
• The FFT coded by "Henko" for SmartBasic and derived from a Fortran program. Start index is 1
See https://kibernetik.pro/forum/viewtopic.php?f=20&t=1587s
Code: Select all
'Test DFT & FFT versions
'by Ton Nillesen, October 2019
'
'Two versions of FFT:
'FFT according to Wikipedia "Cooley–Tukey FFT algorithm"
'and Henko's version in his SmartBasic DSP-utility
'
'Example according to "Rosetta's code"
'at http://www.rosettacode.org/wiki/Fast_Fourier_transform
'
'--- select function DFTx
DFT=0 'Select DFT-function: DFT0=0, DFT1=1
'
'---- Constants
SCRSIZE 575 605
OPTION BASE 0 ' Wiki_FFT requires base 0
OPTION TAB 10
pwr=3
size=2^pwr ' size of FFT
'--- show both versions
PRINT "Input data"
FOR mode=1 TO 3
DIM A(size)
RESTORE
'--- fill array with test-data
first=0 : last=size-1
IF (mode=3) OR (mode=1 AND DFT=1) THEN
first+=1 : last+=1
ENDIF
FOR i=first TO last
READ A(i): IF mode=1 THEN PRINT " ";A(i)
NEXT i
PRINT
'--- Transform and show
IF mode=1
IF DFT=0 THEN A=DFT0(A,size) ELSE A=DFT1(A,size)
PRINT "Result of straight DFT"
ELSE IF mode=2 THEN
A=Wiki_FFT(A,size)
PRINT "FFT result (Wikipedia algorithm)"
ELSE IF mode=3 THEN
A=Henko_FFT(A,size)
PRINT "FFT result of Henko's version"
ENDIF
PRINT "Real","Imag"
FOR i=first TO last
PRINT REAL(A(i)),IMAG(A(i))
NEXT i
NEXT mode
END
'
DATA 1,1,1,1,0,0,0,0 'Rosetta's example of input
'
DEF DFT0(A,n) 'straight DFT base 0
' Input: Array A with complex or real values A(0)-A(n-1)
' Return value: array with the DFT
IF n<2 THEN RETURN
DIM B(n)
FOR k=0 TO n-1
q=-1i*2*PI*k/n : T=0
FOR i=0 TO n-1 : T+=A(i)*EXP(q*i) : NEXT i
B(k)=T
NEXT k
RETURN B
END DEF
'
DEF DFT1(A,n) 'straight DFT base 1
' Input: Array A with complex or real values A(1)-A(n)
' Return value: array with the DFT
IF n<2 THEN RETURN
DIM B(n)
FOR k=1 TO n
q=-1i*2*PI*(k-1)/n : T=0
FOR i=0 TO n-1 : T+=A(i+1)*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
- Test DFT&FFT's screenshot.JPG (41.14 KiB) Viewed 336 times
The example is according to the one given in "Rosetta's code" at http://www.rosettacode.org/wiki/Fast_Fourier_transform
A speed test will follow
