pypret.fourier module

This module implements the Fourier transforms on linear grids.

The following code approximates the continuous Fourier transform (FT) on equidistantly spaced grids. While this is usually associated with ‘just doing a fast Fourier transform (FFT)’, surprisingly, much can be done wrong.

The reason is that the correct expressions depend on the grid location. In fact, the FT can be calculated with one DFT but in general it requires a prior and posterior multiplication with phase factors.

The FT convention we are going to use here is the following:

Ẽ(w) = 1/2pi ∫ E(t) exp(+i w t) dt
E(t) =       ∫ Ẽ(w) exp(-i t w) dw

where w is the angular frequency. We can approximate these integrals by their Riemann sums on the following equidistantly spaced grids:

t_k = t_0 + k Δt, k=0, ..., N-1
w_n = w_0 + n Δw, n=0, ..., N-1

and define E_k = E(t_k) and Ẽ_n = Ẽ(w_n) to obtain:

Ẽ_n = Δt/2pi ∑_k E_k exp(+i w_n t_k)
E_k = Δw     ∑_n Ẽ_n exp(-i t_k w_n).

To evaluate the sum using the FFT we can expand the exponential to obtain:

Ẽ_n = Δt/2pi exp(+i n t_0 Δw) ∑_k [E_k exp(+i t_k w_0) ] exp(+i n k Δt Δw)
E_k = Δw     exp(-i t_k w_0)  ∑_n [Ẽ_n exp(-i n t_0 Δw)] exp(-i k n Δt Δw)

Additionally, we have to require the so-called reciprocity relation for the grid spacings:

      !
Δt Δw = 2pi / N = ζ     (reciprocity relation)

This is what enables us to use the DFT/FFT! Now we look at the definition of the FFT in NumPy:

 fft[x_m] -> X_k =     ∑_m exp(-2pi i m k / N)
ifft[X_k] -> x_m = 1/N ∑_k exp(+2pi i k m / N)

which gives the final expressions:

Ẽ_n = Δt N/2pi r_n   ifft[E_k s_k  ]
E_k = Δw       s_k^*  fft[Ẽ_n r_n^*]

with r_n = exp(+i n t_0 Δw)
     s_k = exp(+i t_k w_0)

where ^* means complex conjugation. We see that the array to be transformed has to be multiplied with an appropriate phase factor before and after performing the DFT. And those phase factors mainly depend on the starting points of the grids: w_0 and t_0. Note also that due to our sign convention for the FT we have to use ifft for the forward transform and vice versa.

Trivially, we can see that for w_0 = t_0 = 0 the phase factors vanish and the FT is approximated well by just the DFT. However, in optics these grids are unusual. For w_0 = l Δw and t_0 = m Δt, where l, m are integers (i.e., w_0 and t_0 are multiples of the grid spacing), the phase factors can be incorperated into the DFT. Then the phase factors can be replaced by circular shifts of the input and output arrays.

This is exactly what the functions (i)fftshift are doing for one specific choice of l and m, namely for:

t_0 = -floor(N/2) Δt
w_0 = -floor(N/2) Δw.

In this specific case only we can approximate the FT by:

Ẽ_n = Δt N/2pi fftshift(ifft(ifftshift(E_k)))
E_k = Δw       fftshift( fft(ifftshift(Ẽ_n))) (no mistake!)

We see that the ifftshift _always_ has to appear on the inside. Failure to do so will still be correct for even N (here fftshift is the same as ifftshift) but will produce wrong results for odd N.

Additionally you have to watch out not to violate the assumptions for the grid positions. Using a symmetrical grid, e.g.,:

x = linspace(-1, 1, 128)

will also produce wrong results, as the elements of x are not multiples of the grid spacing (but shifted by half a grid point).

The main drawback of this approach is that circular shifts are usually far more time- and memory-consuming than an elementwise multiplication, especially for higher dimensions. In fact I see no advantage in using the shift approach at all. But for some reason it got stuck in the minds of people and you find the notion of having to re-order the output of the DFT everywhere.

Long story short: here we are going to stick with multiplying the correct phase factors. The code tries to follow the notation used above.

Good, more comprehensive expositions of the issues above can be found in [Briggs1995] and [Hansen2014]. For the reason why the first-order approximation to the Riemann integral suffices, see [Trefethen2014].

class pypret.fourier.FourierTransform(N, dt=None, dw=None, t0=None, w0=None)[source]
backward(x, out=None)[source]

Calculates the backward (inverse) Fourier transform of x.

For n-dimensional arrays it operates on the last axis, which has to match the size of x.

Parameters:
  • x (ndarray) – The array of which the Fourier transform will be calculated.
  • out (ndarray or None, optional) – A location into which the result is stored. If not provided or None, a freshly-allocated array is returned.
forward(x, out=None)[source]

Calculates the (forward) Fourier transform of x.

For n-dimensional arrays it operates on the last axis, which has to match the size of x.

Parameters:
  • x (ndarray) – The array of which the Fourier transform will be calculated.
  • out (ndarray or None, optional) – A location into which the result is stored. If not provided or None, a freshly-allocated array is returned.