Fourier Theory and Fourier Transforms

Discrete Fourier Transforms on very simple lattices

Create a simple lattice to upon which a fourier transform will  taken

WhiteSquare =  
        Table[1, {i, 8}, {j, 8}];
BlackSquare = Table[0, {i, 8}, {j, 8}];
Join[WhiteSquare, BlackSquare, BlackSquare] // MatrixForm

Example of construction of a  slightly larger structure

latcell  = 
        Join[
Transpose[
            Join[BlackSquare, BlackSquare, BlackSquare, 
              BlackSquare]], 
Transpose[
            Join[BlackSquare, BlackSquare, WhiteSquare, 
              BlackSquare]], 
Transpose[
            Join[BlackSquare, BlackSquare, WhiteSquare, 
              BlackSquare]], 
Transpose[
            Join[BlackSquare, BlackSquare, BlackSquare, 
              BlackSquare]]
];

Visualize structure:

ListDensityPlot[latcell, MeshStyle -> {Hue[1]}];

The following duplicates an input matrix and create (2^n) copies alligned in a column

ColumnDuplicateNsq[matrix_ ,  nlog2_]  :=  
      Nest[Join[#, #] &, matrix, nlog2]

example:

ListDensityPlot[ColumnDuplicateNsq[latcell, 2], 
      MeshStyle -> {Hue[1]}]

The following duplicates an input matrix and create (2^n) copies alligned in a row

RowDuplicateNsq[matrix_ ,  nlog2_]  := 
      Transpose[ColumnDuplicateNsq[matrix, nlog2]]

example:

ListDensityPlot[RowDuplicateNsq[latcell, 2], 
      MeshStyle -> {Hue[1]}]

Create a 256 by 256 data set

XtalData  = 
        Transpose[
          ColumnDuplicateNsq[RowDuplicateNsq[latcell, 3], 
            3]];

Function to create graphics for delayed display

DisplayLater  =  
        DisplayFunction -> Identity;
DisplayNow  =  
        DisplayFunction ->  $DisplayFunction;

ImagePlot[data_ ] :=  
      ListDensityPlot[data, Mesh -> False, ImageSize -> 144, 
        DisplayLater]

Example

XtalImage  =  ImagePlot[XtalData]

To see the data:

Show[XtalImage, DisplayNow, ImageSize -> 400]

Take Discrete Fourier Transform of constructed simple lattice

FourierData  =  Fourier[XtalData];

Create Image of Fourier Transform (Here a special color function is used (but not really explained). Red is high intensity, blue is low, gray is roughly zero)
Abs is used to get intensities of spots at each of the N×M wavevectors

FourierColorFunctionWhiteBack  :=  
        ColorFunction -> (If[# <  .05, Hue[1, 0, #], 
                Hue[ .66*(1 - #), #, 
                  1]] &);
FourierColorFunctionBlackBack  :=  
        ColorFunction -> (If[# <  .1, Hue[1, 0,  5*#], 
                Hue[ .66*(1 - #),  .75 (1 + #/3), 
                  1]] &);
FourierImagePlot[data_ ]  :=  
      ListDensityPlot[Abs[data], Mesh -> False, ImageSize -> 144, 
        FourierColorFunctionBlackBack, DisplayLater]

FourierImage  =  FourierImagePlot[FourierData]

Show original, FT, and FT or FT:
Use Chop to remove small spurious imaginary values.  Because the lattice is very perfect, the fourier spots are very small, one could use the magnification feature to see the spots

Show[
      GraphicsArray[{XtalImage, FourierImage,  
          ImagePlot[Chop[Fourier[FourierData]]]}], ImageSize -> 1000,  
      DisplayNow]

Create a function to create a defect in the lattice...

HoleFunc[data_ ,  xc_ ,  yc_,  twicew_,  twiceh_] :=  
        Module[{nrows, ncols}, nrows =  Dimensions[data][[1]]; 
          ncols  =  Dimensions[data][[2]]; 
          Table[If[
              And[Abs[j - xc]  <=  twicew, Abs[i  -  yc]  <=  twiceh], 
              0, 1], {i, nrows}, {j, ncols}]];

XtalData  = 
        Transpose[
          ColumnDuplicateNsq[RowDuplicateNsq[latcell, 3], 
            3]];

Create template for defect creation

hole  =  HoleFunc[XtalData, 28, 28, 6, 6];

Recreate Images and Redisplay

XtalData  = 
        Transpose[
          ColumnDuplicateNsq[RowDuplicateNsq[latcell, 3], 
            3]]; 
XtalData  =  hole*XtalData;
XtalImage  =  ImagePlot[XtalData];
FourierData  =  Fourier[XtalData];
FourierImage  =  FourierImagePlot[FourierData];

Note curious result that the position of the defect has changed to a symmetric position in the reconstructed image

Show[
      GraphicsArray[{XtalImage, FourierImage,  
          ImagePlot[Chop[Fourier[FourierData]]]}], ImageSize -> 1000,  
      DisplayNow]

Visualization of the Fourier Transform of a lattice that has a little "thermal" noise in the lattice positions

Function to make a square with a specified size, with a lattive composed of lattice vectors:
MakeLattice[Width, Height, {(Subscript[a, x]), (Subscript[a, y]), repeats}, {(Subscript[b, x]), (Subscript[b, y]), repeats}, AtomSize, {noise_back, noise_forward}]
Function is not ideal, works best when the lattice vectors are perfect divisors of the width and height

MakeLattice[W_,  H_,  latvecA_,  
        latvecB_,  size_ ,  randrange_]  :=  
      Module[{result = Table[0, {i, H}, {j, W}], lata = -1, 
          latb = -1, xpos, 
          ypos,   untouched  =  
            Table[True, {i, H}, {j, W}]}, 
For[
          lata  =  0,  
          lata  ≤  latvecA[[3]], 
For[
            latb = 0, latb ≤  latvecB[[3]], 
            xpos = Mod[lata*latvecA[[1]]  +  latb*latvecB[[1]], H, 
                1] ; ypos  =  
              Mod[lata*latvecA[[2]]  +  latb*latvecB[[2]], W, 1] ; 
            If[untouched[[ypos, xpos]], 
              untouched[[ypos, xpos]] = False; 
              xpos  +=  Random[Integer, randrange];  
              ypos  +=  Random[Integer, randrange]; 
For[
                j  =  1 , j ≤  size, 
                For[i  =  1,  i  ≤  size, 
                  result[[Mod[ypos + j, H, 1], Mod[xpos + i, W, 1]]] = 
                    1; i++]; j++]]; 
latb++]; 
lata++]; result
]

Example of a lattice with no noise:

latdata =  
      MakeLattice[400, 400, {0, 20, 40}, {16, 4, 25}, 4, {0, 0}];  
    fourlat  =  Fourier[latdata];

Show[
      GraphicsArray[{ImagePlot[latdata], FourierImagePlot[fourlat],  
          ImagePlot[Chop[Fourier[fourlat]]]}], ImageSize -> 1000,  
      DisplayNow]

Make identical lattice, but add a little noise to the system:

thermallatdata =  
      MakeLattice[400, 400, {0, 20, 40}, {16, 4, 25}, 4, {-2, 2}];  
    thermalfourlat  =  Fourier[thermallatdata];

Visualization of the original image, its fourier transform, and the fourier transform of the fourier transform..

Show[
      GraphicsArray[{ImagePlot[thermallatdata], 
          FourierImagePlot[thermalfourlat],  
          ImagePlot[Chop[Fourier[thermalfourlat]]]}], 
      ImageSize -> 1000,  DisplayNow]

Using an Aperature to look at a particular region of reciprocal space and visualizing its effect

The following is a fairly baroque function to do something what is conceptually straightforward:
    The function takes original data from a lattice and its fourier transform and graphically compares those to data from a "noised-up" lattice. 
        The function allows the user to specify the center of the aperature in reciprocal space as well as (twice) the aperature width and height.
        The function will display eight images in two columns. The left column of graphics illustrates (from top to bottom)the "clean" input image, the entire fourier transform with the rectangular aperature illustrated, the "reconstructed image" that derives from the fourier transform of the aperature region, and finally a magnified image of the fourier transform within the aperature only.
    The right column is the same sequence of images for the "noised-up" initial data

\(Compare[sharpdata_\ , \ sharpfourierdata_\ , \ diffusedata_, \ 
        diffusefourierdata_, \ ApCenterx_\ , 
        ApCentery_, \ \ ApTwicewidth_\ , \ ApTwiceheight_\ ]\  := \ 
      Module[{sharpfourieraperature, diffusefourieraperature, \ 
          dims\  = \ Dimensions[sharpdata], nrows, ncols, sharpfourimage, 
          sharpapimage, diffusefourierimage, \ diffuseapimage, \ 
          sharprevfourimage, diffusefourimage, \ sharpfourmagimage, 
          diffusefourmagimage, aperature, xll, yll, xur, 
          yur}, 
nrows = dims[\([1]\)]; \ 
        ncols = \ 
          dims[\([2]\)]; 

sharpfourimage = \ ListDensityPlot[Abs[sharpfourierdata], Mesh -> False, 
            FourierColorFunctionBlackBack, 
 (*\(ColorFunctionScaling -> 
                False\)\(,\)*) 
DisplayLater]; 
diffusefourimage = \ 
          ListDensityPlot[Abs[diffusefourierdata], Mesh -> False, 
            FourierColorFunctionBlackBack, 
 (*\(ColorFunctionScaling -> 
                False\)\(,\)*) 
DisplayLater]; 

xll\  = \ 
          ApCenterx - ApTwicewidth; 
yll = \ 
          ApCentery\  - \ ApTwiceheight; 
xur\  = \ 
          ApCenterx\  + \ ApTwicewidth; 
yur\  = \ 
          ApCentery\  + \ 
            ApTwicewidth; 

If[
          And[xll\  < \ 1, yll < 1, xur > ncols, yur > nrows], 
          sharpfourieraperature = sharpfourierdata; 
          diffusefourieraperature = 
            diffusefourierdata, 
aperature = 
            Table[If[
                And[Abs[i - ApCentery]\  <= \ ApTwiceheight, 
                  Abs[j\  - \ ApCenterx]\  <= \ ApTwicewidth], 1, 0], {i, 
                nrows}, {j, 
                ncols}]; 
sharpfourieraperature = 
            aperature*sharpfourierdata; 
          diffusefourieraperature = 
            aperature*
              diffusefourierdata]; 


sharpapimage\  = \ 
          Show[sharpfourimage, 
            Graphics[{
 (*Hue[ .5,  .5, 1]*) \ 
                Hue[ .1667, 1, 1], 
                Thickness[
                  2/nrows], 
Line[{{xll - 1, 
                      yll - 1}, {xur + 1, yll - 1}, {xur + 1, \ 
                      yur + 1}, {xll - 1, yur + 1}, {xll - 1, 
                      yll - 1}}]}]]; 
diffuseapimage\  = \ 
          Show[diffusefourimage, 
            Graphics[{
 (*Hue[ .5,  .5, 1]*) \ 
                Hue[ .1667, 1, 1], 
                Thickness[
                  2/nrows], 
Line[{{xll - 1, 
                      yll - 1}, {xur + 1, yll - 1}, {xur + 1, \ 
                      yur + 1}, {xll - 1, yur + 1}, {xll - 1, 
                      yll - 1}}]}]]; 

xll\  = \ If[xll < 1, 1, xll]; \ 
        xur\  = \ If[xur > nrows, nrows, xur]; 
yll\  = \ 
          If[yll < 1, 1, yll]; \ 
        yur\  = If[yur > ncols, ncols, 
            yur]; 

sharpfourmagimage = \ 
          ListDensityPlot[
            Abs[sharpfourierdata[\([Range[yll, yur], Range[xll, xur]]\)]], 
            Mesh -> False, 
            FourierColorFunctionBlackBack,  (*\(ColorFunctionScaling -> 
                False\)\(,\)*) \ 
            DisplayLater]; 
diffusefourmagimage = \ 
          ListDensityPlot[
            Abs[diffusefourierdata[\([Range[yll, yur], Range[xll, xur]]\)]], 
            Mesh -> False, 
            FourierColorFunctionBlackBack,  (*\(ColorFunctionScaling -> 
                False\)\(,\)*) \ 
            DisplayLater]; 

sharprevfourimage = \ 
          ListDensityPlot[Abs[Chop[Fourier[sharpfourieraperature]]], 
            Mesh -> False, 
            DisplayLater]; 
diffuserevfourimage = \ 
          ListDensityPlot[Abs[Chop[Fourier[diffusefourieraperature]]], 
            Mesh -> False, 
            DisplayLater]; 

Show[
          GraphicsArray[{
{ImagePlot[sharpdata], 
                ImagePlot[diffusedata]}, 
{sharpapimage, 
                diffuseapimage}, 
{sharprevfourimage, \ 
                diffuserevfourimage}, 
{sharpfourmagimage, 
                diffusefourmagimage}
}
], 
          ImageSize -> 1000, GraphicsSpacing -> { .001,  .0}, 
          DisplayNow]; 
Print["data is ", \ 
          nrows, \ " wide, and ", \ 
          ncols, \ " high
", \ "Aperature is x∈(", 
          xll, ",", \ xur, ")  y∈(", yll, \ ",", \ 
          yur, ")"];]\)

An example, use MakeLattice to produce a "clean" and a "noised-up" lattice

latdata =  
        MakeLattice[400, 400, {0, 20, 40}, {16, 4, 25}, 
          4, {0, 0}]; 
fourlat  =  Fourier[latdata];

thermallatdata =  
      MakeLattice[400, 400, {0, 20, 40}, {16, 4, 25}, 4, {-1, 1}];  
    thermalfourlat  =  Fourier[thermallatdata];

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 50, 50, 50, 
      50]

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 350, 350, 50, 
      50]

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 20, 30, 15, 
      15]

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 30, 30, 15, 
      15]

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 35, 25, 15, 
      15]

This suggests the (potentially) interesting exercise of modifying compare to take two aperature specifications to "pick out" other periodicities in the lattice--i.e., something like the following:
      
      Compare[cleandata_ , fouriercleandata_ , noisydata_, fouriernoisydata_,  {FirstApLowX_, FirstApHiX_ , FirstApLowY_, FirstApHiY_}, {SecondApLowX_,SecondApHiX_ , SecondApLowY_, SecondApHiY_}] :=

Consider what happens to the image when the noise in anisotropic... Modify the function MakeLattice to take two "noise arguments"

MakeLattice[W_,  H_,  latvecA_,  
        latvecB_,  size_ ,  Xrandrange_,  Yrandrange_]  :=  
      Module[{result = Table[0, {i, H}, {j, W}], lata = -1, 
          latb = -1, xpos, 
          ypos,   untouched  =  
            Table[True, {i, H}, {j, W}]}, 
For[
          lata  =  0,  
          lata  ≤  latvecA[[3]], 
For[
            latb = 0, latb ≤  latvecB[[3]], 
            xpos = Mod[lata*latvecA[[1]]  +  latb*latvecB[[1]], H, 
                1] ; ypos  =  
              Mod[lata*latvecA[[2]]  +  latb*latvecB[[2]], W, 1] ; 
            If[untouched[[ypos, xpos]], 
              untouched[[ypos, xpos]] = False; 
              xpos  +=  Random[Integer, Xrandrange];  
              ypos  +=  Random[Integer, Yrandrange]; 
For[
                j  =  1 , j ≤  size, 
                For[i  =  1,  i  ≤  size, 
                  result[[Mod[ypos + j, H, 1], Mod[xpos + i, W, 1]]] = 
                    1; i++]; j++]]; 
latb++]; 
lata++]; result
]

thermallatdata =  
      MakeLattice[400, 400, {0, 20, 40}, {16, 4, 25}, 
        4, {0, 0}, {-4, 4}];  
    thermalfourlat  =  Fourier[thermallatdata];

Compare[latdata, fourlat, thermallatdata, thermalfourlat, 200, 200, 
      200, 200]

Fourier Transforms on Images

Importing an image into Mathematica, png is some of many graphics data types that Mathematica can process.

AnImage  =  
        Import["~/3016/Mathematica/fourier_xtal_data.gif"];

Show[AnImage, DisplayNow]

The gray values of this image are stored in the (1,1) position of the image with gray values in the range (0,255)

ImageData  =  AnImage[[1, 1]]/255;

Dimensions[ImageData]

FourierImageData =  Fourier[ImageData];

Show[
      GraphicsArray[{ImagePlot[ImageData], 
          FourierImagePlot[FourierImageData],  
          ImagePlot[Chop[Fourier[FourierImageData]]]}], 
      ImageSize -> 1000,  DisplayNow]

Write a function that takes a file as input, as well as an aperature specification, then displays the image, its transform and its reverse transform of the aperature fraction of the image

\(ImageFourierAperature[imagefile_, \ Apxmin_\ , 
        Apxmax_, \ \ Apymin_\ , \ Apymax_\ ]\  := \ 
      Module[{theimage\  = \ Import[imagefile], \ dims\ , nrows, ncols, 
          fourierdata, fourimage, fourieraperature, apimage, \ fourmagimage, 
          revfourierimage, aperature, xll, yll, xur, 
          yur}, 
	
fourierdata = 
          Fourier[theimage[\([1, 1]\)]/255]; 
fourimage = \ 
          FourierImagePlot[fourierdata]; 
dims = 
          Dimensions[fourierdata]; 
nrows = 
          dims[\([1]\)]; \ 
        ncols = \ dims[\([2]\)]; 
xll\  = \ 
          Round[Apxmin*ncols]; 
yll = \ 
          Round[Apymin*nrows]; 
xur\  = \ 
          Round[Apxmax*ncols]; 
yur\  = \ 
          Round[Apymax*nrows]; 

If[
          And[xll\  ≤ \ \ 1, yll ≤ \ 1, 
            xur ≥ \ ncols, yur ≥ \ nrows], 
          fourieraperature = fourierdata, 
aperature = 
            Table[If[
                And[i\  ≥ \ xll, \ i\  ≤ \ xur, \ 
                  j\  ≥ \ yll, \ j\  ≤ \ yur]\ , 1, 
                0], {i, nrows}, {j, 
                ncols}]; 
fourieraperature = 
            aperature*
              fourierdata]; 


apimage\  = \ 
          Show[fourimage, 
            Graphics[{
\ Hue[ .1667, 1, 1], 
                Thickness[
                  2/nrows], 
Line[{{xll - 1, 
                      yll - 1}, {xur + 1, yll - 1}, {xur + 1, \ 
                      yur + 1}, {xll - 1, yur + 1}, {xll - 1, 
                      yll - 1}}]}]]; 

xll\  = \ If[xll < 1, 1, xll]; \ 
        xur\  = \ If[xur > nrows, nrows, xur]; 
yll\  = \ 
          If[yll < 1, 1, yll]; \ 
        yur\  = If[yur > ncols, ncols, 
            yur]; 

fourmagimage = \ 
          ListDensityPlot[
            Abs[fourierdata[\([Range[yll, yur], Range[xll, xur]]\)]], 
            Mesh -> False, FourierColorFunctionBlackBack, \ 
            DisplayLater]; 


revfourimage = \ 
          ListDensityPlot[Abs[Chop[Fourier[fourieraperature]]], 
            Mesh -> False, 
            DisplayLater]; 

Show[
          GraphicsArray[{
{theimage, 
                apimage}, 
{revfourimage, 
                fourmagimage}
}
], 
          ImageSize -> 1000, GraphicsSpacing -> { .001,  .0}, 
          DisplayNow];
]\)

The final image gets inverted through the center, I am not sure why this is happening...

ImageFourierAperature["~/3016/Mathematica/pentagon1.gif", 0, 1, 0, 
      1]

ImageFourierAperature["~/3016/Mathematica/pentagon2.gif", 0, 1, 0, 
      1]

ImageFourierAperature["~/3016/Mathematica/pentagon3.gif", 0, 1, 0, 
      1]

\(ImageFourierAperature["~/3016/Mathematica/pentagon1.gif",  .46,  .54,  .15,  .18]\)

\(ImageFourierAperature["~/3016/Mathematica/pentagon1.gif",  .3,  .325,  .1,  .125]\)

\(ImageFourierAperature["~/3016/Mathematica/polycrystal.gif",  .02,  .07,  .05,  .10]\)

ImageFourierAperature["~/3016/Mathematica/handler.gif", 0, 1, 0, 
      1]

The gray values of this image are stored in the (1,1) position of the image with gray values in the range (0,255)

\(ImageFourierAperature["~/3016/Mathematica/handler.gif",  .05,  .45,  .05,  .45]\)

\(ImageFourierAperature["~/3016/Mathematica/handler.gif",  .25,  .75,  .25,  .75]\)

ImageFourierAperature["~/3016/Mathematica/handler.gif",  .1,  .2, 
      0, 1]