2-Dimensional Diffraction- FFT in 2D
In my student days we didn't have lasers to give coherent monochromatis light. It was still magic to be able to flip from frequency domain to real world. The beauty of the images, and the rapidly developing application of X-rays to crystallography were very attractive to me. A year or two later and I built a mains powered He-Ne gas laser. Now we can buy high-intensity solid-state laser cheaply on the 'net.
The screen shots below show the analysis and re-synthesis of four waves.
Beats shows 'sin 30 +sin 33'. Pulse is a low mark/space rectangular wave. Square wave approximation was generated from the known series of terms.
. . 
. . 
Example LB code- uses a supplied colour-scale.
' ***********************************************************
' ** **
' ** 1DFFT.bas 27/02/2017 tenochtitlanuk **
' ** beats **
' ***********************************************************
global n, pi ' globally available pi, number of data values- a multiple of 2-
n =1024
' globally available real and imaginary data for main and sub sections
dim realData( n), imagData( n) ' the FFT operates in-place, changing the data to its transform.
pi = 4 *atn( 1)
WindowWidth =1080
WindowHeight = 720
nomainwin
open "FFT demo." for graphics_nsb as #wg
#wg "trapclose quit"
#wg "goto 20 100 ; down ; size 1 ; fill cyan ; color white"
for i =1 to n
'if i >508 and i <516 then realData( i) =1 else realData( i) =0.05
for k = 1 to 3 step 2
realData( i) =realData( i) +1 /k *sin( i /n *k *10 *2 *pi)
next k
imagData( i) =0
#wg "color red ; line "; 20 +1024 *i /n; " 80 "; 20 +1024 *i /n; " "; 80 -int( 25 *realData( i))
next i
call fft -1 ' ie compute FFT not inverse- 1 for fwd FFT, -1 for reverse
#wg "size 5"
for i =1 to n ' but symmetrical on x axis so really only want left-hand half.
#wg "color red ; line "; 20 +1024 *i /n; " 250 "; 20 +1024 *i /n; " "; 250 -int( 0.18 *realData( i))
#wg "color darkblue ; line "; 20 +1024 *i /n; " 350 "; 20 +1024 *i /n; " "; 350 -int( 0.18 *imagData( i))
next i
call fft 1 ' ie compute inverse
#wg "size 1"
for i =1 to n ' but symmetrical on x axis so really only want left-hand half.
#wg "color red ; line "; 20 +1024 *i /n; " 550 "; 20 +1024 *i /n; " "; 550 -int( 25 *realData( i))
next i
#wg "color black ; up ; goto 10 80 ; down ; goto 1034 80"
for v =250 to 550 step 100
if v <>450 then #wg "up ; goto 10 "; v; " ; down ; goto 1034 "; v
next v
wait
sub fft isi
j =1
for i =1 to n
if i <j then call exchangeTerms i, j
m =n /2
while j >m
j =j -m
m =int( ( m +1) /2)
wend
j =j +m
next i
mmax =1
while mmax <n
istep =2 *mmax
for m =1 to mmax
theta =pi *isi *( m -1) /mmax
wr =cos( theta)
wi =sin( theta)
for i =m to n step istep
j =i +mmax
tempr =wr *realData( j) -wi *imagData( j)
tempi =wr *imagData( j) +wi *realData( j)
realData( j) =realData( i) -tempr
imagData( j) =imagData( i) -tempi
realData( i) =realData( i) +tempr
imagData( i) =imagData( i) +tempi
scan
next i
next m
mmax =istep
wend
if isi =-1 then exit sub ' else normalise values
for i =1 to n
realData( i) =realData( i) /n
imagData( i) =imagData( i) /n
next i
end sub
sub exchangeTerms i, j
temp =realData( i)
realData( i) =realData( j)
realData( j) =temp
temp =imagData( i)
imagData( i) =imagData( j)
imagData( j) =temp
end sub
sub quit h$
close #h$
end
end sub
..which needs this colour-scale bmp saved in the same directory folder.
Starting with a circular aperture the result of the two series of 1D transforms is the 2D result shown. Note that to agree with the way we see diffraction the results need re-assembling by re-ordering the four quarters, as done below . . .
. . .
' ***********************************************************
' ** **
' ** 2DFourier9.bas 05/03/2017 tenochtitlanuk **
' ** beats **
' ***********************************************************
nomainwin
global n, pi ' globally a<ailable pi, number of data <alues- a multiple of 2-
n = 128
' globally a<ailable real and imaginary data for main and sub sections
dim realData( n), imagData( n) ' the FFT operates in-place, changing the data array to its transform.
dim imageArray$( n, n) ' holds R and I <alues as a cs< string
dim ColLookUp$( 256)
pi = 4 *atn( 1)
UpperLeftX = 10
UpperLeftY = 10
WindowWidth =1300
WindowHeight = 700
graphicbox #w.gb1, 20, 10, 600, 570
graphicbox #w.gb2, 640, 10, 600, 570
open "2D Fourier Transform demo." for window as #w
#w "trapclose quit"
#w.gb1 "down ; fill 40 40 180 ; size 1 ; color darkgray"
#w.gb2 "down ; fill 40 40 40 ; size 1 ; color darkgray"
loadbmp "scr", "colRange.bmp"
#w.gb1 "drawbmp scr 0 20 ; size 1"
gosub [fillColLookUp]
gosub [initImageArray]
for row =1 to n
for col =1 to n ' copy a row of data to array for computing its FFT
realData( col) =<al( word$( imageArray$( row, col), 1, ","))
imagData( col) =<al( word$( imageArray$( row, col), 2, ","))
next col
call fft -1
for col =1 to n ' do FFT and copy back to image array
imageArray$( row, col) =str$( realData( col)) +"," +str$( imagData( col))
amp =amp( ( realData( col)^2 +imagData( col)^2)^0.5) ' int( 100 *log( 1 +( realData( col)^2 +imagData( col)^2)^0.5))
#w.gb2 "backcolor "; ColLookUp$( amp)
#w.gb2 "color "; ColLookUp$( amp)
call sqDraw 2, row, col
next col
next row
call screenGrab "X"
for col =1 to n
for row =1 to n ' copy a row of data to array for computing its FFT
realData( row) =<al( word$( imageArray$( row, col), 1, ","))
imagData( row) =<al( word$( imageArray$( row, col), 2, ","))
next row
call fft -1
for row =1 to n
imageArray$( row, col) =str$( realData( row)) +"," +str$( imagData( row))
amp =amp( ( realData( row)^2 +imagData( row)^2)^0.5) ' int( 100 *log( 1 +( realData( row)^2 +imagData( row)^2)^0.5))
#w.gb2 "backcolor "; ColLookUp$( amp)
#w.gb2 "color "; ColLookUp$( amp)
call sqDraw 2, row, col
next row
next col
'call screenGrab "Y"
wait
sub screenGrab i$
for x =0 to 256 step 256
for y =0 to 256 step 256
#w.gb2 "getbmp scr "; x +50; " "; y +20; " 256 256"
bmpsa<e "scr", i$ +str$( x) +str$( y) +".bmp"
next y
next x
end sub
sub fft isi
j =1
for i =1 to n
if i
Code will be up shortly as a 2D_FFT.zip