Using xrft#

xrft is a package to do Fourier transforms of xarray data easily. One can use the numpy-like functions, which are not exact (phase and amplitude shifts), but deliver the same results as a numpy-FT would do. There is also the possibility to correct these shifts with xrft. Another handy advantage of xrft is that the visualization of the results can be done with HoloViews in a very simple manner.

(Note the generation of each plot takes about 5 secs., you can add the dynamic=True parameter, if you want to load the plots faster. This is thoroughly explained in the holoviews tutorial)

[1]:

import xrft

%config InlineBackend.figure_formats = ['svg']


import holoviews as hv

from postopus.octopus_run import Run

hv.extension("bokeh", "matplotlib")  # Allow for interactive plots

input file is already defined in the folder (s. GitLab repo), otherwise we recommend defining it in the notebook

[2]:

cd ../octopus_data/interference/

/builds/octopus-code/postopus/docs/octopus_data/interference

Assuming you have octopus in your PATH:

[3]:

!octopus > out_gs.log 2>&1

[4]:

run = Run(".")

Original Data#

[5]:

xa = run.Maxwell.td.b_field(source="z=0").vx

[6]:

xa

[6]:

<xarray.DataArray 'vx' (t: 17, x: 41, y: 41)> Size: 229kB
[28577 values with dtype=float64]
Coordinates:
    step     (t) int64 136B 0 10 20 30 40 50 60 ... 100 110 120 130 140 150 160
  * t        (t) float64 136B 0.0 0.02107 0.04213 0.0632 ... 0.2949 0.316 0.337
  * x        (x) float64 328B -10.0 -9.5 -9.0 -8.5 -8.0 ... 8.0 8.5 9.0 9.5 10.0
  * y        (y) float64 328B -10.0 -9.5 -9.0 -8.5 -8.0 ... 8.0 8.5 9.0 9.5 10.0
Attributes:
    units:    au

Visualization with HoloViews#

[7]:

hv_dst_orig = hv.Dataset(xa)
hv_imt_orig = hv_dst_orig.to(hv.Image, ["x", "y"])
hv.output(max_frames=3000)
hv_imt_orig.opts(
    colorbar=True, width=500, height=400, cmap="seismic", clim=(-(10**-3), 10**-3)
)

[7]:

Numpy-like fft#

[8]:

# `xa.drop_vars("step")` is needed here as "step" and "t" represent the same dimension.
# Fxa = xrft.fft(xa.drop_vars("step"), dim="t") # numpy.fft-like behaviour
Fxa = xrft.fft(xa.drop_vars("step"), dim="t", real_dim="t")  # numpy.rfft-like behaviour

[9]:

Fxa

[9]:

<xarray.DataArray (freq_t: 9, x: 41, y: 41)> Size: 242kB
array([[[-1.19793340e-05+0.00000000e+00j,
         -1.51782559e-05+0.00000000e+00j,
         -1.99394843e-05+0.00000000e+00j, ...,
         -1.86336875e-05+0.00000000e+00j,
         -1.37752569e-05+0.00000000e+00j,
         -1.14321021e-05+0.00000000e+00j],
        [-1.06707295e-05+0.00000000e+00j,
         -9.83375078e-06+0.00000000e+00j,
         -1.21946988e-05+0.00000000e+00j, ...,
         -2.41516278e-05+0.00000000e+00j,
         -1.72796354e-05+0.00000000e+00j,
         -1.13740378e-05+0.00000000e+00j],
        [-1.14535267e-05+0.00000000e+00j,
         -8.13124846e-06+0.00000000e+00j,
         -8.11934843e-06+0.00000000e+00j, ...,
         -2.33994956e-05+0.00000000e+00j,
         -2.10661716e-05+0.00000000e+00j,
         -1.49552563e-05+0.00000000e+00j],
        ...,
        [-2.46157413e-05+0.00000000e+00j,
...
         -7.41306307e-05-3.31533773e-05j],
        ...,
        [ 2.66486388e-05-2.92516002e-05j,
          2.96726524e-06+9.66818951e-06j,
         -1.81344937e-05+4.73515077e-05j, ...,
          2.98016204e-05+1.18512323e-04j,
          2.31442888e-05+1.00987975e-04j,
          1.23292186e-05+7.05357745e-05j],
        [ 5.04332157e-05-6.16614708e-05j,
          2.73593815e-05-2.77036632e-05j,
          2.62716076e-06+1.28677947e-05j, ...,
          2.59060830e-05+1.20600298e-04j,
          2.41800082e-05+1.18674163e-04j,
          1.71752975e-05+9.94970001e-05j],
        [ 6.78150243e-05-8.26997336e-05j,
          4.92095458e-05-5.69909680e-05j,
          2.52443767e-05-2.23129418e-05j, ...,
          1.71446115e-05+1.03976907e-04j,
          1.95656689e-05+1.17210937e-04j,
          1.77978986e-05+1.13644828e-04j]]], shape=(9, 41, 41))
Coordinates:
  * x        (x) float64 328B -10.0 -9.5 -9.0 -8.5 -8.0 ... 8.0 8.5 9.0 9.5 10.0
  * y        (y) float64 328B -10.0 -9.5 -9.0 -8.5 -8.0 ... 8.0 8.5 9.0 9.5 10.0
  * freq_t   (freq_t) float64 72B 0.0 2.792 5.585 8.377 ... 16.75 19.55 22.34

Visualization with HoloViews#

[10]:

Fxa.name = "testFT"  # needed to use holoviews
hv_dst = hv.Dataset(Fxa.real)
hv_imt = hv_dst.to(hv.Image, ["x", "y"])

[11]:

hv_imt.opts(
    colorbar=True, width=500, height=400, cmap="seismic", clim=(-(10**-3), 10**-3)
)

[11]:

Inverse transform#

There is a shift when inverse transforming.

[12]:

# iFxa = xrft.ifft(Fxa, dim="freq_t") # Inverse of np.fft
iFxa = xrft.ifft(Fxa, dim="freq_t", real_dim="freq_t")  # Inverse of np.rfft

Visualization with HoloViews#

[13]:

iFxa.name = "iFtest"
hv_dst_ift = hv.Dataset(iFxa.real)
hv_imt_ift = hv_dst_ift.to(hv.Image, ["x", "y"])
hv.output(max_frames=3000)
hv_imt_ift.opts(
    colorbar=True, width=500, height=400, cmap="seismic", clim=(-(10**-3), 10**-3)
)

[13]:

FFT with true amplitude and true phase#

This is good for the cases in which the data is not centered at 0, see https://xrft.readthedocs.io/en/latest/DFT-iDFT_example.html. As one can see we go back to the original data after the inverse transform. There is no shift.

[14]:

Fxa2 = xrft.fft(xa.drop_vars("step"), dim="t", true_amplitude=True, true_phase=True)

Visualization with HoloViews#

[15]:

Fxa2.name = "testFT2"
hv_dst2 = hv.Dataset(Fxa2.real)
hv_imt_2 = hv_dst2.to(hv.Image, ["x", "y"])
hv_imt_2.opts(
    colorbar=True, width=500, height=400, cmap="seismic", clim=(-(10**-3), 10**-3)
)

[15]:

Inverse Transform#

[16]:

iFxa2 = xrft.ifft(Fxa2, dim="freq_t", true_amplitude=True, true_phase=True)

Visualization with HoloViews#

[17]:

iFxa2.name = "iFtest2"
hv_dst_ift2 = hv.Dataset(iFxa2.real)
hv_imt_ift2 = hv_dst_ift2.to(hv.Image, ["x", "y"])
hv_imt_ift2.opts(
    colorbar=True, width=500, height=400, cmap="seismic", clim=(-(10**-3), 10**-3)
)

[17]:

Spectral density#

[18]:

sd = xrft.power_spectrum(xa.drop_vars("step"), dim="t", real_dim="t", scaling="density")

Visualization with HoloViews#

[19]:

sd.name = "Spectral_density"
hv_dst_sd = hv.Dataset(sd)
hv_imt_sd = hv_dst_sd.to(hv.Image, ["x", "y"])
hv_imt_sd.opts(colorbar=True, width=500, height=400, cmap="seismic")

[19]: