Fourier Transform and Convolution

This module provides some functions to compute fourier transform, convolution and correlation numerically.

Fourier transform

These functions compute fourier transform or its inverse for signals defined on an interval or region which is symmetric w.r.t. the origin. They differ from vanilla Discrete Fourier Transform (DFT) provided by torch.fft in that correct scale and shift are considered and thus serve as numerical approximation to continuous fourier transform.

For example, the function ft1() computes the 1D Fourier transform of signal \(f\). Given signal length \(N\), assuming \(f[k]\) represents sampled values on points \(\{k\delta_f\}_{k=-\lfloor N/2\rfloor}^{\lfloor(N-1)/2\rfloor}\) where \(\delta_f\) is sampling interval, then ft1() computes

\[\ft\{f\}[l]=h[l]=\sum_{k=-\lfloor N/2\rfloor}^{\lfloor(N-1)/2\rfloor} f[k]\e^{-\i 2\pi k\delta_f\times l\delta_h}\delta_f, l=-\lfloor\frac{N}{2}\rfloor,\cdots,\lfloor\frac{N-1}{2}\rfloor\]

where \(\delta_h=1/(N\delta_f)\) is sampling interval in frequency domain. Indices \(-\lfloor\frac{N}{2}\rfloor, \cdots,\lfloor\frac{N-1}{2}\rfloor\) for \(k\) and \(l\) correspond to 0, …, N-1 in the given array. In other words, this function works like

>>> from torch.fft import fft, fftshift, ifftshift
>>>
>>> f = torch.rand(9)
>>> h1 = ft1(f, 0.1)
>>> h2 = fftshift(fft(ifftshift(f))) * 0.1
>>> torch.allclose(h1, h2)
True

Multiplication with sampling interval is necessary if consistent scale with continuous Fourier transform is desired, which is called DFT-scale. Otherwise, it can be set to None.

The following example illustrates its approximation to continuous Fourier transform, as well as its precision in a simple occasion for 32-bit float.

>>> x = dnois.utils.interval(1000, 4e-3)  # [-2, 2]
>>> y = torch.exp(-torch.pi * x.square())  # Gaussian function
>>> fx = dnois.utils.interval(1000, 0.25)  # frequency
>>> g = torch.exp(-torch.pi * fx.square())  # ground truth spectrum
>>> torch.allclose(ft1(y, 4e-3).real, g, atol=1e-6)
True

Note

The sampling interval (like delta argument for ft1()) will be multiplied to transformed signal if given, so it can be a float or a tensor with shape broadcastable to original signal. But if it is not a 0D tensor, its size on the dimensions to be transformed must be 1.

ft1(f[, delta, dim])	Computes the 1D Fourier transform of signal f.
ift1(h[, delta, dim])	Computes the 1D inverse Fourier transform of spectrum h.
ft2(f[, dx, dy, dims])	Computes the 2D Fourier transform of signal f.
ift2(h[, dx, dy, dims])	Computes the 2D inverse Fourier transform of signal h.

Convolution

These functions compute convolution between two signals f1 and f2. dconv() computes discrete convolution. Numerical optimization and correct scale are considered in conv() so it serve as a numerical approximation to continuous convolution.

>>> f1, f2 = torch.tensor([1., 2., 3., 4.]), torch.tensor([1., 2., 3.])
>>> g_full_linear = torch.tensor([1., 4., 10., 16., 17., 12.])
>>> torch.allclose(dconv(f1, f2), g_full_linear)
True
>>> g_same_linear = torch.tensor([4., 10., 16., 17.])
>>> torch.allclose(dconv(f1, f2, out='same'), g_same_linear)
True
>>> g_full_circular = torch.tensor([18., 16., 10., 16.])
>>> torch.allclose(dconv(f1, f2, padding='none'), g_full_circular)
True

>>> f1, f2 = torch.rand(5), torch.rand(6)
>>> torch.allclose(conv(f1, f2, spacing=0.1), dconv(f1, f2) * 0.1)
True

Note

If all the signals involved is real-valued, set real to True to improve the computation and return a real-valued tensor.

conv(f1, f2[, dims, spacing, out, padding, ...])	Computes continuous convolution for two tensors \(f_1\) and \(f_2\) utilizing FFT.
conv1(f1, f2[, dx, dim, out, padding, ...])	1D version of conv().
conv2(f1, f2[, dx, dy, dims, out, padding, ...])	2D version of conv().
dconv(f1, f2[, dims, out, padding, real])	Computes n-D discrete convolution for two tensors \(f_1\) and \(f_2\) utilizing FFT.
dconv1(f1, f2[, dim, out, padding, real])	1D version of dconv().
dconv2(f1, f2[, dims, out, padding, real])	2D version of dconv().

Multiple convolutions

Sometimes one may need to compute convolutions between a group of tensors and a single common tensor. In that case, it is sensible to compute the Fourier transform of the common tensor only once. The following functions serve for this purpose.

dconv_mult(f1s, f2, dims[, out, padding, real])	Discrete version of conv_mult().
dconv_partial(f1, f2_ft, f2_shape, dims, real)	Discrete version of conv_partial().
ft4dconv(f1_shape, f2, dims, real[, padding])	Discrete version of ft4conv().
conv_mult(f1s, f2, dims[, spacing, out, ...])	Compute the convolutions between each signal in f1s and f2.
conv_partial(f1, f2_ft, f2_shape, dims, real)	See conv_mult().
ft4conv(f1_shape, f2, dims, real[, spacing, ...])	See conv_mult().