Current issue

Vol.26 No.1

Vol.26 No.1

Volumes

© 1984-2014
British APL Association
All rights reserved.

Archive articles posted online on request: ask the archivist.

archive/16/3

Volume 16, No.3

Note: This article contains figures which may take some time to download (the largest are 190K files). Since the details in the pictures are important we did not feel we should compress them further. Ed

Fast Fourier Transforms, Diffraction Patterns, and J

by Ned W. Allis (allisn@lafayette.edu), Jeffrey P. Dumont (dumontj@lafayette.edu),
Flynn J. Heiss (flynnheiss@yahoo.com), Clifford A. Reiter (reiterc@lafayette.edu)

The fast Fourier transform package available for J is very useful for simulating diffraction patterns. Here we use it to replicate a classic optical transform and explore its power in creating diffraction patterns of fractal images and chaotic attractors. These types of images are already visually stimulating, and we found that the diffraction patterns could also yield remarkable images. We will define some J functions which are helpful for viewing the results of the J Fourier transform package.

The diffraction patterns of crystals have been the focus of research since the discovery of x-ray diffraction in 1912. Diffraction patterns and their relationship to Fourier transforms are discussed in detail in [6]. Renewed interest in diffraction patterns has arisen due to the discovery of quasicrystals which exhibit classically forbidden symmetry in their diffraction patterns. Senechal gives a good introduction to quasicrystals [11].

The Fast Fourier Script in J:

The fast Fourier transform package can be obtained from the J website [4]. Our J examples assume that the fftw add-on has been installed and the fftw.ijs script has been loaded. The package is based on the code written originally by Frigo and Johnson [2]. The script includes both the forward and reverse transforms as predefined functions fftw and ifftw. A detailed description of the mathematics involved can be found at [2], but here we offer a basic explanation. The fftw function defined in the script takes a complex-valued array X of any dimension d as its right hand argument and returns another d-dimensional complex-valued array Y according to the following definition: equation ,

where equation, for all equation.
The reverse transform is computed using:
equation.
Notice that the only difference between the forward and the reverse transforms is the minus signs within the exponents. We now give a simple example of the use of the fftw function.

   load'c:\j403\system\packages\fftw\fftw.ijs'
   ]m=:2 2$1 2 1j1 3j1
  1   2
1j1 3j1
   fftw m
  7j2 _3
_1j_2  1

Computations using fftw, although rapid, can be very memory intensive. Most of the input arrays we work with are approximately 1000 by 1000 pixels. The following function

   ts=.6!:2 ,7!:2@:]

measures time and space usage. In Table 1, we list the time and memory that it took to perform an fftw operation on our 400 MHz PC. We note that the transform is significantly faster on arrays whose size is a power of 2, although the transform on a slightly larger array uses the same amount of memory. It also seems that the limit on our system of fftw on square matrices is 2047 by 2047.

Size Time Memory
64 by 64 0.016 sec 1.42 megs
128 by 128 0.062 sec 5.65 megs
256 by 256 0.390 sec 22.55 megs
512 by 512 1.578 sec 90.19 megs
513 by 513 3.141 sec 90.19 megs
1024 by 1024 9.531 sec 360.725 megs
1920 by 1920 30.74 sec 721.435 megs
2047 by 2047 51.656 sec 721.435 megs
Table 1: Time and memory involved in fftw computations

Diffraction Patterns

As a first application of fftw we replicate one of the experiments in the Atlas of Optical Transforms [3]. The input, seen in Figure 1, consists of three circular apertures arranged horizontally, and can be found on the left side of plate 10 in the Atlas of Optical Transforms. The corresponding diffraction pattern can be found in the same position on the right side of plate 10. To recreate that input, we use

   a=:100 100$0
   b=:_100 _100{.66 66{.255*25>+/"1*:,"0/~i:15
   B=:_700 _700{.500 500{.(a,.a,.a) ([,],[)b,.b,.b

Here, a gives a 100 by 100 array of 0’s, b gives a circular region of 255’s embedded in a 100 by 100 array of zeros, and B gives us the desired array displayed in Figure 1.

Figure 1
Figure 1: Three apertures

We then compute the magnitude of the complex-valued Fourier transform. Next we use the dis function to make these values discrete in order to create an 8-bit greyscale bitmap.

   max=:>./@,@]
   min=:<./@,@]
   maxout=:0.9999&*@[
   dis=:[:<.(maxout*]-min)%max-min

Here we demonstrate the use of this function:

   x=:(1+i.20)%19r5p1 
   Sin=:1&o.
   data=:(Sin x)%x
   10 dis data
9 9 9 9 9 9 8 8 7 7 6 6 5 4 4 3 2 1 0 0

Applying the dis function to the magnitude of our transform, we write out the bitmap using writebmp8 as defined in raster3.ijs [10] using a greyscale palette, g256.

   tr=:256 dis |fftw B
   g256=:3#"0 i.256
   (g256;tr) writebmp8 'figure2.bmp'

Figure 2 almost gives us the expected diffraction pattern. It turns out that the corners of our images look similar to what we should see as the centre and vice versa.

Figure 2
Figure 2: The magnitude of the Fourier transform of Figure 1
as given by the fftw function with dis

Thus, we use the following functions to re-centre our image.

   hr=: <.@-:@# |. ]
   HR=:hr"1@:hr

The function hr rotates an array half way along its length while HR applies hr to both axes of a matrix. Applying HR to the data shown in Figure 2 yields Figure 3, which closely resembles the expected diffraction pattern found in the Atlas of Optical Transforms; while this image seems distressingly black, it actually shows slightly more contrast than the image in the atlas. In fact, we can see an extra ring.

Figure 3
Figure 3: Figure 2 as recentered by the HR function

At this point, we maximise the contrast to better understand the diffraction pattern. This leads us to the use of the following function from [8]:

   cile=:$@] $ ((/:@/:@] <.@:* (%#)),)

Here we demonstrate the use of the cile function:

   10 cile data
9 9 8 8 7 7 6 6 5 5 4 4 3 3 2 2 1 1 0 0

Notice the difference between the outputs of the cile and dis functions applied on the same data. The output of the cile function distributes the data into groups of 2 while the dis function has an abundance of 9s but only one 2. Using this function on the output of the Fourier transform, HR 256 cile |fftw B, provides us a detailed look at the fringes in the diffraction pattern (see Figure 4).

Figure 4
Figure 4: The magnitude of the Fourier transform of Figure 1 as
given by the cile function and recentered according to the HR function

Transforms of Sierpinski Fractals

Other visually attractive diffraction patterns are created by using Sierpinski fractals as input. We created both the Sierpinski triangle and carpet with the following function given in [7].

   ggl=:+./"1@([(+./@(+./)@(</~@i.
      @[*.*."1/~@]))"2>@>@{@(#<@(<"1@])))

Figure 5 displays a Sierpinski triangle on the left created with

   st=:2 ggl (n#3)#:i.2^n=.9

Figure 5
Figure 5: A Sierpinski triangle and the magnitude of its transform

Its Fourier transform is located on the right side in Figure 5. Here we maximise contrast with the cile function. The palette is black at one end, white at the other end, and blends between red and magenta in the middle. Note that the Fourier transform, like the input array, displays self-similarity, a common feature of fractals. Figure 6 displays a Sierpinski carpet on the left created by

   sc=:2 ggl (n#3)#:i.3^n=.5

and shows its Fourier transform on the right, which also displays fractal-like qualities.

Figure 6
Figure 6: A Sierpinski carpet and the magnitude of its transform

To give an example of a 3-dimensional Fourier transform, we create a Menger Sponge as follows.

   ms=:3 ggl (n#3)#:i.3^n=.4

We then take the Fourier transform of the sponge and display each plane of both the Menger Sponge and the transform in an animation which may be found at [9].

Transforms of Chaotic Attractors

We will see that chaotic attractors and their diffraction patterns may be visually appealing. We start with an attractor with rotational symmetry. Such images have diffraction patterns that also contain rotational symmetry which is either equal to or twice that of the original image. For example, the left side of Figure 7 displays a chaotic attractor with 5-fold rotational symmetry, which was created by the methods discussed in [5]. Its transform, located on the right, has 10-fold symmetry, arising from the centrosymmetric property of the Fourier transform [11] combined with the original symmetry type carried over by the transform.

Figure 7
Figure 7: A chaotic attractor with dihedral 5-fold symmetry
and the magnitude of its transform

In search of quasicrystal-like chaotic attractors, we also tried taking transforms of chaotic attractors with almost forbidden symmetry created with techniques discussed in [1]. Figure 8 gives an example of a chaotic attractor with almost pentagonal local symmetry and its transform, which consists of points on a square lattice. These points are so focused that they need to be enlarged. We use the following J function to make them more visible.

   tess=:4 4&((>./@,);._3

This results in a 7-fold magnification of the lattice points on the right side of Figure 8.

Figure 8
Figure 8: A chaotic attractor with forbidden symmetry
and the magnitude of its transform

Acknowledgement: This work was supported in part by NSF grant DMS-9805507. The chaotic attractor in Figure 7 was joint work with a co-author from [5].

References

  1. Dumont, J. and Reiter C. Chaotic Attractors near Forbidden Symmetry, Chaos Solitons and Fractals, to appear
  2. Frigo, M. and Johnson, S. FFTW Homepage, http://www.fftw.org
  3. Harburn, G. , Taylor, C. A. and Welberry, T. R. Atlas of Optical Transforms, G. Bells & Sons, Ltd, (1972)
  4. J Software, http://www.jsoftware.com
  5. Jones, K. and Reiter, C. A. Chaotic Attractors with Cyclic Symmetry Revisited, Computers and Graphics, to appear
  6. Lipson, S. G. , Lipson H. and Tannhauser D. S., Optical Physics, Cambridge University Press, (1995)
  7. Reiter, C. A. Sierpinski Fractals and GCDs, Computers and Graphics, 18 (1994) 885–891
  8. Reiter, C. A. Fractals, Visualization, and J, Iverson Software, (1995)
  9. Reiter, C. A. http://www.lafayette.edu/~reiterc/fft/
  10. Reiter, C. A. Reiter, C. A. http://www.lafayette.edu/~reiterc/j
  11. Senechal, M. Quasicrystals and Geometry, Cambridge University Press, (1995)
Ned W. Allis
Lafayette College, Box 8210
Easton, Pa 18042, U.S.A.
Jeffrey P. Dumont
Lafayette College, Box 7756
Easton, Pa 18042, U.S.A.
Flynn J. Heiss
57 Main Street
Geneseo, NY 14454, U.S.A.
Clifford A. Reiter
Department of Mathematics, Lafayette College
Easton, Pa 18042, U.S.A.

script began 6:23:30
caching off
debug mode off
cache time 3600 sec
indmtime not found in cache
cached index is fresh
recompiling index.xml
index compiled in 0.4533 secs
read index
read issues/index.xml
identified 26 volumes, 99 issues
array (
  'id' => '10007200',
)
regenerated static HTML
article source is 'HTML'
source file encoding is ''
read as 'Windows-1252'
URL: #ref6 => art10007200#ref6
URL: #ref11 => art10007200#ref11
URL: #ref4 => art10007200#ref4
URL: #ref2 => art10007200#ref2
URL: #ref2 => art10007200#ref2
URL: cliffeq1.gif => trad/v163/cliffeq1.gif
URL: cliffeq2.gif => trad/v163/cliffeq2.gif
URL: cliffeq3.gif => trad/v163/cliffeq3.gif
URL: cliffeq4.gif => trad/v163/cliffeq4.gif
URL: #ref3 => art10007200#ref3
URL: cliff1.gif => trad/v163/cliff1.gif
URL: #ref10 => art10007200#ref10
URL: cliff2.jpg => trad/v163/cliff2.jpg
URL: cliff3.jpg => trad/v163/cliff3.jpg
URL: #ref8 => art10007200#ref8
URL: cliff4.jpg => trad/v163/cliff4.jpg
URL: #ref7 => art10007200#ref7
URL: cliff5.gif => trad/v163/cliff5.gif
URL: cliff6.gif => trad/v163/cliff6.gif
URL: #ref9 => art10007200#ref9
URL: #ref5 => art10007200#ref5
URL: #ref11 => art10007200#ref11
URL: cliff7.gif => trad/v163/cliff7.gif
URL: #ref8 => art10007200#ref8
URL: cliff8.gif => trad/v163/cliff8.gif
URL: #ref5 => art10007200#ref5
URL: http://www.fftw.org => http://www.fftw.org
URL: http://www.jsoftware.com => http://www.jsoftware.com
URL: http://www.lafayette.edu/~reiterc/fft/ => http://www.lafayette.edu/~reiterc/fft/
URL: http://www.lafayette.edu/~reiterc/j/ => http://www.lafayette.edu/~reiterc/j/
completed in 0.5259 secs