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
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
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