pyoptics.fft¶

Helper module to deal with FFTs.

Convention: fields vary as exp(i(k.r - omega*t)). A field is a superposition of plane waves exp(+i*(k_x*x + k_y*y)); the forward transform FT() has the kernel exp(-i*(k_x*x + k_y*y)) and yields their amplitudes, the inverse transform inv_FT() sums them up again:

FT(f)[m, n] = sum_jl f[j, l]*exp(-i*(k_x[n]*x[l] + k_y[m]*y[j])),
inv_FT(F)[j, l] = 1/(Nx*Ny) sum_mn F[m, n]*exp(+i*(k_x[n]*x[l] + k_y[m]*y[j])).

Both take and return centered arrays: x = 0 at index N//2 in position space and k = 0 at index N//2 in frequency space, the latter matching utils.frequencies() and utils.freq_grid() with normal_order=True. For even N, utils.grid1d(-L/2, L/2, N) gives such a centered x. For odd N use e.g. dx*(np.arange(N) - N//2). (If x is not centered, FT() equals the sum above up to the constant phase factor exp(+i*k*x[N//2]), which inv_FT() removes again.)

FT() and inv_FT() are unnormalized like numpy.fft.fft2() and ifft2(); multiply FT() by dx*dy to approximate the continuous Fourier integral. FT_unitary() and inv_FT_unitary() are scaled by 1/sqrt(Nx*Ny) and sqrt(Nx*Ny) such that FT_unitary() conserves the sum of |f|**2.

rFT() and inv_rFT() are the versions for real input. Their spectra only hold the non-negative k_x (the columns N//2, N//2 + 1, … of the FT() spectrum, for even N followed by the Nyquist frequency); only the k_y axis is centered.

czt() is Bluestein’s chirp-z transform. It evaluates sums of the above kind for arbitrary, equidistant output frequencies (zoomed or shifted spectra) with FFTs, in O((N + M)*log(N + M)) operations. zoom_FT() and zoom_inv_FT() use it to evaluate FT() and inv_FT() on any equidistant grid of output frequencies or positions, e.g. a fine grid over a small part of a spectrum.

The FFTs are scipy.fft’s (same sign convention and normalization as numpy.fft), run with the number of threads in the module variable workers (default -1: all CPU cores). On multi-core machines they are 2 to 4 times faster than numpy.fft for grids of 512**2 and more. Set pyoptics.fft.workers = 1 for single threaded transforms, e.g. when running several simulations in parallel processes.

References¶

L. R. Rabiner, R. W. Schafer and C. M. Rader, “The chirp z-transform algorithm”, IEEE Trans. Audio Electroacoust. 17, 86-92 (1969).

L. I. Bluestein, “A linear filtering approach to the computation of discrete Fourier transform”, IEEE Trans. Audio Electroacoust. 18, 451-455 (1970).

Functions

`FT`
`FT_unitary`
`czt`	Chirp-z transform (Bluestein's algorithm).
`fft`
`fft2`
`ifft`
`ifft2`
`inv_FT`
`inv_FT_unitary`
`inv_rFT`	Inverse of rFT().
`inv_rFT_unitary`	Unitary inverse of rFT() (scaled by the square root of the array size).
`irfft2`
`rFT`	FT() for real x.
`rFT_unitary`
`rfft2`
`zoom_FT`	FT() evaluated on an arbitrary equidistant frequency grid.
`zoom_inv_FT`	inv_FT() evaluated on an arbitrary equidistant grid of positions.

pyoptics.fft.rFT(x)[source]¶: FT() for real x. Returns the spectrum for k_x >= 0 only (shape (Ny, Nx//2 + 1)); k_y is centered.

pyoptics.fft.inv_rFT(X, shape=None)[source]¶

Inverse of rFT().

Parameters:

X2d array: Half spectrum as returned by rFT().
shapetuple, optional: Shape of the real array; needs to be given for an odd number of x samples.

Returns:

x2d array: Real array.

pyoptics.fft.inv_rFT_unitary(X, shape=None)[source]¶

Unitary inverse of rFT() (scaled by the square root of the array size).

Parameters:

X2d array: Half spectrum.
shapetuple, optional: Shape of the real array, see inv_rFT().

Returns:

x2d array: Real array.

pyoptics.fft.czt(x, M=None, w=None, a=1.0, axis=-1)[source]¶

Chirp-z transform (Bluestein’s algorithm).

Computes

X[m] = sum_{n=0}^{N-1} x[n]*a**(-n)*w**(n*m),   m = 0, ..., M-1,

i.e. the z transform of x on the points z_m = a*w**(-m), with FFTs of length >= N + M - 1. With a = 1, w = exp(-2j*pi/N) and M = N, this is the DFT numpy.fft.fft(x). For a = exp(1j*theta_0) and w = exp(-1j*dtheta), the points z_m lie on the unit circle at the angles theta_0 + m*dtheta (a zoomed DFT); this is the case used in pyoptics. (scipy.signal.czt() uses the same definition.)

Parameters:

xarray: Input data; the transform runs along axis.
Mint, optional: Number of output points. Defaults to the length N of x along axis.
wcomplex, optional: Ratio between subsequent points. Defaults to exp(-2j*pi/M).
acomplex, optional: Starting point. Defaults to 1.
axisint, optional: Axis to transform.

Returns:

Xarray: Complex array with M points along axis.

Notes

Bluestein’s identity n*m = (n**2 + m**2 - (m - n)**2)/2 turns the sum into a convolution with the chirp w**(-n**2/2), evaluated cyclically by FFT. The chirps are computed as exp(log(w)*n**2/2) from a single logarithm, which keeps the product of the three chirps equal to w**(n*m) for any branch of the logarithm. For |w| != 1 the chirps grow or decay exponentially; this is accurate only for moderate N and M.

References

L. Bluestein, “A linear filtering approach to the computation of discrete Fourier transform,” IEEE Trans. Audio Electroacoust. 18, 451-455 (1970) L. Rabiner, R. Schafer and C. Rader, “The chirp z-transform algorithm,” IEEE Trans. Audio Electroacoust. 17, 86-92 (1969)

pyoptics.fft.zoom_FT(f, x, y, kx_out, ky_out)[source]¶

FT() evaluated on an arbitrary equidistant frequency grid.

Computes

F[m, n] = sum_jl f[j, l]*exp(-i*(kx_out[n]*x[l] + ky_out[m]*y[j]))

with two chirp-z transforms. The frequencies may be finer or coarser than 2*pi/L, shifted and of any number; the sum is periodic in k_x with the period 2*pi/dx. For centered x, y (see the module docstring) and the frequencies kx_out, ky_out of utils.frequencies(), F equals FT(f). Unlike FT(), the sum always uses the coordinates as given: no grid centering is needed and the phase of the result refers to the origin x = y = 0.

Parameters:

farray: Field of shape (Ny, Nx).
x, yarrays: Linear, equidistant coordinates of f.
kx_out, ky_outarrays: Linear, equidistant output wavenumbers (angular spatial frequencies).

Returns:

Farray: Spectrum of shape (len(ky_out), len(kx_out)). Multiply by dx*dy to approximate the Fourier integral.

pyoptics.fft.zoom_inv_FT(F, kx, ky, x_out, y_out)[source]¶

inv_FT() evaluated on an arbitrary equidistant grid of positions.

Computes the plane wave sum

f[j, l] = 1/(Nx*Ny) sum_mn F[m, n]*exp(+i*(kx[n]*x_out[l] + ky[m]*y_out[j]))

with two chirp-z transforms, i.e. the trigonometric interpolation of inv_FT(F) on any equidistant grid (the result has the period 2*pi/dk_x in x). For x_out, y_out equal to the centered grid that belongs to kx, ky, f equals inv_FT(F).

Parameters:

Farray: Plane wave amplitudes of shape (Ny, Nx), e.g. from FT().
kx, kyarrays: Linear, equidistant wavenumbers of F, e.g. utils.frequencies(x).
x_out, y_outarrays: Linear, equidistant output coordinates.

Returns:

farray: Field of shape (len(y_out), len(x_out)).