Goals

  • Determine how Julia's FFTW implementation works for some simple 1D and 2D Gaussian test cases, including units.
  • [Separate notebook] Read in RADMC3D image and try 2D FFT

Conclusions

It seems as though the fftshift behaviour is the same as that required in python.


In [45]:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
using constants
using visibilities
import PyPlot.plt

In [66]:
# Create an x-axis
#dx = 0.1 
#xx = Float64[i*dx for i=-10:9];
n = 64
xx = fftspace(5, n)
dx = xx[2] - xx[1];
ss = fftshift(fftfreq(n, dx));

In [48]:
function Gauss(x::Vector{Float64})
    exp(-pi * x.^2)
end

function GaussFT(s::Vector{Float64})
    exp(-pi * s.^2)
end


Out[48]:
GaussFT (generic function with 1 method)

In [67]:
yy = Gauss(xx);

In [68]:
plt.plot(xx,yy, "o")


Out[68]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ffd86ae4550>

In [69]:
YY = dx * fftshift(fft(fftshift(yy)))


Out[69]:
64-element Array{Complex{Float64},1}:
          2.13024e-14+0.0im
 7.86871e-14+1.51788e-17im 
 5.25691e-13+8.67362e-18im 
  3.35426e-12+2.1684e-18im 
 2.01043e-11+8.67362e-18im 
  1.13161e-10+1.0842e-17im 
 5.98158e-10-8.67362e-18im 
  2.96926e-9-3.46945e-17im 
           1.38418e-8+0.0im
  6.05967e-8-1.83995e-17im 
           2.49126e-7+0.0im
  9.61834e-7+2.02597e-17im 
  3.48734e-6-8.67362e-18im 
                     ⋮     
  3.48734e-6-8.67362e-18im 
  9.61834e-7-1.51788e-17im 
  2.49126e-7+1.73472e-17im 
  6.05967e-8+2.60209e-17im 
           1.38418e-8+0.0im
  2.96926e-9+3.06659e-17im 
          5.98158e-10+0.0im
 1.13161e-10-2.45965e-17im 
 2.01043e-11+8.67362e-18im 
 3.35423e-12-1.80275e-17im 
 5.25621e-13-3.39504e-17im 
  7.8583e-14+2.84865e-18im 

In [70]:
plt.plot(ss, YY, "o")
plt.plot(ss, GaussFT(ss))


Out[70]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ffd86a3ea90>

In [71]:
plt.plot(ss, YY - GaussFT(ss), "o")


Out[71]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ffd869af7f0>

2D model

Because the disk is not symmetric on the sky (position angle warp), the output will have imaginary components. But, because the actual image only contains real values (sky intensities) then we can still use rfft.


In [72]:
# 2D gaussian
function Gauss(x::Float64, y::Float64)
    exp(-pi * (x.^2 + y.^2))
end


Out[72]:
Gauss (generic function with 3 methods)

In [76]:
# 2D Gaussian image vectorized over input arrays
function Gauss(x::Vector{Float64}, y::Vector{Float64})
    
    nx = length(x)
    ny = length(y)
    
    out = Array(Float64, (ny, nx))

    for j=1:ny
        for i=1:nx
            out[j,i] = Gauss(x[i], y[j])
        end
    end
    
    return out
end

function GaussFT(u::Vector{Float64}, v::Vector{Float64})
    
    nx = length(u)
    ny = length(v)
    
    out = Array(Float64, (ny, nx))

    for j=1:ny
        for i=1:nx
            out[j,i] = Gauss(u[i], v[j])
        end
    end
    
    return out
end


Out[76]:
GaussFT (generic function with 2 methods)

In [77]:
img = Gauss(xx, xx)


Out[77]:
64x64 Array{Float64,2}:
 6.04202e-69  7.58072e-67  8.15865e-65  …  8.15865e-65  7.58072e-67
 7.58072e-67  9.51127e-65  1.02364e-62     1.02364e-62  9.51127e-65
 8.15865e-65  1.02364e-62  1.10168e-60     1.10168e-60  1.02364e-62
 7.53194e-63  9.45007e-61  1.01705e-58     1.01705e-58  9.45007e-61
 5.96451e-61  7.48347e-59  8.05399e-57     8.05399e-57  7.48347e-59
 4.05157e-59  5.08337e-57  5.47091e-55  …  5.47091e-55  5.08337e-57
 2.36076e-57  2.96196e-55  3.18778e-53     3.18778e-53  2.96196e-55
 1.17994e-55  1.48043e-53  1.5933e-51      1.5933e-51   1.48043e-53
 5.05881e-54  6.34712e-52  6.83101e-50     6.83101e-50  6.34712e-52
 1.86045e-52  2.33424e-50  2.5122e-48      2.5122e-48   2.33424e-50
 5.86902e-51  7.36366e-49  7.92505e-47  …  7.92505e-47  7.36366e-49
 1.58816e-49  1.99261e-47  2.14452e-45     2.14452e-45  1.99261e-47
 3.68641e-48  4.62521e-46  4.97782e-44     4.97782e-44  4.62521e-46
 ⋮                                      ⋱                          
 3.68641e-48  4.62521e-46  4.97782e-44     4.97782e-44  4.62521e-46
 1.58816e-49  1.99261e-47  2.14452e-45     2.14452e-45  1.99261e-47
 5.86902e-51  7.36366e-49  7.92505e-47     7.92505e-47  7.36366e-49
 1.86045e-52  2.33424e-50  2.5122e-48   …  2.5122e-48   2.33424e-50
 5.05881e-54  6.34712e-52  6.83101e-50     6.83101e-50  6.34712e-52
 1.17994e-55  1.48043e-53  1.5933e-51      1.5933e-51   1.48043e-53
 2.36076e-57  2.96196e-55  3.18778e-53     3.18778e-53  2.96196e-55
 4.05157e-59  5.08337e-57  5.47091e-55     5.47091e-55  5.08337e-57
 5.96451e-61  7.48347e-59  8.05399e-57  …  8.05399e-57  7.48347e-59
 7.53194e-63  9.45007e-61  1.01705e-58     1.01705e-58  9.45007e-61
 8.15865e-65  1.02364e-62  1.10168e-60     1.10168e-60  1.02364e-62
 7.58072e-67  9.51127e-65  1.02364e-62     1.02364e-62  9.51127e-65

In [75]:
plt.imshow(img, interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))


Out[75]:
PyObject <matplotlib.image.AxesImage object at 0x7ffd869988d0>

In [78]:
plt.imshow(fftshift(img), interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))


Out[78]:
PyObject <matplotlib.image.AxesImage object at 0x7ffd868fdf98>

In [89]:
imout = dx * dx * fftshift(fft(fftshift(img)))


Out[89]:
64x64 Array{Complex{Float64},2}:
 -1.18146e-17+0.0im  -1.25118e-17+0.0im  …  -1.25118e-17+0.0im
  2.26856e-17+0.0im   2.19832e-17+0.0im      2.19832e-17+0.0im
 -7.63232e-18+0.0im  -8.51314e-18+0.0im     -8.51314e-18+0.0im
 -4.17086e-17+0.0im  -4.14881e-17+0.0im     -4.14881e-17+0.0im
 -1.11336e-17+0.0im  -1.14391e-17+0.0im     -1.14391e-17+0.0im
  -5.6599e-18+0.0im  -5.78659e-18+0.0im  …  -5.78659e-18+0.0im
 -7.09361e-18+0.0im  -7.04296e-18+0.0im     -7.04296e-18+0.0im
 -1.24138e-17+0.0im   -1.2317e-17+0.0im      -1.2317e-17+0.0im
 -2.04189e-17+0.0im   -1.9944e-17+0.0im      -1.9944e-17+0.0im
 -1.69986e-17+0.0im  -1.68548e-17+0.0im     -1.68548e-17+0.0im
 -2.29704e-18+0.0im    -2.143e-18+0.0im  …    -2.143e-18+0.0im
 -3.89974e-18+0.0im  -3.23982e-18+0.0im     -3.23982e-18+0.0im
 -7.72421e-18+0.0im  -7.16338e-18+0.0im     -7.16338e-18+0.0im
             ⋮                           ⋱                    
 -7.72421e-18+0.0im  -7.16338e-18+0.0im     -7.16338e-18+0.0im
 -3.89974e-18+0.0im  -3.23982e-18+0.0im     -3.23982e-18+0.0im
 -2.29704e-18+0.0im    -2.143e-18+0.0im       -2.143e-18+0.0im
 -1.69986e-17+0.0im  -1.68548e-17+0.0im  …  -1.68548e-17+0.0im
 -2.04189e-17+0.0im   -1.9944e-17+0.0im      -1.9944e-17+0.0im
 -1.24138e-17+0.0im   -1.2317e-17+0.0im      -1.2317e-17+0.0im
 -7.09361e-18+0.0im  -7.04296e-18+0.0im     -7.04296e-18+0.0im
  -5.6599e-18+0.0im  -5.78659e-18+0.0im     -5.78659e-18+0.0im
 -1.11336e-17+0.0im  -1.14391e-17+0.0im  …  -1.14391e-17+0.0im
 -4.17086e-17+0.0im  -4.14881e-17+0.0im     -4.14881e-17+0.0im
 -7.63232e-18+0.0im  -8.51314e-18+0.0im     -8.51314e-18+0.0im
  2.26856e-17+0.0im   2.19832e-17+0.0im      2.19832e-17+0.0im

In [92]:
plt.imshow(real(imout), interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))
plt.colorbar()


Out[92]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7ffd865a35c0>

In [93]:
plt.imshow(imag(imout), interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))
plt.colorbar()


Out[93]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7ffd86377c88>

In [94]:
imoutFT = GaussFT(ss, ss)


Out[94]:
64x64 Array{Float64,2}:
 1.14175e-28  8.26284e-28  5.61564e-27  …  5.61564e-27  8.26284e-28
 8.26284e-28  5.9798e-27   4.06402e-26     4.06402e-26  5.9798e-27 
 5.61564e-27  4.06402e-26  2.76201e-25     2.76201e-25  4.06402e-26
 3.58411e-26  2.59381e-25  1.76282e-24     1.76282e-24  2.59381e-25
 2.1482e-25   1.55465e-24  1.05658e-23     1.05658e-23  1.55465e-24
 1.20916e-24  8.75064e-24  5.94715e-23  …  5.94715e-23  8.75064e-24
 6.39149e-24  4.6255e-23   3.14361e-22     3.14361e-22  4.6255e-23 
 3.17274e-23  2.2961e-22   1.56049e-21     1.56049e-21  2.2961e-22 
 1.47903e-22  1.07037e-21  7.27453e-21     7.27453e-21  1.07037e-21
 6.47493e-22  4.68589e-21  3.18465e-20     3.18465e-20  4.68589e-21
 2.66198e-21  1.92647e-20  1.30928e-19  …  1.30928e-19  1.92647e-20
 1.02775e-20  7.43778e-20  5.05491e-19     5.05491e-19  7.43778e-20
 3.72632e-20  2.69673e-19  1.83277e-18     1.83277e-18  2.69673e-19
 ⋮                                      ⋱                          
 3.72632e-20  2.69673e-19  1.83277e-18     1.83277e-18  2.69673e-19
 1.02775e-20  7.43778e-20  5.05491e-19     5.05491e-19  7.43778e-20
 2.66198e-21  1.92647e-20  1.30928e-19     1.30928e-19  1.92647e-20
 6.47493e-22  4.68589e-21  3.18465e-20  …  3.18465e-20  4.68589e-21
 1.47903e-22  1.07037e-21  7.27453e-21     7.27453e-21  1.07037e-21
 3.17274e-23  2.2961e-22   1.56049e-21     1.56049e-21  2.2961e-22 
 6.39149e-24  4.6255e-23   3.14361e-22     3.14361e-22  4.6255e-23 
 1.20916e-24  8.75064e-24  5.94715e-23     5.94715e-23  8.75064e-24
 2.1482e-25   1.55465e-24  1.05658e-23  …  1.05658e-23  1.55465e-24
 3.58411e-26  2.59381e-25  1.76282e-24     1.76282e-24  2.59381e-25
 5.61564e-27  4.06402e-26  2.76201e-25     2.76201e-25  4.06402e-26
 8.26284e-28  5.9798e-27   4.06402e-26     4.06402e-26  5.9798e-27 

In [95]:
plt.imshow(imoutFT, interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))
plt.colorbar()


Out[95]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7ffd862cee48>

In [96]:
plt.imshow(real(imout) - imoutFT, interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"))
plt.colorbar()


Out[96]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7ffd863c9438>

In [ ]: