Nonequispaced Fourier Transform

Sseziwa Mukasa,
Software Engineer
JEOL USA
mukasa@jeol.com

Problem Description

The nonequispaced Fourier transform solves the sum

(Subscript[f, j])=f((Subscript[ν, j]))=(∑_(k∈(Subscript[I, N]^d)))(Subscript[(OverHat[f]), k])(e^(-2 π i k (Subscript[ν, j]))) (j∈(Subscript[I, M]^1))

or its transpose

(Subscript[(OverHat[h]), k])=(∑_(j∈(Subscript[I, M]^1)))(Subscript[f, j]) (e^(2 πi k (Subscript[ν, j]))) (k∈(Subscript[I, N]^d))

where (Subscript[ν, j])are M d-dimensional frequencies with ((-1)/2)≤(Subscript[ν, j])<(1/2)for each component of (Subscript[ν, j]) and (Subscript[I, N]^d)={k∈(ℤ^d);-(N/2)≤k<(N/2)}.  For (N^d) = M, and (Subscript[ν, j])=(j/N)(j∈(Subscript[I, N]^d)) the first sum is equivalent to the discrete Fourier transform.  For arbitrary (Subscript[ν, j]), algorithms for approximating the sums to accuracy ε in 𝒪((N^d)logN+(log^d)((1/ε)) M) time are given in [1].  This notebook demonstrates the use of a Mathematica package nfft.m which implements both the direct and approximate methods of solving the two sums.  The approximate methods will hereafter be refered to as nfft and nfftTransposed for the first and second sums respectively.

NFFT

The functions nfft and nfftTranspose are contained in a Mathematica package file nfft.m which should be placed in the appropriate location for your Mathematica environment.

<<nfft`

One application of the nfft if is for solving Fourier like sums which occur in the Filter Diagonalization Method [2].  Each of the diagonal entries of the evolution matrices in the Filter Diagonalization method is a sum of the form:

(∑_((OverVector[n])=(OverVector[0]))^(2 (OverVector[M])))(e^(i (OverVector[n]).(OverVector[φ]))) c((OverVector[n])+(OverVector[p]))

where c is a regularly sampled d-dimensional data set and (OverVector[φ]) is a basis point in a d-dimensional space with -π ≤ φ < π componentwise.  Letting ν=(φ/(2 π)) and padding c with zeroes to the length 4 (OverVector[M])+1 so the sum can be taken over [-(OverVector[M]),(OverVector[M])] then the sum is of the form of a nonequispaced Fourier transform assuming the basis functions do not form a regular lattice, in which case a regular fast Fourier transform can be employed.  As a concrete example the diagonal elements of a random 8x8 data set at the four basis points (-π,-π), (-π,0), ((π/2),-π), ((π/2),0) can be computed as follows:

ν={{-0.5`,-0.5`},{-0.5`,0},{0.25`,-0.5`},{0.25`,0}};
c=PadLeft[RandomReal[{0,100000},{7,7}],{13,13}];
result=nfft[c,ν]

The messages are due to singularities in the default window function and can be safely ignored, window functions will be described in greater detail later.  It should be emphasized that nfft is an approximation, which can be compared with a direct calculation of the sum (∑_(k∈(Subscript[I, N]^d)))(Subscript[(OverHat[f]), k])(e^(-2 π i k (Subscript[ν, j]))) (j∈(Subscript[I, M]^1)) by using the routine ndft.  For example

With[{x=result-ndft[c,ν]},Sqrt[x.Conjugate[x]]]

(NFFT^T)

The routine nfftTranspose computes the value of the sum (∑_(j∈(Subscript[I, M]^1)))(Subscript[f, j]) (e^(2 πi k (Subscript[ν, j]))) (k∈(Subscript[I, N]^d)).  The NFFT sum can be seen as a matrix multiplication f=A (OverHat[f]) where (Subscript[A, jk])=(e^(-2 π i (j-(N/2)) (Subscript[ν, k]))) (j∈0..N-1).  In that case the sum (∑_(j∈(Subscript[I, M]^1)))(Subscript[f, j]) (e^(2 πi k (Subscript[ν, j]))) (k∈(Subscript[I, N]^d)) can then be seen as the matrix equation (OverHat[f])=(A^T) f.  As an example of the use of nfftTranspose consider a decaying signal sampled at 39 exponentially increasing intervals.

samplepoints=FoldList[Plus,0,Table[(δ^i),{i,0,39}]/.FindRoot[Sum[(δ^(i )),{i,0,n-1}]==255/.{n->40},{δ,2}]];
data=Exp[2 π I (0.25+0.01 I) samplepoints];

relabeling the sample points as

ν=-(samplepoints/510);

and computing the transposed nonequispaced Fourier transform at 512 points

result=nfftTranspose[data,Transpose[{ν}],{512}];

results in an approximation of the sum result[[j]]=(∑_(j∈(Subscript[I, 512]^1)))data[[i]] (e^(-2 π i j (Subscript[ν, k]))) which is similar to the Fourier transform of the same signal sampled at regular intervals: (Subscript[(OverHat[f]), k])=(∑_(j=0)^255)(Subscript[f, j]) (e^(-2 π i j (k/512))) (k∈{-128..127}).  The two can be compared with:

<<"ErrorBarPlots`";<<"PlotLegends`"
ListPlot[{Transpose[{(1/512) Range[-256,254,2],Re[RotateRight[Fourier[Exp[2 π I (0.25`+0.01` I) Range[0,255]],FourierParameters->{1,-1}],128]]}],Transpose[{(1/512) Range[-256,255],Re[result]}]},Joined->True,PlotStyle->{Hue[0,1,0],Hue[0]},PlotMarkers->None,PlotRange->All]

where the red line is the result of nfftTranspose.  It should be noted in general that the transposed NFFT is not the inverse of the NFFT unlike the relationship of the discrete Fourier transform and inverse discrete Fourier transform.  This is because when viewed as a matrix multiplication, the matrix A of the nonequispaced Fourier transform is rectangular and thus the expression (A^T) fis generally over or underdetermined.  A strategy for computing inverse nonequispaced Fourier transforms is discussed in [3] but not implemented in the current Mathematica package.

For the sake of comparison a version of the transposed NFFT which does not use an approximation scheme is also implemented as ndftTranspose, the difference between the approximation and explicit calculation for the transposed NFFT example can be computed as follows:

Sqrt[Total[# Conjugate[#]&[result-ndftTranspose[data,Transpose[{ν}],{512}]]]]

Approximation Windows

nfft and nfftTranspose are both approximate routines.  The approximation itself is accomplished by using a window function which is localized in both the time and frequency domains [1].  The quality of the approximation is determined by the options to nfft and nfftTranspose: approximationWidth, oversamplingFactor, precomputedψ, and precomputedc.  approximationWidth and oversamplingFactor control the truncation in the time domain of the window function and amount of oversampling in the time domain used to compute the NFFT.  The precomputedψ and precomputedc contain lists representing the value of the window function in the time and frequency domain at the appropriate nodes.  They are best computed using the precomputeψ and precomputec functions.

precomputedψ contains the values of the window function ψ evaluated at the nodes (Subscript[ν, j])-(l/n) where n (Subscript[ν, j])-m≤l≤n (Subscript[ν, j])+m, n is Max[Dimensions[data]*oversamplingFactor,2*m+2] where m is the approximation width.  precomputedψ also contains the bounds on l for each (Subscript[ν, j]).  precomputedc contains the value of the Fourier transform of ψ evaluated at nodes ⌊((-N)/2)⌋...⌊(N/2)⌋.  The functions precomputeψ and precomputec can be used to compute appropriate lists as follows where the utility function precomputen is used to determine the size of the oversampled data set as follows:

nfftTranspose[data,Transpose[{ν}],{512},precomputedψ->precomputeψ[Transpose[{ν}],precomputen[{512}]],precomputedc->precomputec[precomputen[{512}]]];

The window functions used in the approximation can be overriden with the windowFunction and windowFunctionTransformed options of precomputeψ and precomputec.  The windowFunction is expressed as a function of ν, the size of the approximation grid, half the width of the approximation grid minus one, and the oversampling factor.  windowFunctionTransformed is a function of the index k, the size of the approximation grid, half the width of the approximation grid minus one, and the oversampling factor.  For example to use a Guassian window function, define

gaussian[ν_,n_,m_,α_]:=With[{b=((2 α m)/(π (2 α-1)))},((Exp[((-(((n ν))^2))/b)])/(Sqrt[π b]))]
gaussianTransformed[k_,n_,m_,α_]:=With[{b=((2 α m)/(π (2 α-1)))},Exp[-(((((π k)/n)))^2) b]]

This difference between using this window and the default for the same window width and oversampling rate can be computed as:

result2=nfftTranspose[data,Transpose[{ν}],{512},precomputedψ->precomputeψ[Transpose[{ν}],precomputen[{512}],windowFunction->gaussian],precomputedc->precomputec[precomputen[{512}],windowFunctionTransformed->gaussianTransformed]];
ListPlot[Transpose[{(1/512) Range[-256,255],Re[result-result2]}],Joined->True,PlotRange->All]
With[{x=result-result2},(Sqrt[x.Conjugate[x]])]

Bibliography

[1] D. Potts,  G. Steidl, and  M. Tasche, "Fast Fourier transforms for nonequispaced data: A tutorial" in Modern Sampling Theory: Mathematics and Applications, J.J. Benedetto and P. Ferreira (Eds.), Chapter 12, pages 249-274.

[2] V.A.Mandelshtam,The Multidimensional Filter Diagonalization Method,I.Theory and Numerical Implementation,J.Mag.Res.144 (2000) 343-356.

[3] S. Kunis and D. Potts.Stability results for scattered data interpolation by trigonometric polynomials.  Preprint, Univ. Lübeck, A-04-12, 2004.