FIELD OF THE INVENTION

This disclosure generally relates to signal processing and, in particular, to determination of a spectrum, i.e., frequency components, of a signal that is sparse in the frequency domain.

BACKGROUND

The Fast Fourier Transform (FFT) computation reduces the computational cost of computing the Discrete Fourier Transform (DFT) of a general complex N-vector from O(N2) to O(N log(N)). The general complex N-vector may include N discrete values of a signal, such as a signal received from a sensor such as a radar, and infra-red sensor, etc., a signal received from or transmitted using an antenna, etc. Since its popularization in the 1960s, the FFT algorithm has played a crucial role in multiple areas including scientific computing, signal processing, and computer science. In the general case, such scaling from O(N2) to O(N log(N)) is at most a factor log(N) from optimality. In more restricted cases however, such as when the signal/vector to be recovered is sparse, it is possible, at least in theory, to significantly improve on the latter.

Indeed, the past decade or so has seen the design and study of various algorithms that can compute the DFT of sparse vectors using significantly less time and measurements than those required by a traditional FFT. That is, if ƒ is an N×1 vector corresponding to the DFT of an N×1 vector {circumflex over (ƒ)} containing at most R<<N nonzero elements, it is possible to recover {circumflex over (ƒ)} using significantly fewer samples than the traditional “Nyquist rate” («O(N)) and in computational complexity much lower than that of the FFT («O(N log(N))). These schemes are generally referred to as “sparse Fast Fourier Transform” (sFFT) algorithms, and they generally fall within two categories: 1) deterministic versus randomized algorithms, and 2) exactly sparse versus approximately sparse algorithms.

In general, for a periodic signal,

f⁡(x)=∑j=0N-1⁢e-2⁢π⁢⁢ixj⁡(f^j+η⁢v^j),(1)
we say that the signal is R-sparse if the support of the spectrum {circumflex over (ƒ)}, i.e., the set of indices for which a frequency component (also called frequency coefficient) {circumflex over (ƒ)}jis “non-zero,” has cardinality smaller than or equal to R. Equation (1) describes a time-domain signal ƒ(x) in terms of its frequency coefficients {circumflex over (ƒ)} and, as such, may be referred to as inverse DFT. According to one convention, the sign of the exponent representing an inverse DFT or inverse FFT is positive.

Recovery of the location and magnitude of the nonzero coefficients is referred to as the exactly R-sparse FFT problem if n=0 and as the noisy R-sparse FFT problem if 0<n is relatively small. The value of a frequency coefficient is considered to be “zero” if the absolute value is less than or equal to a specified small number, such as 0.0002, 0.1, 0.6, 1, 1.5, 3, etc. Frequency coefficients having a value greater than or equal to a specified threshold are called “non-zero” coefficients. When the signal is not R-sparse but nonetheless compressible, the problem of finding the best R-sparse approximation is referred to as the approximately sparse FFT problem and takes the following form: Given an accuracy parameter e, find ƒ* such that,

f^*-f^a≤(1+ɛ)⁢miny^:R-sparse⁢y^-f^b+η⁢f^c(2)
where ∥·∥a, ∥·∥band ∥·∥carep-norms (generally p=1 or 2).

Of the former category, randomized algorithms have had the most success in practice thus far; although some deterministic algorithms do exhibit quasilinear complexity in R and polylogarithmic complexity in N, the algorithmic constant and the exponents are often so large that the methods are only competitive when N is impractically large or when the sparsity satisfies stringent conditions. On the other hand, both the algorithmic constant and the complexity of randomized sparse FFTs are in general much smaller. The most recent results on the topic are shown in Table 1 inFIG. 1.

Among approximate sparse FFT algorithms, some randomized techniques have achieved the best complexities in the average and worst-case scenarios. They use different techniques to locate the support of the important frequencies, such as orthogonal frequency division multiplexing (OFDM) where the problem is “lifted” from 1D to 2D and sparsity along each column/row is leveraged. Another technique uses binning techniques (shifting, shuffling and filtering) to produce a support-locating estimator. These techniques ultimately rely in general on estimating some complex phase containing information with regard to the support of the signal.

Among exact sparse FFT algorithms, one technique achieved the best complexities in the worst case, and a deterministic technique and a randomized technique with constant probability of success achieved the best complexity in the average case. We note however that these techniques can be very unstable, even to low-level noise, being based upon the Chinese Remainder Theorem. The performance of the other technique can degrade significantly when a worst case scenario is considered. In addition, to achieve such low complexity, one technique makes use of adaptive sampling, where samples for a subsequent iterations are chosen based on information available at the current iteration. These results show an interesting trade-off: in an R-sparse worst-case scenario, algorithms exhibiting the best sampling/computational complexity often possess scaling much inferior than their counter-part per regards to the probability of failure p, and vice-versa. Such characteristics can be detrimental in settings where high-performance is needed and failure cannot be tolerated. Finally, we also observe from Table 1 that whenever the algorithm is generalizable to multiple dimensions, the algorithmic constant exhibits a scaling which is at least exponential in the dimension d.

SUMMARY

Various techniques are described for treating a noisy R-sparse FFT problem when the spectrum {circumflex over (ƒ)} is nonnegative. Various embodiments or implementations for determining the spectrum {circumflex over (ƒ)} of a signal ƒ(x), i.e., the frequency coefficients of the signal {circumflex over (ƒ)}(x), exhibit a low sampling complexity of(R log(R)log(N)) and a low computational complexity of(R log2(R)log(N)) (whereindicates the presence of logc(log(·)) factors), and further scales according to log(p) with respect to the probability of failure p, thus alleviating the issues underlined above. In addition, various embodiments can perform with constant scaling with respect to more than one dimensions, i.e., in d dimensions with N=Mdunknowns, the scaling is of the form log(N)=d log(M) without further dependence. In this context, a low sampling complexity may mean 1/10th, 1/25th, 1/100th, 1/1000thetc., of the traditional O(N) Nyquist-rate samples, or even fewer samples. N can also be the total number of frequency coefficients of the signal ƒ(x). Similarly, a low computational complexity may mean requiring 1/10th, 1/25th, 1/100th, 1/1000thetc., of the number of arithmetic and/or data-access operations that are required in a traditional FFT computation or in other sFFT computations. The algorithm is also non-adaptive. Therefore, sampling and computations can be interleaved, but this is not necessary. The required sampling can be performed in a data-collection phase, and the sets of signal samples can be used subsequently to determine the frequency coefficients of the signal.

Various embodiments use techniques such as shuffling and filtering. However, our support-locating scheme does not rely on sample-based phase estimation; it works in Fourier space rather than in sample space. Furthermore, rather than necessarily performing simultaneously the operations locating coefficients and computing the values of the coefficients of the “heavy” (i.e., non-zero) frequencies, we can separate both processes. According to various embodiments, locating the support without attempting to estimate values accurately is a less complex task that can be carried out more rapidly and with higher probability of success than attempting to perform both tasks at once. As for the treatment of the high-dimensional problem, we reduce any d-dimensional problem to a 1-dimensional problem without substantial overhead cost through the use of rank-1 lattices, and demonstrate the correctness of this approach.

Computation of Values when Support Setis Known

Accordingly, in one aspect, a method is provided for identifying magnitudes of frequency components of a signal ƒ(x). The method includes obtaining a first group of sets of signal samples, where each set contains a number of samples of the signal ƒ(x). The cardinality of each set is equal to or greater than an upper bound (R) of a number of non-zero frequency components of the signal ƒ(x), and less than a number of samples N of the signal ƒ(x) according to a Nyquist rate. N can also be the total number of frequency coefficients of the signal ƒ(x). The cardinality can be prime number. The method also includes computing by a processor an initial set of frequency coefficients {circumflex over (ƒ)}0using a support setand the first group of sets of signal samples. The support setidentifies indices of the non-zero frequency coefficients of the signal ƒ(x).

Obtaining the first group of sets of signal samples may include accessing from memory T sets of signal samples, where T is on the order of O(1) or on the order of O(R logRN). A t-th set of signal samples, wherein index t ranges from 1 through T, may include P(t)signal samples, where P(t)is greater than or equal to R and smaller than N. In some cases, P(t)can be a prime number, and/or a number less than R. Samples collected and stored in memory prior to the execution of various computations steps may be accessed from the memory.

Obtaining the first group of sets of signal samples may include determining the number of samples (N) of the signal ƒ(x) according to the Nyquist rate, and obtaining the upper bound (R) of the number of non-zero frequency components of the signal. The number N itself is determined or estimated; N actual samples of the signal ƒ(x) may not be obtained. In addition, the method includes selecting a number of sets (T), where T is on the order of O(1) or on the order of O(R logRN), and selecting T recovery parameters P(t), wherein index t ranges from 1 through T. One or more recovery parameters P(t)can be prime numbers. The method further includes, for each recovery parameter P(t), sampling the signal ƒ(x) at P(t)distinct locations, providing the first group of T sets of signal samples, where each sample set has P(t)samples.

In some embodiments, computing the initial set of frequency coefficients {circumflex over (ƒ)}0includes multiplying the first group of sets of signal samples by: (i) a discrete Fourier transform (DFT) matrix F having T DFT submatrices, and (ii) a binary matrix B based on the support set. The size of each submatrix is P(t). Each P(t)can be a prime number, and index t ranges from 1 through T. In some embodiments, the method includes generating a set of frequency coefficients {circumflex over (ƒ)} using a binary matrix B based on the support setand the initial set of frequency coefficients {circumflex over (ƒ)}0.

The method may also include determining that a norm based on a binary matrix B based on the support setand the initial set of frequency coefficients {circumflex over (ƒ)}0is not less than a selected threshold. If the norm is greater than or equal to the selected threshold, the method includes obtaining a second group of sets of signal samples, each set in the second group also having a number of samples of the signal ƒ(x). Cardinality of each set in the second group is a number less than N; the cardinality can be greater than R and/or a prime number. If the norm is greater than or equal to the selected threshold, the method includes re-computing the initial set of frequency coefficients {circumflex over (ƒ)}0using the support setand the second group of sets of signal samples.

In some embodiments, the method includes determining that a norm based on the binary matrix B and the re-computed initial set of frequency coefficients {circumflex over (ƒ)}0is less than or equal to the selected threshold, and generating a set of frequency coefficients {circumflex over (ƒ)} using the binary matrix B and the re-computed initial set of frequency coefficients {circumflex over (ƒ)}0. The signal ƒ(x) may be received from one or more signal source including a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

In another aspect, a system is provided for identifying magnitudes of frequency components of a signal ƒ(x). The system includes a first processor and a first memory in electrical communication with the first processor. The first memory includes instructions that can be executed by a processing unit including the first processor or a second processor, or both. The processing unit may be in electronic communication with a memory module that includes the first memory or a second memory or both.

The instructions in the first memory program the processing unit to: obtain a first group of sets of signal samples, where each set contains a number of samples of the signal ƒ(x). The cardinality of each set is equal to or greater than an upper bound (R) of a number of non-zero frequency components of the signal ƒ(x), and less than a number of samples N of the signal ƒ(x) according to a Nyquist rate. N can also be the total number of frequency coefficients of the signal ƒ(x). The cardinality can be prime number. The instructions also program the processing unit to compute an initial set of frequency coefficients {circumflex over (ƒ)}0using a support setand the first group of sets of signal samples. The support setidentifies indices of the non-zero frequency coefficients of the signal ƒ(x). In various embodiments, the instructions can program the processing unit to perform one or more of the method steps described above. In some embodiments, the system includes one or more signal samplers to obtain signal samplers, and one or more sample shufflers and/or filters to filter a shuffled set of samples.

In another aspect, an article of manufacture is provided that includes a non-transitory storage medium having stored therein instructions which, when executed by a processing unit program the processing unit, which is in electronic communication with a memory module, to identifying magnitudes of frequency components of a signal ƒ(x). The instructions program the processing unit to: obtain a first group of sets of signal samples, where each set contains a number of samples of the signal ƒ(x). The cardinality of each set is equal to or greater than an upper bound (R) of a number of non-zero frequency components of the signal ƒ(x), and less than a number of samples N of the signal ƒ(x) according to a Nyquist rate. N can also be the total number of frequency coefficients of the signal ƒ(x). The cardinality can be prime number. The instructions also program the processing unit to compute an initial set of frequency coefficients {circumflex over (ƒ)}0using a support setand the first group of sets of signal samples. The support setidentifies indices of the non-zero frequency coefficients of the signal ƒ(x). In various embodiments, the instructions can program the processing unit to perform one or more of the method steps described above.

Finding Aliased Supportkwhenk−1ork−1are Known

In another aspect, a method is provided for identifying frequency components of a signal ƒ(x). The method includes performing by a processor the following steps (a) through (d):

(a) initially designating a current set of candidate support coefficients (k) as a current set of aliased support coefficientsk;

(b) obtaining a first set of K shuffled samples of the signal ƒ(x) corresponding to a first sampling interval based on a first shuffling parameter, where K is a fraction of the number of samples (N) of the signal ƒ(x) according to a Nyquist rate;

(c) filtering the shuffled samples in the first set, and computing a first group of frequency coefficients of the shuffled samples in the first set. Each frequency coefficient corresponds to a respective candidate support coefficient in the current setk; and

(d) removing from the current set of aliased support coefficientska subset of candidate support coefficients wherein, for each candidate support coefficient in the subset a value of a corresponding computed frequency coefficient in the first group of frequency coefficients is less than a threshold.

In some embodiments, the method further includes obtaining a second set of shuffled samples of the signal ƒ(x), corresponding to a second sampling interval based on a second shuffling parameter. In these embodiments, the method also includes filtering the shuffled samples in the second set, computing a second group of frequency coefficients of the shuffled samples in the second set, and removing from the current set of aliased support coefficientska subset of candidate support coefficients where for each candidate support coefficient in the subset, i.e., each candidate support that this removed, a value of a corresponding computed frequency coefficient in the second plurality of frequency coefficients is less than the threshold. These steps may be repeated for additional sets of shuffled samples of the signal ƒ(x), corresponding to respective sampling intervals based on respective shuffling parameters.

In some embodiments, the first shuffling parameter is a first coprime of an index limit (Mk) of the current set of candidate support coefficients (k), and the second shuffling parameter is a second, different coprime of the index limit (Mk). The first and/or the second shuffling parameters can be prime numbers. An index limit (Mk) of the current set of candidate support coefficients (k) may be associated with a support set growth factor ρkthat is a ratio of an index limit (Mk) and the number of samples in the first set of shuffled samples (K). The method may further include obtaining a prior set of aliased support coefficientsk−1, and dealiasing the prior set of aliased support coefficients to obtain the current set of candidate support coefficients (k).

In some embodiments, the method further includes the following steps:

(e) selecting a next index limit (Mk+1) of a next set of candidate support coefficients (k+1), where the next index limit is greater than a current index limit (Mk) of the current set of candidate support coefficients (k);

(f) after the step (d), dealiasing the current set of aliased support coefficients (k) using the next index limit (Mk+1), to obtain the next set of candidate support coefficients (k+1);

(g) updating the current set of candidate support coefficients (k) by designating the next set of candidate support coefficients (k+1) as the current set of candidate support coefficients (k); and

(h) repeating the steps (a) through (d) using the updated current set of candidate support coefficients (k).

The method may also include determining that a next index limit (Mk+1) is not less than a number of samples (N) of the signal ƒ(x) according to a Nyquist rate; and after the step (h), designating the current set of aliased support coefficientskas a final set of aliased support coefficients.

Obtaining the first set of shuffled samples of the signal ƒ(x) may include selecting the first sampling interval using the first shuffling parameter (Qk(l)), wherein the first shuffling parameter (Qk(l)) corresponds to an index limit (Mk) of the current set of candidate support coefficients (k). Obtaining the first set of shuffled samples may also include sampling the signal ƒ(x) using the first sampling interval, to obtain a first sampled signal, and shuffling the first sampled signal using the first shuffling parameter and the index limit. Additionally or in the alternative, obtaining the first set of shuffled samples of the signal ƒ(x) may include selecting the first shuffling parameter (Qk(l)) corresponding to an index limit (Mk) of the current set of candidate support coefficients (k), and accessing from memory the first set of shuffled samples corresponding to the first shuffling parameter.

Filtering the shuffled samples in the first set may include low-pass filtering. In some embodiments, the low-pass filtering includes filtering the first set using a Gaussian filter having a standard deviation (a) that is based on an upper bound of a number of non-zero frequency components (R) of the signal ƒ(x), to obtain a filtered shuffled signal.

In some embodiments, the method include generating the current set of candidate support coefficients by selecting a first index limit (M1), where M1is less than the number of samples (N) of the signal ƒ(x) according to the Nyquist rate. Generating the current set of candidate support coefficients may also include obtaining several samples of the signal ƒ(x), determining M1frequency coefficients of a signal formed using the several samples, and including each non-zero frequency coefficient from the M1frequency coefficients as a respective candidate support coefficient in the current set of candidate support coefficients. A non-zero frequency coefficient may include a frequency coefficient having a magnitude (also called an absolute value) greater than a specified threshold value. A frequency coefficient designated as a zero coefficient may have a magnitude less than or equal to the specified threshold value.

In another aspect, a system is provided for identifying frequency components of a signal ƒ(x). The system includes a first processor and a first memory in electrical communication with the first processor. The first memory includes instructions that can be executed by a processing unit including the first processor or a second processor, or both. The processing unit may be in electronic communication with a memory module that includes the first memory or a second memory or both.

The instructions in the first memory program the processing unit to:

(a) initially designate a current set of candidate support coefficients (k) as a current set of aliased support coefficientsk,

(b) obtain a first set of K shuffled samples of the signal ƒ(x) corresponding to a first sampling interval based on a first shuffling parameter, where K is a fraction of the number of samples (N) of the signal ƒ(x) according to a Nyquist rate;

(c) filter the shuffled samples in the first set, and computing a first group of frequency coefficients of the shuffled samples in the first set. Each frequency coefficient corresponds to a respective candidate support coefficient in the current setk; and

(d) remove from the current set of aliased support coefficientka subset of candidate support coefficients wherein, for each candidate support coefficient in the subset a value of a corresponding computed frequency coefficient in the first group of frequency coefficients is less than a threshold.

In various embodiments, the instructions can program the processing unit to perform one or more of the method steps described above. In some embodiments, the system includes one or more signal samplers to obtain signal samplers, and one or more sample shufflers and/or filters to filter a shuffled set of samples.

In another aspect, an article of manufacture is provided that includes a non-transitory storage medium having stored therein instructions which, when executed by a processing unit program the processing unit, which is in electronic communication with a memory module, to identify frequency components of a signal ƒ(x). To this end, the instructions program the processing unit to perform one or more operations according to the method steps and/or the system described above.

Identifying Non-Zero Frequency Components of One and Multi-Dimensional Signals and Determining Values of the Non-Zero Frequency Components

In another aspect, a method is provided for generating non-zero frequency coefficients of a signal ƒ(x). The method includes performing by a processor the steps of: obtaining a group of sets of samples of the signal ƒ(x), where the group includes a first set and a different second set. The method also includes determining indices of the non-zero frequency coefficients of the signal ƒ(x) using at least the first set, and determining values of the non-zero frequency coefficients of the signal ƒ(x) using the indices thereof and at least the second set.

Determining the indices of the non-zero frequency coefficients may include obtaining: (i) a number of samples (N) of the signal ƒ(x) according to a Nyquist rate, and (ii) an upper bound of a number of non-zero frequency components (R) of the signal ƒ(x). Only the number N may be determined or obtained, without actually collecting N samples. Determining the indices of the non-zero frequency coefficients may also include selecting a sample size K, a number of iterations P, one or more support set growth factors ρk, and one or more index limits Mk, where each index limit represents an index limit of a respective set of candidate support coefficientsk, where Mk=ρkK and

N=K⁢∏k=1P⁢ρk.
Determining the indices of the non-zero frequency coefficients may further include determining a current set of aliased support coefficientkusing the first set of samples of the signal ƒ(x), the first set having less than N samples.

Determining the indices of the non-zero frequency coefficients may include performing one or more iterations, where at least one iteration includes obtaining a current set of candidate support coefficientskby dealiasing the current set of aliased support coefficientsk, and determining a next set of aliased support coefficients using the current set of candidate support coefficientskand one or more sets of samples from the group of sets of samples of the signal ƒ(x). Each of the one or more sets of samples may be obtained using a sampling interval based on Mk. Determining the indices of the non-zero frequency coefficients may also include designating the next set of aliased support coefficients as the current set of aliased support coefficientsk, for the next iteration.

Determining the indices of the non-zero frequency coefficients may include performing at least one iteration that includes obtaining a current set of candidate support coefficientskby dealiasing the current set of aliased support coefficientsk, and determining a next set of aliased support coefficients using the current set of candidate support coefficientskand one or more sets of samples from the group of sets of samples of the signal ƒ(x). Each one of the one or more sets of samples may be obtained using a sampling interval based on Mk. Determining the indices of the non-zero frequency coefficients may also include designating the next set of aliased support coefficients as final set of aliased support coefficients, where the final set of aliased support coefficients represents the indices of the non-zero frequency coefficients.

The signal ƒ(x) may include a one-dimensional signal obtained from one or more signal sources including a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor. In some embodiments, the method includes mapping a multi-dimensional signal y(x) into a one-dimensional signal ƒ(x), and storing an inverse mapping identifying a sample index for y(x) corresponding to a sample index for ƒ(x). Obtaining the group of sets of samples of the signal ƒ(x) may include selecting a number of sample indices of ƒ(x), and obtaining samples of y(x) using sample indices of y(x), each one being determined via the inverse mapping of the respective selected sample index of ƒ(x).

The non-zero frequency coefficients of the signal ƒ(x) may represent non-zero frequency coefficients of the multi-dimensional signal y(x). The multi-dimensional signal y(x) may be obtained from one or more signal sources including a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor. A non-zero frequency coefficient may include a frequency coefficient having a magnitude (also called an absolute value) greater than a specified threshold value. A frequency coefficient designated as a zero coefficient may have a magnitude less than or equal to the specified threshold value.

In some embodiments, determining the indices of the non-zero frequency coefficients is performed using a first value (N1) of a number of samples of the signal ƒ(x) according to a Nyquist rate, to obtain a first support set1. The step of determining the indices of the non-zero frequency coefficients is repeated using a second value (N2) of the number of samples, to obtain a second support set2. Only the values N1and N2are used; N1and/or N2samples of the signal ƒ(x) are not collected or used. In these embodiments, the method may also include generating a final support set as a union of the first and second support sets, and designating coefficients of the final support set as the indices of the non-zero frequency coefficients. The step of determining the indices of the non-zero frequency coefficients may be repeated using one or more additional values (Nz) of the number of samples, to obtain one or more additional support sets, and by computing the union of all support sets.

In another aspect, a system is provided for generating non-zero frequency coefficients of a signal ƒ(x). The system includes a first processor and a first memory in electrical communication with the first processor. The first memory includes instructions that can be executed by a processing unit including the first processor or a second processor, or both. The processing unit may be in electronic communication with a memory module that includes the first memory or a second memory or both. The instructions in the first memory program the processing unit to: obtain a group of sets of samples of the signal ƒ(x), where the group includes a first set and a different second set, determine indices of the non-zero frequency coefficients of the signal ƒ(x) using at least the first set, and determine values of the non-zero frequency coefficients of the signal ƒ(x) using the indices thereof and at least the second set.

In various embodiments, the instructions can program the processing unit to perform one or more of the method steps described above. In some embodiments, the system includes one or more signal samplers to obtain signal samplers, and one or more sample shufflers and/or filters to filter a shuffled set of samples.

In another aspect, an article of manufacture is provided that includes a non-transitory storage medium having stored therein instructions which, when executed by a processing unit program the processing unit, which is in electronic communication with a memory module, to generate non-zero frequency coefficients of a signal ƒ(x). To this end, the instructions program the processing unit to perform one or more operations according to the method steps and/or the system described above.

