Sparse Fourier transform

The sparse Fourier transform (SFT) is a kind of discrete Fourier transform (DFT) for handling big data signals. Specifically, it is used in GPS synchronization, spectrum sensing and analog-to-digital converters.:

The fast Fourier transform (FFT) plays an indispensable role on many scientific domains, especially on signal processing. It is one of the top-10 algorithms in the twentieth century. However, with the advent of big data era, the FFT still needs to be improved in order to save more computing power. Recently, the sparse Fourier transform (SFT) has gained a considerable amount of attention, for it performs well on analyzing the long sequence of data with few signal components.

Definition
Consider a sequence xn of complex numbers. By Fourier series, xn can be written as

x_n=(F^*X)_n=\sum_{k=0}^{N-1}X_k e^{j\frac{2\pi}{N}kn}. $$ Similarly, Xk can be represented as

X_k=\frac{1}{N}(Fx)_k=\frac{1}{N}\sum_{n=0}^{N-1}x_n e^{-j\frac{2\pi}{N}kn}. $$ Hence, from the equations above, the mapping is $$F:\mathbb C^N\to \mathbb C^N$$.

Single frequency recovery
Assume only a single frequency exists in the sequence. In order to recover this frequency from the sequence, it is reasonable to utilize the relationship between adjacent points of the sequence.

Phase encoding
The phase k can be obtained by dividing the adjacent points of the sequence. In other words,

\frac{x_{n+1}}{x_n}=e^{j\frac{2\pi}{N}k} = \cos\left(\frac{2\pi k}{N}\right)+j \sin\left(\frac{2\pi k}{N}\right). $$ Notice that $$x_n \in \mathbb C^N$$.

An aliasing-based search
Seeking phase k can be done by Chinese remainder theorem (CRT).

Take $$k=104{,}134$$ for an example. Now, we have three relatively prime integers 100, 101, and 103. Thus, the equation can be described as

k=104{,}134\equiv \left\{ \begin{array}{rl} 34 & \bmod 100, \\ 3 & \bmod 101, \\ 1 & \bmod 103. \end{array} \right. $$ By CRT, we have

k=104{,}134\bmod (100\cdot101\cdot103)=104{,}134\bmod 1{,}040{,}300 $$

Randomly binning frequencies
Now, we desire to explore the case of multiple frequencies, instead of a single frequency. The adjacent frequencies can be separated by the scaling c and modulation b properties. Namely, by randomly choosing the parameters of c and b, the distribution of all frequencies can be almost a uniform distribution. The figure Spread all frequencies reveals by randomly binning frequencies, we can utilize the single frequency recovery to seek the main components.



x_n'=X_k e^{j\frac{2\pi}{N}(c\cdot k+b)}, $$ where c is scaling property and b is modulation property.

By randomly choosing c and b, the whole spectrum can be looked like uniform distribution. Then, taking them into filter banks can separate all frequencies, including Gaussians, indicator functions, spike trains,    and Dolph-Chebyshev filters. Each bank only contains a single frequency.

The prototypical SFT
Generally, all SFT follows the three stages

Identifying frequencies
By randomly bining frequencies, all components can be separated. Then, taking them into filter banks, so each band only contains a single frequency. It is convenient to use the methods we mentioned to recover this signal frequency.

Estimating coefficients
After identifying frequencies, we will have many frequency components. We can use Fourier transform to estimate their coefficients.

X_k'=\frac 1 L \sum_{l=1}^L x_n'e^{-j\frac{2\pi}{N}n'\ell} $$

Repeating
Finally, repeating these two stages can we extract the most important components from the original signal.

x_n-\sum_{k'=1}^k X_k' e^{j\frac{2\pi}{N}k'n} $$

Sparse Fourier transform in the discrete setting
In 2012, Hassanieh, Indyk, Katabi, and Price proposed an algorithm that takes $$ O ( k \log n \log (n/k) ) $$ samples and runs in the same running time.

Sparse Fourier transform in the high dimensional setting
In 2014, Indyk and Kapralov proposed an algorithm that takes $$ 2^{O(d \log d)} k \log n $$ samples and runs in nearly linear time in $$ n $$. In 2016, Kapralov proposed an algorithm that uses sublinear samples $$ 2^{O(d^2)} k \log n \log \log n $$ and sublinear decoding time $$ k \log^{O(d)} n $$. In 2019, Nakos, Song, and Wang introduced a new algorithm which uses nearly optimal samples $$ O( k \log n \log k ) $$ and requires nearly linear time decoding time. A dimension-incremental algorithm was proposed by Potts, Volkmer  based on sampling along rank-1 lattices.

Sparse Fourier transform in the continuous setting
There are several works about generalizing the discrete setting into the continuous setting.

Implementations
There are several works based on MIT, MSU, ETH and Universtity of Technology Chemnitz [TUC]. Also, they are free online.
 * MSU implementations
 * ETH implementations
 * MIT implementations
 * GitHub
 * TUC implementations