Non-uniform discrete Fourier transform

Definition

The nonuniform discrete Fourier transform (NUDFT) transforms a sequence of N complex numbers x_0, \ldots, x_{N-1} into another sequence of complex numbers X_0, \ldots, X_{N-1} defined by where p_0, \ldots, p_{N-1} \in [0,1] are sample points and f_0, \ldots, f_{N-1} \in [0,N] are frequencies. Note that if p_n = n/N and f_k = k, then equation () reduces to the discrete Fourier transform.

The NUDFT defines a linear operator. Fast numerical algorithms for evaluating this operator are known as the nonuniform fast Fourier transform (NUFFT).

Terminology and conventions

Two related but distinct conventions are used in the literature to classify nonuniform Fourier transforms:

• In harmonic analysis and signal processing, NUDFTs are commonly described from a series-evaluation viewpoint, in which transform types are distinguished by whether the sample locations or the frequencies are nonuniform.

• In numerical analysis and scientific computing, particularly in the context of NUFFT algorithms and software libraries, transforms are classified by the mapping direction between uniform and nonuniform representations, reflecting linear-operator structure and adjoint relationships.

NUDFT types

• The nonuniform discrete Fourier transform of type I (NUDFT-I) uses uniform sample points p_n = n/N but nonuniform (i.e. non-integer) frequencies f_k. This corresponds to evaluating a generalized Fourier series at equispaced points. It is also known as NDFT.. It is sometimes called forward NDFT and corresponds to NUFFT Type II (uniform to nonuniform) in the operator-based convention.
• The nonuniform discrete Fourier transform of type II (NUDFT-II) uses uniform (i.e. integer) frequencies f_k = k but nonuniform sample points p_n. This corresponds to evaluating a Fourier series at nonequispaced points. It is also known as adjoint NDFT or NUFFT Type I (nonuniform to uniform).
• The nonuniform discrete Fourier transform of type III (NUDFT-III) uses both nonuniform sample points p_n and nonuniform frequencies f_k. This corresponds to evaluating a generalized Fourier series at nonequispaced points. It is also known as NNDFT or NUFFT Type III (nonuniform to nonuniform).

Adjoint relationship

In the operator-based convention in NUFFT, the Type I and Type II transforms are adjoint to each other in the sense of linear operators. Specifically, the NUFFT Type I transform is the Hermitian adjoint of the NUFFT Type II transform with respect to the standard inner products on the corresponding discrete spaces, up to scaling factors.

A similar family of NUDFTs can be defined by substituting -i for +i in equation (). Unlike in the uniform case, however, this substitution is unrelated to the inverse Fourier transform. The inversion of the NUDFT is a separate problem, discussed below.

Multidimensional NUDFT

The multidimensional NUDFT converts a d-dimensional array of complex numbers x_{\mathbf{n}} into another d-dimensional array of complex numbers X_{\mathbf{k}} defined by X_{\mathbf{k}} = \sum_{\mathbf{n}=\mathbf{0}}^{\mathbf{N}-1} x_{\mathbf{n}} e^{-2\pi i \mathbf{p}{\mathbf{n}} \cdot \boldsymbol{f}{\mathbf{k}}} where \mathbf{p}{\mathbf{n}} \in [0,1]^d are sample points, \boldsymbol{f}{\mathbf{k}} \in [0,N_1] \times [0,N_2] \times \cdots \times [0,N_d] are frequencies, and \mathbf{n} = (n_1,n_2,\ldots,n_d) and \mathbf{k} = (k_1,k_2,\ldots,k_d) are d-dimensional vectors of indices from 0 to The multidimensional NUDFTs of types I, II, and III are defined analogously to the 1D case.

Relationship to Z-transform

The NUDFT-I can be expressed as a Z-transform. The NUDFT-I of a sequence x[n] of length N is X(z_k)=X(z)|{z=z_k}=\sum{n=0}^{N-1}x[n]z_k^{-n},\quad k=0, 1, \dots, N-1,

where X(z) is the Z-transform of x[n], and {z_i}_{i=0}^{N-1} are arbitrarily distinct points in the z-plane. Note that the NUDFT reduces to the DFT when the sampling points are located on the unit circle at equally spaced angles.

Expressing the above as a matrix, we get \mathbf{X} = \mathbf{D} \mathbf{x}

where \mathbf{X} = \begin{bmatrix} X(z_0) \ X(z_1) \ \vdots \ X(z_{N-1}) \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x[0] \ x[1] \ \vdots \ x[N-1]\end{bmatrix}, and \mathbf{D} = \begin{bmatrix} 1 & z_0^{-1} & z_0^{-2} & \cdots & z_0^{-(N-1)} \ 1 & z_1^{-1} & z_1^{-2} & \cdots & z_1^{-(N-1)} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & z_{N-1}^{-1} & z_{N-1}^{-2} & \cdots & z_{N-1}^{-(N-1)} \end{bmatrix}.

Direct inversion of the NUDFT-I

As we can see, the NUDFT-I is characterized by \mathbf{D} and hence the N {z_k} points. If we further factorize \det(\mathbf{D}), we can see that \mathbf{D} is nonsingular provided the N {z_k} points are distinct. If \mathbf{D} is nonsingular, we can get a unique inverse NUDFT-I as follows: \mathbf{x}=\mathbf{D^{-1}}\mathbf{X}.

Given \mathbf{X}\text{ and }\mathbf{D}, we can use Gaussian elimination to solve for \mathbf{x}. However, the complexity of this method is O(N^3). To solve this problem more efficiently, we first determine X(z) directly by polynomial interpolation: \hat X [k] = X(z_k), \quad k=0, 1, \dots, N-1.

Then x[n] are the coefficients of the above interpolating polynomial.

Expressing X(z) as the Lagrange polynomial of order N-1, we get X(z) = \sum_{k=0}^{N-1} \frac{L_k(z)}{L_k(z_k)} \hat X [k],

where \left{L_i(z)\right}_{i=0}^{N-1} are the fundamental polynomials:

L_k(z) = \prod_{i \ne k} \left(1 - z_i z^{-1}\right), \quad k = 0, 1, \dots, N-1.

Expressing X(z) by the Newton interpolation method, we get X(z) = c_0 + c_1 \left(1 - z_0 z^{-1}\right) + c_2 \left(1 - z_0 z^{-1}\right) \left(1 - z_1 z^{-1}\right) + \cdots + c_{N-1} \prod_{k=0}^{N-2} \left(1 - z_k z^{-1}\right),

where c_j is the divided difference of the j-th order of \hat X [0], \hat X [1], \dots, \hat X [j] with respect to z_0, z_1, \dots, z_j: \begin{align} c_0 &= \hat X [0], & c_1 &= \frac{\hat X [1] - c_0}{1 - z_0 z_1^{-1}}, \ c_2 &= \frac{\hat X [2] - c_0 - c_1 \left(1 - z_0 z^{-1}\right)}{\left(1 - z_0 z_2^{-1}\right) \left(1 - z_1 z_2^{-1}\right)}, & & \cdots \end{align}

The disadvantage of the Lagrange representation is that any additional point included will increase the order of the interpolating polynomial, leading to the need to recompute all the fundamental polynomials. However, any additional point included in the Newton representation only requires the addition of one more term.

We can use a lower triangular system to solve {c_j}: \mathbf{L}\mathbf{c} = \mathbf{X} where \mathbf{X} = \begin{bmatrix} \hat X [0] \ \hat X [1] \ \vdots \ \hat X [N-1] \end{bmatrix}, \qquad \mathbf{c} = \begin{bmatrix} c_0 \ c_1 \ \vdots \ c_{N-1} \end{bmatrix}, and \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \ 1 & \left(1 - z_0 z_1^{-1}\right) & 0 & \cdots & 0 \ 1 & \left(1 - z_0 z_2^{-1}\right) & \left(1 - z_0 z_2^{-1}\right) \left(1 - z_1 z_2^{-1}\right) & \cdots & 0 \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & \left(1 - z_0 z_{N-1}^{-1}\right) & \left(1 - z_0 z_{N-1}^{-1}\right) \left(1 - z_1 z_{N-1}^{-1}\right) & \cdots & \prod_{k=0}^{N-2} \left(1 - z_k z_{N-1}^{-1}\right) \end{bmatrix}.

By the above equation, {c_j} can be computed within O(N^2) operations. In this way Newton interpolation is more efficient than Lagrange interpolation unless the latter is modified by L_{k+1}(z) = \frac{1 - z_{k+1} z^{-1}}{1 - z_k z^{-1}} L_k(z),\quad k=0, 1, \dots, N-1.

Nonuniform fast Fourier transform

While a naive application of results in an O(N^2) algorithm for computing the NUDFT, O(N \log N) algorithms based on the fast Fourier transform (FFT) do exist. Such algorithms are referred to as NUFFTs or NFFTs and have been developed based on oversampling and interpolation, min-max interpolation, In general, NUFFTs leverage the FFT by converting the nonuniform problem into a uniform problem (or a sequence of uniform problems) to which the FFT can be applied. Software libraries for performing NUFFTs are available in 1D, 2D, and 3D.

Applications

The applications of the NUDFT include:

• Digital signal processing
• Magnetic resonance imaging
• Numerical partial differential equations
• Semi-Lagrangian schemes
• Spectral methods
• Spectral analysis
• Digital filter design
• Antenna array design
• Detection and decoding of dual-tone multi-frequency (DTMF) signals

References

References

1. (1999). "The Nonuniform Discrete Fourier Transform and Its Applications in Signal Processing". Springer US.
2. (February 2003). "Nonuniform fast fourier transforms using min-max interpolation". IEEE Transactions on Signal Processing.
3. (June 2005). "The type 3 nonuniform FFT and its applications". Journal of Computational Physics.
4. (January 2004). "Accelerating the Nonuniform Fast Fourier Transform". SIAM Review.
5. (2019). "Numerical Fourier Analysis". Birkhäuser.
6. "Basic use of PyNUFFT — PyNUFFT 2023.2.2 documentation".
7. "Mathematical definitions of transforms — finufft 2.2.0 documentation".
8. (2001). "Nonuniform Sampling: Theory and Practice". Springer.
9. Dutt, Alok. (May 1993). "Fast Fourier Transforms for Nonequispaced Data". Yale University.
10. (November 1993). "Fast Fourier Transforms for Nonequispaced Data". SIAM Journal on Scientific Computing.
11. (January 2003). "Fast Summation at Nonequispaced Knots by NFFT". SIAM Journal on Scientific Computing.
12. (December 1992). "A fast algorithm for Chebyshev, Fourier, and sinc interpolation onto an irregular grid". Journal of Computational Physics.
13. (20 February 2018). "A Nonuniform Fast Fourier Transform Based on Low Rank Approximation". SIAM Journal on Scientific Computing.
14. "NUFFT page".
15. "NFFT".
16. (2019-02-13). "MikaelSlevinsky/FastTransforms.jl".
17. (2019-02-07). "chebfun/chebfun".