FFT

ducc0.fft

Fast Fourier, sine/cosine, and Hartley transforms.

This module supports
  • single, double, and long double precision

  • complex and real-valued transforms

  • multi-dimensional transforms

For two- and higher-dimensional transforms the code will use SSE2 and AVX vector instructions for faster execution if these are supported by the CPU and were enabled during compilation.

ducc0.fft.c2c(a: numpy.ndarray, axes: object = None, forward: bool = True, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a complex FFT.

Parameters:
  • a (numpy.ndarray (any complex or real type)) – The input data. If its type is real, a more efficient real-to-complex transform will be used.

  • axes (list of integers) – The axes along which the FFT is carried out. If not set, all axes will be transformed.

  • forward (bool) – If True, a negative sign is used in the exponent, else a positive one.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed axes.

  • out (numpy.ndarray (same shape as a, complex type with same accuracy as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data.

Return type:

numpy.ndarray (same shape as a, complex type with same accuracy as a)

ducc0.fft.c2r(a: numpy.ndarray, axes: object = None, lastsize: int = 0, forward: bool = True, inorm: int = 0, out: object = None, nthreads: int = 1, allow_overwriting_input: bool = False) numpy.ndarray

Performs an FFT whose output is strictly real.

Parameters:
  • a (numpy.ndarray (any complex type)) – The input data

  • axes (list of integers) – The axes along which the FFT is carried out. If not set, this is assumed to be list(range(a.ndim)). The complex-to-real transform will be executed along axes[-1], and will be executed last.

  • lastsize (the output size of the last axis to be transformed.) – If the corresponding input axis has size n, this can be 2*n-2 or 2*n-1.

  • forward (bool) – If True, a negative sign is used in the exponent, else a positive one.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed output axes.

  • out (numpy.ndarray (real type with same accuracy as a)) – For the required shape, see the Returns section. Must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

  • allow_overwriting_input (bool) – If True, the input array may be overwritten with some intermediate data. This can avoid allocating temporary variables and improve performance.

Returns:

The transformed data. The shape is identical to that of the input array, except for axes[-1], which has now lastsize entries.

Return type:

numpy.ndarray (real type with same accuracy as a)

ducc0.fft.convolve_axis(in: numpy.ndarray, out: numpy.ndarray, axis: int, kernel: numpy.ndarray, nthreads: int = 1) numpy.ndarray

Performs a circular convolution along one axis.

The result is equivalent to

import scipy.ndimage
import scipy.signal
import scipy.fft
kernel = scipy.fft.fftshift(kernel)
tmp = scipy.ndimage.convolve1d(in, kernel, axis, mode='wrap')
out[()] = scipy.signal.resample(tmp, out.shape[axis], axis=axis)
return out
Parameters:
  • in (numpy.ndarray (any real or complex type)) – The input data

  • out (numpy.ndarray (same type as in)) – The output data. Must have the same shape as in except for the axis to be convolved

  • axis (integer) – The axis along which the convolution is carried out.

  • kernel (one-dimensional numpy.ndarray (same type as in)) – The kernel to be used for convolution The length of this array must be equal to in.shape[axis]

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The convolved input

Return type:

numpy.ndarray (identical to out)

Notes

The main purpose of this routine is efficiency: the combination of the above operations can be carried out more quickly than running the individual operations in succession.

If in.shape[axis]!=out.shape[axis], the appropriate amount of zero-padding or truncation will be carried out after the convolution step.

in and out may overlap in memory. If they do, their first elements must be at the same memory location, and all their strides must be equal.

ducc0.fft.dct(a: numpy.ndarray, type: int, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a discrete cosine transform.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • type (integer) – the type of DCT. Must be in [1; 4].

  • axes (list of integers) – The axes along which the transform is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • inorm (integer) –

    the normalization type
    0 : no normalization
    1 : make transform orthogonal and divide by sqrt(N)
    2 : divide by N

    where N is the product of n_i for every transformed axis i. n_i is 2*(<axis_length>-1 for type 1 and 2*<axis length> for types 2, 3, 4. Making the transform orthogonal involves the following additional steps for every 1D sub-transform:

    Type 1

    multiply first and last input value by sqrt(2); divide first and last output value by sqrt(2)

    Type 2

    divide first output value by sqrt(2)

    Type 3

    multiply first input value by sqrt(2)

    Type 4

    nothing

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.dst(a: numpy.ndarray, type: int, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a discrete sine transform.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • type (integer) – the type of DST. Must be in [1; 4].

  • axes (list of integers) – The axes along which the transform is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : make transform orthogonal and divide by sqrt(N)
    2 : divide by N

    where N is the product of n_i for every transformed axis i. n_i is 2*(<axis_length>+1 for type 1 and 2*<axis length> for types 2, 3, 4. Making the transform orthogonal involves the following additional steps for every 1D sub-transform:

    Type 1

    nothing

    Type 2

    divide last output value by sqrt(2)

    Type 3

    multiply last input value by sqrt(2)

    Type 4

    nothing

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.genuine_fht(a: numpy.ndarray, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a full Hartley transform. A full forward Fourier transform is carried out over the requested axes, and the real part minus the imaginary part of the result is stored in the output array. For a single transformed axis, this is identical to separable_fht, but when transforming multiple axes, the results are different.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the transform is carried out. If not set, all axes will be transformed.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed axes.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.genuine_hartley(a: numpy.ndarray, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a full Hartley-like transform. A full forward Fourier transform is carried out over the requested axes, and the sum of real and imaginary parts of the result is stored in the output array. For a single transformed axis, this is identical to separable_hartley, but when transforming multiple axes, the results are different.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the transform is carried out. If not set, all axes will be transformed.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed axes.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

Notes

This function uses a nonstandard Hartley convention. Only use if you know exactly what you are doing!

ducc0.fft.good_size()

Returns a good length to pad an FFT to.

Parameters:
  • n (int) – Minimum transform length

  • real (bool, optional) – True if either input or output of FFT should be fully real.

Returns:

out – The smallest fast size >= n

Return type:

int

ducc0.fft.r2c(a: numpy.ndarray, axes: object = None, forward: bool = True, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs an FFT whose input is strictly real.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the FFT is carried out. If not set, this is assumed to be list(range(a.ndim)). The real-to-complex transform will be executed along axes[-1], and will be executed first.

  • forward (bool) – If True, a negative sign is used in the exponent, else a positive one.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed input axes.

  • out (numpy.ndarray (complex type with same accuracy as a)) – For the required shape, see the Returns section. Must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data. The shape is identical to that of the input array, except for axes[-1]. If the length of that axis was n on input, it is n//2+1 on output.

Return type:

numpy.ndarray (complex type with same accuracy as a)

ducc0.fft.r2r_fftpack(a: numpy.ndarray, axes: object, real2hermitian: bool, forward: bool, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a real-valued FFT using FFTPACK’s halfcomplex storage scheme.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the FFT is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • real2hermitian (bool) – if True, the input is purely real and the output will have Hermitian symmetry and be stored in FFTPACK’s halfcomplex ordering, otherwise the opposite.

  • forward (bool) – If True, a negative sign is used in the exponent, else a positive one.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the length of axis.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data. The shape is identical to that of the input array.

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.r2r_fftw(a: numpy.ndarray, axes: object, forward: bool, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a real-valued FFT using FFTW’s halfcomplex storage scheme.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the FFT is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • forward (bool) – If True, a negative sign is used in the exponent, else a positive one.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the length of axis.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data. The shape is identical to that of the input array.

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.separable_fht(a: numpy.ndarray, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a separable Hartley transform. For every requested axis, a 1D forward Fourier transform is carried out, and the intermediate result is its real part minus its imaginary part. Then the next axis is processed.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the transform is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed axes.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

ducc0.fft.separable_hartley(a: numpy.ndarray, axes: object = None, inorm: int = 0, out: object = None, nthreads: int = 1) numpy.ndarray

Performs a separable Hartley-like transform. For every requested axis, a 1D forward Fourier transform is carried out, and the real and imaginary parts of the result are added before the next axis is processed.

Parameters:
  • a (numpy.ndarray (any real type)) – The input data

  • axes (list of integers) – The axes along which the transform is carried out. If not set, this is assumed to be list(range(a.ndim)). Axes will be transformed in the specified order.

  • inorm (int) –

    Normalization type
    0 : no normalization
    1 : divide by sqrt(N)
    2 : divide by N

    where N is the product of the lengths of the transformed axes.

  • out (numpy.ndarray (same shape and data type as a)) – May be identical to a, but if it isn’t, it must not overlap with a. If None, a new array is allocated to store the output.

  • nthreads (int) – Number of threads to use. If 0, use the system default (typically the number of hardware threads on the compute node).

Returns:

The transformed data

Return type:

numpy.ndarray (same shape and data type as a)

Notes

This function uses a nonstandard Hartley convention. Only use if you know exactly what you are doing!