kspace.FourierAnalysis#

class kspace.FourierAnalysis(width, ddims)[source]#

Bases: object

FFT-based analysis of scalar and vector fields on a regular grid: power spectra, divergence/curl, and vector-potential inversion.

Fourier transforms are always taken with the k = 2*pi*f convention and normalized by the cell volume dV (see fftn/ifftn) rather than scipy.fft’s default normalization.

Parameters:

  • width (array-like) – The physical size of the box [kpc] along each dimension.

  • ddims (array-like) – The number of grid cells along each dimension.

__init__(width, ddims)[source]#

Methods

__init__(width, ddims)

curl_of_field(data_vec[, periodic])

Compute the curl of a vector field using real-space finite differences.

divergence_component(data_vec[, diff_type, ...])

Project out the compressional (curl-free, longitudinal) component of a vector field: v_compressional_hat(k) = k_hat (k_hat .

divergence_of_field(data_vec[, periodic])

Compute the divergence of a vector field using real-space finite differences.

fftn(x, **kwargs)

Forward FFT of real-space data, normalized by the cell volume dV (F(k) = dV * sum_x f(x) exp(-i k.x)) rather than scipy.fft's default normalization.

generate_waves(diff_type)

Generate the wavenumber array used for spectral differentiation, for a choice of finite-difference approximation.

ifftn(x, **kwargs)

Inverse FFT of Fourier-space data, normalized by the cell volume dV (the exact inverse of fftn).

integrate_kspace(x[, axis])

Integrate a gridded quantity over k-space (e.g. to recover a total variance from a power spectrum), sum(x) * dVk / (2*pi)**ndim.

make_binned_powerspec(data, bins)

Compute the power spectrum of a field and bin it radially in the wavenumber magnitude into a 1D power spectrum.

make_powerspec(data)

Compute the gridded power spectrum P(k) = abs(data_hat(k))**2 / V of a scalar or vector field.

potential_of_field(data_vec[, diff_type, ...])

Invert B = curl(A) to recover the vector potential A for a divergence-free vector field B, in the Coulomb gauge (div(A) = 0): A_hat(k) = i(k x B_hat(k)) / (k .

solenoidal_component(data_vec[, diff_type, ...])

Project out the divergence-free (solenoidal) component of a vector field: v_solenoidal = v - v_compressional, where v_compressional is the component returned by divergence_component.

window_data(data[, filter_function])

https://stackoverflow.com/questions/27345861/extending-1d-function-across-3-dimensions-for-data-windowing

Attributes

khat

The unit wavenumber vector at each grid point, kvec / kmag, with the k=0 mode (division by zero) set to zero.

kmag

The magnitude of the wavenumber vector at each grid point, sqrt(kvec .

kvec

The wavenumber vector k at each grid point, of shape (ndim, *ddims), using the k = 2*pi*f convention.

kx

The x-component of the wavenumber vector at each grid point (kvec[0]).

ky

The y-component of the wavenumber vector at each grid point (kvec[1]).

kz

The z-component of the wavenumber vector at each grid point (kvec[2]).
curl_of_field(data_vec, periodic=False)[source]#

Compute the curl of a vector field using real-space finite differences.

Parameters:

  • data_vec (ndarray) – The vector field, of shape (ndim, *ddims).

  • periodic (boolean, optional) – If True, use periodic (wraparound) central differences instead of one-sided differences at the domain boundary. Only appropriate for data on a periodic domain (e.g. a field generated via FFT) – for general, non-periodic data the default (False) is correct. Default is False.

divergence_component(data_vec, diff_type='central', return_fft=False)[source]#

Project out the compressional (curl-free, longitudinal) component of a vector field: v_compressional_hat(k) = k_hat (k_hat . v_hat(k)).

Parameters:

  • data_vec (ndarray or FFTArray) – The vector field to project, of shape (ndim, *ddims).

  • diff_type (str, optional) – Which wavenumbers to use for the projection (passed to generate_waves). Default is “central”.

  • return_fft (boolean, optional) – If True, return the compressional component in Fourier space. Default is False.

divergence_of_field(data_vec, periodic=False)[source]#

Compute the divergence of a vector field using real-space finite differences.

Parameters:

  • data_vec (ndarray) – The vector field, of shape (ndim, *ddims).

  • periodic (boolean, optional) – If True, use periodic (wraparound) central differences instead of one-sided differences at the domain boundary. Only appropriate for data on a periodic domain (e.g. a field generated via FFT) – for general, non-periodic data the default (False) is correct. Default is False.

fftn(x, **kwargs)[source]#

Forward FFT of real-space data, normalized by the cell volume dV (F(k) = dV * sum_x f(x) exp(-i k.x)) rather than scipy.fft’s default normalization.

Parameters:

  • x (ndarray) – The data to transform, either a scalar field of shape ddims or a vector field of shape (ndim, *ddims).

  • **kwargs – Additional keyword arguments passed to scipy.fft.fftn.

Returns:

The transformed data.

Return type:

FFTArray

generate_waves(diff_type)[source]#

Generate the wavenumber array used for spectral differentiation, for a choice of finite-difference approximation.

Parameters:

diff_type (str) –

Which differencing scheme’s wavenumber to use:

  • ”continuum”: the exact FFT wavenumbers (kvec), i.e. an exact spectral derivative.

  • ”central”: the wavenumber corresponding to a periodic central-difference derivative, sin(k*dx)/dx.

  • ”forward”: the (complex-valued) wavenumber corresponding to a periodic forward-difference derivative.

Returns:

  • k (ndarray) – The wavenumber array, of shape (ndim, *ddims).

  • kmag (ndarray) – The Hermitian modulus sqrt(k . conj(k)) at each grid point.

ifftn(x, **kwargs)[source]#

Inverse FFT of Fourier-space data, normalized by the cell volume dV (the exact inverse of fftn).

Parameters:

  • x (FFTArray) – The Fourier-space data to transform, either a scalar field of shape ddims or a vector field of shape (ndim, *ddims). Must be an FFTArray (not a plain ndarray).

  • **kwargs – Additional keyword arguments passed to scipy.fft.ifftn.

Returns:

The real part of the transformed, real-space data.

Return type:

ndarray

integrate_kspace(x, axis=None)[source]#

Integrate a gridded quantity over k-space (e.g. to recover a total variance from a power spectrum), sum(x) * dVk / (2*pi)**ndim.

Parameters:

  • x (FFTArray) – The gridded quantity to integrate.

  • axis (int or tuple of ints, optional) – The axis or axes to integrate over. If None, integrates over all spatial axes. Default is None.

Returns:

The integrated quantity.

Return type:

ndarray or float

property khat#

The unit wavenumber vector at each grid point, kvec / kmag, with the k=0 mode (division by zero) set to zero.

Type:

ndarray

property kmag#

The magnitude of the wavenumber vector at each grid point, sqrt(kvec . kvec). Computed lazily on first access.

Type:

ndarray

property kvec#

The wavenumber vector k at each grid point, of shape (ndim, *ddims), using the k = 2*pi*f convention. Computed lazily on first access.

Type:

ndarray

property kx#

The x-component of the wavenumber vector at each grid point (kvec[0]).

Type:

ndarray

property ky#

The y-component of the wavenumber vector at each grid point (kvec[1]).

Type:

ndarray

property kz#

The z-component of the wavenumber vector at each grid point (kvec[2]).

Type:

ndarray

make_binned_powerspec(data, bins)[source]#

Compute the power spectrum of a field and bin it radially in the wavenumber magnitude into a 1D power spectrum.

Parameters:

  • data (ndarray or FFTArray) – The field, either real-space or already Fourier-transformed.

  • bins (int or ndarray) – If an int, the number of log-spaced bins to use between the smallest and largest wavenumbers resolved by the grid. If an ndarray, the bin edges to use directly.

Returns:

  • kbins (ndarray) – The bin edges used.

  • Pk (numpy.ma.MaskedArray) – The binned power spectrum, with empty bins masked.

make_powerspec(data)[source]#

Compute the gridded power spectrum P(k) = abs(data_hat(k))**2 / V of a scalar or vector field.

Parameters:

data (ndarray or FFTArray) – The field, either real-space or already Fourier-transformed.

Returns:

The power spectrum, same shape as the (Fourier-space) input.

Return type:

FFTArray

potential_of_field(data_vec, diff_type='continuum', return_fft=False)[source]#

Invert B = curl(A) to recover the vector potential A for a divergence-free vector field B, in the Coulomb gauge (div(A) = 0): A_hat(k) = i(k x B_hat(k)) / (k . k).

In 3D, A is returned as a 3-component vector field. In 2D, k and B_hat are both 2-vectors, so their cross product is a scalar and A is just the out-of-plane component.

Parameters:

  • data_vec (ndarray or FFTArray) – The divergence-free vector field B, of shape (ndim, *ddims).

  • diff_type (str, optional) – Which wavenumbers to use for the inversion (passed to generate_waves). Default is “continuum”, the exact FFT wavenumbers – this is the exact algebraic inverse of a curl computed the same way, not of curl_of_field, which uses real-space finite differences and will only agree with this at low k.

  • return_fft (boolean, optional) – If True, return the potential in Fourier space. Default is False.

solenoidal_component(data_vec, diff_type='central', return_fft=False)[source]#

Project out the divergence-free (solenoidal) component of a vector field: v_solenoidal = v - v_compressional, where v_compressional is the component returned by divergence_component.

Parameters:

  • data_vec (ndarray or FFTArray) – The vector field to project, of shape (ndim, *ddims).

  • diff_type (str, optional) – Which wavenumbers to use (passed to generate_waves). Default is “central”.

  • return_fft (boolean, optional) – If True, return the solenoidal component in Fourier space. Default is False.

window_data(data, filter_function='tukey', **kwargs)[source]#

https://stackoverflow.com/questions/27345861/extending-1d-function-across-3-dimensions-for-data-windowing

Performs an in-place windowing on N-dimensional spatial-domain data. This is done to mitigate boundary effects in the FFT.

Parameters:

  • data (ndarray) – Input data to be windowed, modified in place. Must be either the same shape as the FourierAnalysis instance, or be a 2-D or 3-D vector field with the same shape

  • filter_function (str) – Function can accept one argument: the window length, but all other keyword arguments are passed to the function. Default: tukey