An Efficient High-Dimensional Sparse Fourier Transform¹¹1This manuscript has been submitted to IEEE Transactions on Aerospace and Electronic Systems for reviewing on Sept. 1st, 2016.

We use lower-case (upper-case) bold letters to denote vectors (matrices). $(\cdot)^{T}$ and $(\cdot)^{H}$ respectively denote the transpose and conjugate transpose of a matrix or a vector. $\|\cdot\|$ is Euclidean norm for a vector. $[\mathbf{a}]_{i}$ is the $i_{th}$ element of vector $\mathbf{a}$. All operations on indices in this paper are taken modulo $N$ denoted by $[\cdot]_{N}$. We use $\lfloor\cdot\rfloor$ to denote rounding to floor. $[S]$ refers to the set of indices $\{0,...,S-1\}$, and $[S]\setminus a$ is for eliminating element $a$ from set $[S]$. We use $\{0,1\}^{B}$ to denote the set of $B$-dimensional binary vectors. We use $\mathop{\mathrm{diag}}\nolimits(\cdot)$ to denote forming a diagonal matrix from a vector and use $\mathbb{E}\{\cdot\}$ to denote expectation. The DFT of signal ${\mathbf{s}}$ is denoted as $\hat{{\mathbf{s}}}$. We also assume that the signal length in each dimension is an integer power of 2.

This paper is organized as follows. A brief background on the SFT algorithm is given in Section II. Details of the proposed RSFT algorithm are given in Section III. Section IV presents the derivation of the optimal threshold design for the RSFT algorithm. Then in Section V we provide some numerical results to verify the theoretical findings. An application of the RSFT algorithm in radar signal processing is presented in Section VI and finally, concluding remarks are made in Section VII.

As opposed to the FFT which computes the coefficients of all $N$ discrete frequency components of an $N$-sample long signal, the SFT[2] computes only the $K$ frequency components of a $K$-sparse signal. Before outlining the SFT algorithm we provide some key definitions and properties, which is extracted and reformulated based on [2] .

The permutation has the following property.

At a high level, the SFT algorithm in [2] runs two loops, namely the Location loop and the Estimation loop. The former finds the indices of the $K$ most significant frequencies from the input signal, while the latter estimates the corresponding Fourier coefficients. Here, we emphasize on Location more than Estimation, since the former is more relevant to the radar application that we consider. The Location step provides frequency locations, which in the radar case is directly related to target parameters.

In the Location loop, a permutation procedure reorders the input data in the time domain, causing the frequencies to also reorder. The permutation causes closely spaced frequencies to appear in well separated locations with high probability. Then, a flat-window[2] is applied on the permuted signal for the purpose of extending a single frequency into a (nearly) boxcar, for a reason that will become apparent in the following. The windowed data are aliased, as in Definition 2. The frequency domain equivalent of this aliasing is undersampling by $N/B$ (see Property 2). The flat-window used at the previous step ensures that no peaks are lost due to the effective undersampling in the frequency domain. After this stage, a FFT of length $B$ is employed.

The permutation and the aliasing procedure effectively map the signal frequencies from $N$-dimensional space into a reduced $B$-dimensional space, where the first-stage-detection procedure finds the significant frequencies’ peaks, and then their indices are reverse mapped into the original $N$-dimensional frequency space. However, the reverse mapping yields not only the true location of the significant frequencies, but also ${N}/{B}$ ambiguous locations for each frequency. To remove the ambiguity, multiple iterations of Location with randomized permutation are performed. Finally, the second-stage-detection procedure locates the $K$ most significant frequencies from the accumulated data for each iteration. More details about the SFT algorithm can be found in [2].

In real world applications, the frequencies are continuous and can take any value in $[0,2\pi)$. When fitting a grid on these frequencies, leakage occurs from off-grid frequencies, which can diminish the sparsity of the signal. As the leakage due to strong frequency components can mask the contributions of weak frequency components, it becomes difficult to determine the frequency domain peaks after permutation. (see Fig. 1 (c)). To address this problem, we propose to multiply the received time domain signal with a window before permutation. We call this procedure pre-permutation windowing. The idea is to confine the leakage within a finite number of frequency bins, as illustrated in Fig. 1.

The choice of the pre-permutation window is determined by the required resolution, computational complexity, and degree of leakage suppression. Specifically, the side-lobe level should be lower than the noise level after windowing (see Fig. 4). However, the larger the attenuation of the side-lobes, the wider the main lobe would be, thus lowering the frequency resolution. Meanwhile, a broader main lobe results in increased computational load, which will be discussed in Section IV-D.

In the SFT, detection of the significant frequencies is needed in two stages. With knowing the number of the significant frequencies and assuming they are all on-grid, the detection of the signal can be accomplished by finding the $K$ highest spectral amplitude values.

In reality however, we usually do not have the knowledge of the exact sparsity, i.e., $K$. Moreover, even if we knew $K$, due to the leakage caused by the off-grid signals, the $K$ highest spectral peaks might not be the correct representation of the signal frequencies. Finally, the additive noise could generate false alarms, which would add more difficulties in signal detection.

In order to solve the detection problem, we propose to use NP detection in the two stages of detection, which does not require knowing $K$. However, the optimal thresholds still rely on the number of significant frequencies as well as their SNR. In light of this, we provide an optimal design based on a bound for the sparsity and signal SNR level. The optimal threshold design is presented in Section IV.

In the following, we elaborate on the high dimensional extension of the RSFT for its main stages.

We model the signal in continuous time domain as a superposition of $K$ sinusoids as well as additive white noise. We then sample the signal uniformly both in I and Q channels with a sampling frequency above the Nyqvist rate. Assume the total sampling time is divided into $T$ consecutive equal length segments, each of which contains $N$ samples, and $K<<N$ (i.e., the signal is sparse in frequency domain). Then for each time segment, i.e., $s\in[T]$, we have

where $\mathbf{v}(\omega_{i})$ denotes for the $i_{th}$ complex sinusoid, with $\omega_{i}\in[0,2\pi)$ as its frequency, i.e.,

We further assume that $\omega_{i}$ is unknown deterministic quantity and is constant during the whole process, while the complex amplitude of the sinusoid, i.e., $b_{i,s}$ takes random value for each segment. More specifically, we model $b_{i,s}$ as a circularly symmetric Gaussian process with the distribution $b_{i,s}\sim\mathcal{CN}(0,\sigma_{bi}^{2})$. Likewise, the noise $\mathbf{n}_{s}$ is distributed as $\mathbf{n}_{s}\sim\mathcal{CN}(\mathbf{0},\sigma_{n}^{2}\mathbf{I})$, where $\mathbf{0}$ is $N$-dimensional zero vector, and $\mathbf{I}\in\mathbb{R}^{N\times N}$ is unit matrix. We also assume each sinusoid and the noise are uncorrelated. In addition, the neighboring sinusoids are resolvable in the frequency domain, i.e., the frequency spacing of neighboring sinusoids is greater than $\eta_{m}\frac{2\pi}{N}$, where $\eta_{m}$ is the 6.0-dB bandwidth[16] of a window that applies on $\mathbf{r}_{s}$. Note the signal model in (10) is also commonly used in array signal processing literature for Uniform Linear Array (ULA) settings (See e.g., [17]), in which the time samples in each segment is replaced by spatial samples from array elements, and a sample segment is referred as a time snapshot.

We want to detect and estimate each $\omega_{i}$ from the input signal. From a non-parametric and data-independent perspective, this is a classic spectral analysis and detection problem that can be solved by any FFT-based spectrum estimation method, for example, the Bartlett method followed by a NP detection procedure (see Appendix A). In that case, the FFT computes the signal spectrum on $N$ frequency bins, and the detection is carried on each bin to determine whether there exists a significant frequency. In what follows, we define the detection and estimation problem related to the design of the RSFT.

Let $SNR_{i}=\sigma_{bi}^{2}/\sigma_{n}^{2}$ be the SNR of the $i_{th}$ sinusoid. Let us define the worst case SNR, i.e., $SNR_{min}$, as

Let $P_{d}$ denote for the probability of detection for the sinusoid with $SNR_{min}$, and $P_{fa}$ the corresponding probability of false alarm on each frequency bin.

In the following, we investigate the two stages of detection separately, then summarize the solution into an optimization problem.

The first stage detection is performed on each data segment. After pre-permutation windowing, permutation and flat-windowing, the input signal can be expressed as

where $\sigma_{s}$ is the permutation parameter for the $s_{th}$ segment, which has an uniform random distribution; $\mathbf{P}_{\sigma_{s}}$ is the permutation matrix, which functions as (2) with $\tau=0$; $\mathbf{W}=\mathop{\mathrm{diag}}\nolimits(\mathbf{w})$, $\overline{\mathbf{W}}=\mathop{\mathrm{diag}}\nolimits(\bar{\mathbf{w}})$, where $\mathbf{w}$ and $\bar{\mathbf{w}}$ are pre-permutation window and flat-window, respectively.

Regarding the design of the flat-window, we take its frequency domain main-lobe to have width $2\pi/B$, and choose its length in time domain as $N$. As indicated in [2], it is possible to use less data in flat-windowing by choosing the length of $\bar{\mathbf{w}}$ less than $N$, i.e., dropping some samples in each segment after the permutation. However, a reduced length window in the time domain would result in longer transition regions in the frequency domain and as a result the detection performance of the system would degrade, since a larger transition region would allow more noise to enter the estimation process.

The time domain aliasing can be described as

where $L=N/B$; $\overline{\mathbf{W}}_{i}$ is the $i_{th}$ sub-matrix of $\overline{\mathbf{W}}$, which is comprised of the $iB_{th}$ to the $((i+1)B-1)_{th}$ rows of $\overline{\mathbf{W}}$. And $\mathbf{V}_{\sigma_{s}}=\sum_{i\in[L]}\overline{\mathbf{W}}_{i}\mathbf{P}_{\sigma_{s}}\mathbf{W}$.

The FFT operation on the aliased data $\mathbf{f}$ can be expressed as

where $\mathbf{D}\in\mathbb{C}^{B\times B}$ is the DFT matrix. For the $k_{th}$ entry of $\hat{\mathbf{f}}$, we have

where $\mathbf{u}_{k}$ is the $k_{th}$ column of $\mathbf{D}$, i.e., $\mathbf{u}_{k}=[1\quad e^{jk\Delta\omega_{B}}\;\cdots\;e^{jk(B-1)\Delta\omega_{B}}]^{T}$, and $\Delta\omega_{B}=2\pi/B$.

Substituting (10) into (16), and taking the $m_{th}$ sinusoid, which we assume is the weakest sinusoid with its SNR equals to $SNR_{min}$, out of the summation

Since $[\hat{\mathbf{f}}]_{k}$ is a linear combination of $b_{i},[\mathbf{n}]_{j},i\in[K],j\in[N]$, it holds that

where

and

It is easy to see that $\sigma_{fk}^{2}$ is summation of weighted variance from each signal and noise component.

We now investigate the $p_{th}$ bin, where $\omega_{m}$ is mapped to, and we have the following claim.

On assuming that only $\omega_{m}$ maps to bin $p$, and that the side-lobes (leakage) are far below the noise level (owning to the two stages of windowing, which attenuate the leakage down to a desired level), the effect of leakage from other sinusoids can be ignored. Then we can approximate the variance of $[\hat{\mathbf{f}}]_{p}$ as

In case that multiple frequencies are mapped to the same bin (collision), (22) gives a underestimate of the variance. The probability of a collision occurring reduces as $K<<B$.

The bin $u\in[B]$, for which no significant frequency is mapped to, contains only noise, and the corresponding variance for $[\hat{\mathbf{f}}]_{u}$ is

Hence, the hypothesis test for the first-stage-detection on $[\hat{\mathbf{f}}]_{j},j\in[B]$ is formulated as

• •

  $\mathcal{H}0$: no significant frequency is mapped to it.

• •

  $\mathcal{H}1$: at least one significant frequency is mapped to it, with worst case SNR equals to $SNR_{min}$.

The log likelihood ratio test (LLRT) is

where $P_{f_{j}|H1}(x)$ and $P_{f_{j}|H0}(x)$ are the probability density function (PDF) of $[\hat{\mathbf{f}}]_{j}$ under $\mathcal{H}1$ and $\mathcal{H}0$ respectively, and $\gamma^{\prime}$ is a threshold.

Substituting the PDF of $[\hat{\mathbf{f}}]_{j}$ under both hypothesis into (24), and after some manipulations we get

Hence, $|[\hat{\mathbf{f}}]_{j}|^{2}$ is a sufficient statistics for first stage detection. Since $[\hat{\mathbf{f}}]_{j}$ has circularly symmetric Gaussian distribution, $|[\hat{\mathbf{f}}]_{j}|^{2}$ is exponentially distributed with cumulative distribution function (CDF)

where $\zeta^{2}$ equals to $\sigma^{2}_{f_{u}}$ under ${\mathcal{H}0}$ and $\sigma^{2}_{f_{p}}$ under ${\mathcal{H}1}$.

Based on (26), in the first stage of detection, the false alarm rate on each of $B$ bins and the probability of detection of the weakest sinusoid can be derived to be equal to

where $\gamma$ is the detection threshold. Both $\tilde{P}_{fa}$ and $\tilde{P}_{d}$ depend on the permutation $\sigma_{s}$. Taking expectation with respect to $\sigma_{s}$, we have

where $\bar{P}_{d}(\omega_{m})=\mathbb{E}\{\tilde{P}_{d}(\sigma_{s},\omega_{m})\}$, $\bar{P}_{fa}=\mathbb{E}\{\tilde{P}_{fa}(\sigma_{s})\}$, $\bar{\alpha}(p,\omega_{m})=\mathbb{E}\{\tilde{\alpha}(p,\omega_{m},\sigma_{s})\}$, and $\bar{\beta}=\mathbb{E}\{\tilde{\beta}(\sigma_{s})\}$.

Let $\mathbf{c}_{\sigma_{s}}\in\{0,1\}^{B}$ denote the output of the first-stage-detection for the $s_{th}$ segment, with permutation factor $\sigma_{s}$. Each entry in $\mathbf{c}_{\sigma_{s}}$ is a Bernoulli random variable, i.e., for $j\in[B]$,

Note that under $\mathcal{H}1$, we assume that $[\mathbf{c}_{\sigma_{s}}]_{j}$ corresponds to the weakest sinusoid. For the other $K-1$ co-existing sinusoids, since their SNR may be grater than $SNR_{min}$, their probability of detection may also be grater than $\tilde{P}_{d}(\omega_{m},\sigma_{s})$ (see Claim 3).

The reverse-mapping stage hashes the $B$-dimensional $\mathbf{c}_{\sigma_{s}}$ back to the $N$-dimensional $\mathbf{a}_{\sigma_{s}}$. According to Definition 4, it holds that

After accumulation of $T$ iterations, each entry in the accumulated output is summation of $T$ Bernoulli variables with different success rate. Define $\bar{\mathbf{a}}$ as the accumulated output, then for its $i_{th},i\in[N]$ entry, we have

Note that in (31), each term inside the sum corresponds to a different segment, i.e., $[\mathbf{c}_{\sigma_{s}}]_{j}$ is from the $s_{th}$ segment. Since $\sigma_{s}$ is drawn randomly for each segment, $j$ may take different values, and relates to $i$ via a reverse-mapping. Fig. 5 gives a graphical illustration of the mapping and reverse-mapping.

Now, the hypothesis test for the second-stage-detection on $[\bar{\mathbf{a}}]_{i},i\in[N]$ is formulated as

• •

  $\overline{\mathcal{H}}0$: no significant frequency exists.

• •

  $\overline{\mathcal{H}}1$: there exists a significant frequency, whose SNR is at least $SNR_{min}$.

In the following, we investigate the statistics of $[\bar{\mathbf{a}}]_{i}$ under both hypothesis in an asymptotic senses. Before that however, we will take a closer look at the mapping and the reverse mapping by providing the following properties.

The proofs of these properties are provided in Appendix B. The two properties simply reveal the following facts: a mapped location can be recovered by reverse mapping (with ambiguities). Also, when applying reverse mapping to two distinct locations with the same permutation parameter, the resulting locations are also distinct.

Under $\overline{\mathcal{H}}1$, assuming that $[\bar{\mathbf{a}}]_{i}$ corresponds to the $m_{th}$ sinusoid, i.e., the weakest sinusoid, then each term inside the sum of (31) has distribution $[\mathbf{c}_{\sigma_{s}}]_{j}\sim\mathrm{Bernoulli}\left(\tilde{P}_{d}(\omega_{m},\sigma_{s})\right),s\in[T]$. Then we present the following claim.

The distribution of $[\bar{\mathbf{a}}]_{i}$ under $\overline{\mathcal{H}}0$ is more complicated, and we have following claim.

A natural extension of Remark 1 is Remark 2, which gives the condition under which the RSFT will reach its limit.

Based on the above discussion, the optimal threshold design in Problem 1 can be solved by the following optimization problem, i.e.,

where $g_{a_{0}}(u),g_{a_{1}}(u)$ are the asymptotic PDF²²2We take the upper bounds of the variances in both distributions. It is shown in Section V-E that the actual variances is close to their upper bounds. of $[\bar{\mathbf{a}}]_{i}$ (which corresponds the weakest sinusoid) under $\overline{\mathcal{H}}0$ and $\overline{\mathcal{H}}1$, respectively. Since both of them are Normal distributions, with fixed threshold, i.e., $\mu$, we can solve for $\bar{P}_{d}(\omega_{m}),\bar{P}_{fa}$, and then compute the $SNR_{min}$. By enumerating $\mu\in[T]$, the minimum worst case SNR, i.e., $SNR_{min}^{\ast}$ can be found, and the corresponding $\mu^{\ast}$ is the optimal threshold for the second stage of detection. The optimal threshold for the first stage of detection, i.e., $\gamma^{\ast}$, can thus be calculated via (28).

By averaging over the permutation, asymptotically, $SNR_{min}^{\ast}$ does not depend on the permutation. However, it still depends on $\omega_{m}$, and we have the following claim to manifest their relationship.