DFT and FFT's

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

DFT and FFT's

Post by Dutchman » Thu Oct 10, 2019 1:36 pm

Fourier transform is a basic function in digital signal processing (DSP).
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
Test DFT&FFT's screenshot.JPG (41.14 KiB) Viewed 90 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
It is still a long way to go

User avatar
kibernetik
Site Admin
Posts: 103
Joined: Tue Aug 06, 2019 3:03 pm

Re: DFT and FFT's

Post by kibernetik » Thu Oct 10, 2019 1:40 pm

Yes, in-app implementation of FFT function is planned.

User avatar
Dutchman
Posts: 92
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

Re: DFT and FFT's

Post by Dutchman » Thu Oct 10, 2019 1:58 pm

Great :idea:
Please add also the DFT, which is attractive for short series because it doesn't have the requirement for the size to be a power of two.
Adding zero's to the data series in order to fit the size, corrupts the data.
It is still a long way to go

User avatar
kibernetik
Site Admin
Posts: 103
Joined: Tue Aug 06, 2019 3:03 pm

Re: DFT and FFT's

Post by kibernetik » Mon Oct 14, 2019 5:39 pm

As version 1.5 (featuring built-in FFT) is already published, now this FFT program could look like:

Code: Select all

dim m()=1,1,1,1,0,0,0,0
m = math fft m
for i = 0 to 7
  print m(i)
next
With output:

Code: Select all

 4 
 1-2.41421i 
 0 
 1-0.414214i 
 0 
 1+0.414214i 
 0 
 1+2.41421i 

Post Reply