Systems and method for determining frequency coefficients of signals

Abstract

A system for determining the frequency coefficients of a one or multi-dimensional signal that is sparse in the frequency domain includes determining the locations of the non-zero frequency coefficients, and then determining values of the coefficients using the determined locations. If N is total number of frequency coefficients across the one or more dimension of the signal, and if R is an upper bound of the number of non-zero ones of these frequency coefficients, the systems requires up to (O (R log(R) (N))) samples and has a computation complexity of up to O (R log.sup.2(R) log (N). The system and the processing technique are stable to low-level noise and can exhibit only a small probability of failure. The frequency coefficients can be real and positive or they can be complex numbers.

Claims

1. A method for identifying magnitudes of frequency components of a signal ƒ(x), the method comprising: obtaining at a receiver a first plurality of sets of samples of a signal received by an antenna or a sensor, each set comprising a plurality of samples of the signal ƒ(x) obtained using a respective sampling interval based on a respective shuffling parameter, cardinality of each set being 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, wherein obtaining each set of signal samples comprises subsampling by the receiver the signal at a respective rate less than the Nyquist rate, and at least one set of signal samples determines a support set ; and computing by a processor an initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 using the support set and frequency coefficients of the first plurality of sets of signal samples, the support set identifying indices of the non-zero frequency coefficients of the signal ƒ(x).

2. The method of claim 1, wherein obtaining the first plurality of sets of signal samples comprises accessing from memory T sets of signal samples, wherein: T is on the order of O(1) or on the order of O(R log.sub.R N); and a t-th set of signal samples, wherein index t ranges from 1 through T, comprises P.sup.(t) signal samples, P.sup.(t) being greater than R and smaller than N.

3. The method of claim 1, wherein obtaining the first plurality of sets of signal samples comprises: 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; selecting a number of sets (T), wherein T is on the order of O(1) or on the order of O(R log.sub.R N), and selecting T recovery parameters P.sup.(t), wherein index t ranges from 1 through T; for each recovery parameter P.sup.(t), sampling the signal ƒ(x) at P.sup.(t) distinct locations, providing the first plurality of T sets of signal samples, each sample set having P.sup.(t) samples.

4. The method of claim 1, wherein computing the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 comprises multiplying the first plurality of sets of signal samples by: (i) a discrete Fourier transform (DFT) matrix F comprising T DFT submatrices, each submatrix having a size P.sup.(t), and index t ranging from 1 through T, and (ii) a binary matrix B based on the support set .

5. The method of claim 1, further comprising: generating a set of frequency coefficients {circumflex over (ƒ)} using a binary matrix B based on the support set and the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0.

6. The method of claim 1, further comprising: determining that a norm based on a binary matrix B based on the support set and the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 is not less than a selected threshold; obtaining a second plurality of sets of signal samples, each set comprising a plurality of samples of the signal ƒ(x); and re-computing the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 using the support set and the second plurality of sets of signal samples.

7. The method of claim 6, further comprising: determining that a norm based on the binary matrix B and the re-computed initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 is 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 (ƒ)}.sub.0.

8. The method of claim 1, wherein the signal ƒ(x) is received from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

9. A system for identifying magnitudes of frequency components of a signal ƒ(x), the system comprising: a first processor; and a first memory in electrical communication with the first processor, the first memory comprising instructions which, when executed by a processing unit comprising at least one of the first processor and a second processor, and in electronic communication with a memory module comprising at least one of the first memory and a second memory, program the processing unit to: obtain at a receiver a first plurality of sets of samples of a signal received by an antenna or a sensor, each set comprising a plurality of samples of the signal ƒ(x) obtained using a respective sampling interval based on a respective shuffling parameter, cardinality of each set being 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, wherein obtaining each set of signal samples comprises subsampling by the receiver the signal at a respective rate less than the Nyquist rate, and at least one set of signal samples determines a support set ; and compute an initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 using a support set and frequency coefficients of the first plurality of sets of signal samples, the support set identifying indices of the non-zero frequency coefficients of the signal ƒ(x).

10. The system of claim 9, wherein to obtain the first plurality of sets of signal samples the processing unit is configured to access from the memory module T sets of signal samples, wherein: T is on the order of O(1) or on the order of O(R log.sub.R N); and a t-th set of signal samples, wherein index t ranges from 1 through T, comprises P.sup.(t) signal samples, P.sup.(t) being greater than R and smaller than N.

11. The system of claim 9, wherein to obtain the first plurality of sets of signal samples the processing unit is configured to: determine the number of samples (N) of the signal ƒ(x) according to the Nyquist rate, and obtain the upper bound (R) of the number of non-zero frequency components of the signal; select a number of sets (T), wherein T is on the order of O(1) or on the order of O(R log.sub.R N), and select T recovery parameters P.sup.(t), wherein index t ranges from 1 through T; for each recovery parameter P.sup.(t), configure a signal sampler to sample the signal ƒ(x) at P.sup.(t) distinct locations, providing the first plurality of T sets of signal samples, each sample set having P.sup.(t) samples.

12. The system of claim 9, wherein to compute the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 the processing unit is programmed to: multiply the first plurality of sets of signal samples by: (i) a discrete Fourier transform (DFT) matrix F comprising T DFT submatrices, each submatrix having a size P.sup.(t), and index t ranging from 1 through T, and (ii) a binary matrix B based on the support set .

13. The system of claim 9, wherein the instructions further program the processing unit to: generate a set of frequency coefficients {circumflex over (ƒ)} using a binary matrix B based on the support set S and the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0.

14. The system of claim 9, wherein the instructions further program the processing unit to: determine that a norm based on a binary matrix B based on the support set and the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 is not less than a selected threshold; obtain a second plurality of sets of signal samples, each set comprising a plurality of samples of the signal ƒ(x); and re-compute the initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 using the support set and the second plurality of sets of signal samples.

15. The system of claim 14, wherein the instructions further program the processing unit to: determine that a norm based on the binary matrix B and the re-computed initial set of frequency coefficients {circumflex over (ƒ)}.sub.0 is less than or equal to the selected threshold; and generate a set of frequency coefficients {circumflex over (ƒ)} using the binary matrix B and the re-computed initial set of frequency coefficients {circumflex over (ƒ)}.sub.0.

16. The system of claim 9, wherein the signal ƒ(x) is received from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

17. A method for identifying frequency components of a signal ƒ(x), the method comprising performing by a processor the steps of: (a) initially designating a current set of candidate support coefficients (.sub.k) as a current set of aliased support coefficients .sub.k; (b) obtaining at a receiver a plurality of sets of samples comprising a first set of K shuffled samples of the signal ƒ(x), received by an antenna or a sensor and corresponding to a first sampling interval based on a first shuffling parameter, wherein K is a fraction of a number of samples (N) of the signal ƒ(x) according to a Nyquist rate, and wherein obtaining the first set of K shuffled samples comprises subsampling by the receiver the signal at a rate less than the Nyquist rate, at least one set of signal samples determining the set of aliased support coefficients .sub.k; (c) filtering the shuffled samples in the first set, and computing a first plurality of frequency coefficients of the shuffled samples in the first set; and (d) removing from the current set of aliased support coefficients .sub.k a 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 plurality of frequency coefficients is less than a threshold.

18. The method of claim 17, further comprising: obtaining a second set of shuffled samples of the signal ƒ(x), corresponding to a second sampling interval based on a second shuffling parameter; filtering the shuffled samples in the second set, and computing a second plurality of frequency coefficients of the shuffled samples in the second set; and removing from the current set of aliased support coefficients .sub.k a subset of candidate support coefficients wherein, for each candidate support coefficient in the subset a value of a corresponding computed frequency coefficient in the second plurality of frequency coefficients is less than the threshold.

19. The method of claim 18, wherein: the first shuffling parameter is a first coprime of an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); and the second shuffling parameter is a second, different coprime of the index limit (M.sub.k).

20. The method of claim 17, wherein: an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k) is associated with a support set growth factor ρ.sub.k that is a ratio of an index limit (M.sub.k) and the number of samples in the first set of shuffled samples (K).

21. The method of claim 17, further comprising: obtaining a prior set of aliased support coefficients .sub.k−1; and dealiasing the prior set of aliased support coefficients to obtain the current set of candidate support coefficients (.sub.k).

22. The method of claim 17, further comprising: (e) selecting a next index limit (M.sub.k+1) of a next set of candidate support coefficients (.sub.k+1), wherein the next index limit is greater than a current index limit (M.sub.k) of the current set of candidate support coefficients(.sub.k); (f) after the step (d), dealiasing the current set of aliased support coefficients (.sub.k) using the next index limit (M.sub.k+1), to obtain the next set of candidate support coefficients (.sub.k+1); (g) updating the current set of candidate support coefficients (.sub.k) by designating the next set of candidate support coefficients (.sub.k+1) as the current set of candidate support coefficients (M.sub.k); and (h) repeating the steps (a) through (d) using the updated current set of candidate support coefficients (.sub.k).

23. The method of claim 22, further comprising: determining that a next index limit (M.sub.k+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 coefficients .sub.k as a final set of aliased support coefficients .

24. The method of claim 17, wherein obtaining the first set of shuffled samples of the signal ƒ(x) comprises: selecting the first sampling interval using the first shuffling parameter (Q.sub.k.sup.(l)), wherein the first shuffling parameter corresponds to an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); 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.

25. The method of claim 17, wherein obtaining the first set of shuffled samples of the signal ƒ(x) comprises: selecting the first shuffling parameter (Q.sub.k.sup.(l)) corresponding to an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); and accessing from memory the first set of shuffled samples corresponding to the first shuffling parameter.

26. The method of claim 17, wherein filtering the shuffled samples in the first set comprises filtering the first set using a Gaussian filter having a standard deviation (σ) 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.

27. The method of claim 17, further comprising generating the current set of candidate support coefficients by: selecting a first index limit (M.sub.1), wherein M.sub.1 is less than the number of samples (N) of the signal ƒ(x) according to the Nyquist rate; obtaining a plurality of samples of the signal ƒ(x); determining M.sub.1 frequency coefficients of a signal comprising the plurality of samples; and including each non-zero frequency coefficient from the M.sub.1 frequency coefficients as a respective candidate support coefficient in the current set of candidate support coefficients.

28. The method of claim 27, wherein: a non-zero frequency coefficient comprises a frequency coefficient having a magnitude greater than a specified threshold value; and a frequency coefficient designated as a zero coefficient has a magnitude less than or equal to the specified threshold value.

29. A system for identifying frequency components of a signal ƒ(x), the system comprising: a first processor; and a first memory in electrical communication with the first processor, the first memory comprising instructions which, when executed by a processing unit comprising at least one of the first processor and a second processor, and in electronic communication with a memory module comprising at least one of the first memory and a second memory, program the processing unit to: (a) initially designate a current set of candidate support coefficients (.sub.k) as a current set of aliased support coefficients .sub.k; (b) obtain a plurality of sets of samples comprising a first set of K shuffled samples of the signal ƒ(x), received by an antenna or a sensor and corresponding to a first sampling interval based on a first shuffling parameter, wherein K is a fraction of a number of samples (N) of the signal ƒ(x) according to a Nyquist rate, wherein obtaining the first set of K shuffled samples comprises subsampling by the receiver the signal at a rate less than the Nyquist rate, at least one set of signal samples determining the set of aliased support coefficients .sub.k; (c) filter the shuffled samples in the first set, and computing a first plurality of frequency coefficients of the shuffled samples in the first set; and (d) remove from the current set of aliased support coefficients .sub.k a 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 plurality of frequency coefficients is less than a threshold.

30. The system of claim 29, wherein the processing unit is further programmed to: obtain a second set of shuffled samples of the signal ƒ(x), corresponding to a second sampling interval based on a second shuffling parameter; filter the shuffled samples in the second set, and computing a second plurality of frequency coefficients of the shuffled samples in the second set; and remove from the current set of aliased support coefficients .sub.k a subset of candidate support coefficients wherein, for each candidate support coefficient in the subset a value of a corresponding computed frequency coefficient in the second plurality of frequency coefficients is less than the threshold.

31. The system of claim 30, wherein: the first shuffling parameter is a first coprime of an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); and the second shuffling parameter is a second, different coprime of the index limit (M.sub.k).

32. The system of claim 29, wherein: an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k) is associated with a support set growth factor ρ.sub.k that is a ratio of an index limit (M.sub.k) and the number of samples in the first set of shuffled samples (K).

33. The system of claim 29, wherein the processing unit is further programmed to: obtain a prior set of aliased support coefficients .sub.k−1; and dealias the prior set of aliased support coefficients to obtain the current set of candidate support coefficients (.sub.k).

34. The system of claim 29, wherein the processing unit is further programmed to: (e) select a next index limit (M.sub.k+1) of a next set of candidate support coefficients (.sub.k+1), wherein the next index limit is greater than a current index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); (f) after performing the operation (d), dealias the current set of aliased support coefficients (.sub.k) using the next index limit (M.sub.k+1), to obtain the next set of candidate support coefficients (.sub.k+1); (g) update the current set of candidate support coefficients (.sub.k) by designating the next set of candidate support coefficients (.sub.k+1) as the current set of candidate support coefficients (.sub.k); and (h) repeat the operations (a) through (d) using the updated current set of candidate support coefficients (.sub.k).

35. The system of claim 34, wherein the processing unit is further programmed to: determine that a next index limit (M.sub.k+1) is not less than a number of samples (N) of the signal ƒ(x) according to a Nyquist rate; and after the operation (h), designate the current set of aliased support coefficients .sub.k as a final set of aliased support coefficients .

36. The system of claim 29, wherein to obtain the first set of shuffled samples of the signal ƒ(x) the processing unit is programmed to: select the first sampling interval using the first shuffling parameter (Q.sub.k.sup.(l)), wherein the first shuffling parameter corresponds to an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); configure a signal sampler to sample the signal ƒ(x) using the first sampling interval, to obtain a first sampled signal; and shuffle the first sampled signal using the first shuffling parameter and the index limit.

37. The system of claim 29, wherein to obtain the first set of shuffled samples of the signal ƒ(x) the processing unit is programmed to: select the first shuffling parameter (Q.sub.k.sup.(l)) corresponding to an index limit (M.sub.k) of the current set of candidate support coefficients (.sub.k); and access from the memory module the first set of shuffled samples corresponding to the first shuffling parameter.

38. The system of claim 29, wherein to filter the shuffled samples in the first set the processing unit is programmed to filter the first set using a Gaussian function having a standard deviation (σ) 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.

39. The system of claim 29, wherein the processing unit is further programmed to generate the current set of candidate support coefficients by: selecting a first index limit (M.sub.1), wherein M.sub.1 is less than the number of samples (N) of the signal ƒ(x) according to the Nyquist rate; obtaining a plurality of samples of the signal ƒ(x); determining M.sub.1 frequency coefficients of a signal comprising the plurality of samples; and including each non-zero frequency coefficient from the M.sub.1 frequency coefficients as a respective candidate support coefficient in the current set of candidate support coefficients.

40. The system of claim 39, wherein: a non-zero frequency coefficient comprises a frequency coefficient having a magnitude greater than a specified threshold value; and a frequency coefficient designated as a zero coefficient has a magnitude less than or equal to the specified threshold value.

41. A method for generating non-zero frequency coefficients of a signal ƒ(x), the method comprising performing by a processor the steps of: obtaining at a receiver a plurality of sets of samples of the signal ƒ(x) received by an antenna or a sensor and comprising a first set corresponding to a first sampling interval and a different second set, wherein obtaining each of the sets of samples comprises subsampling by the receiver the signal at a rate less than a Nyquist rate corresponding to the signal; determining indices of the non-zero frequency coefficients of the signal ƒ(x) using at least the first set, the first sampling interval corresponding to the indices of non-zero frequency coefficients; and determining values of the non-zero frequency coefficients of the signal ƒ(x) using the indices thereof and at least the second set.

42. The method of claim 41, wherein determining the indices of the non-zero frequency coefficients comprises: 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); selecting a sample size K, a number of iterations P, a plurality of support set growth factors ρ.sub.k, and a plurality of index limits M.sub.k, each index limit representing an index limit of a respective set of candidate support coefficients .sub.k, wherein M.sub.k=ρ.sub.kK and $N=K\underset{k=1}{\overset{P}{.Math.}}{\rho }_k;$ and determining a current set of aliased support coefficients .sub.k using the first set of samples of the signal ƒ(x), the first set comprising less than N samples.

43. The method of claim 41, wherein determining the indices of the non-zero frequency coefficients comprises: performing at least one iteration, comprising: obtaining a current set of candidate support coefficients .sub.k by dealiasing the current set of aliased support coefficients .sub.k; determining a next set of aliased support coefficients using the current set of candidate support coefficients .sub.k and at least one set of samples from the plurality of sets of samples of the signal ƒ(x), each one of the at least one set of samples being obtained using a sampling interval based on M.sub.k; and designating the next set of aliased support coefficients as the current set of aliased support coefficients .sub.k.

44. The method of claim 41, wherein determining the indices of the non-zero frequency coefficients comprises: performing at least one iteration, comprising: obtaining a current set of candidate support coefficients .sub.k by dealiasing the current set of aliased support coefficients .sub.k; determining a next set of aliased support coefficients using the current set of candidate support coefficients .sub.k and at least one set of samples from the plurality of sets of samples of the signal ƒ(x), each one of the at least one set of samples being obtained using a sampling interval based on M.sub.k; and designating the next set of aliased support coefficients as final set of aliased support coefficients S, the final set of aliased support coefficients representing the indices of the non-zero frequency coefficients.

45. The method of claim 41, wherein the signal ƒ(x) comprises a one-dimensional signal obtained from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

46. The method of claim 41, further comprising: 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), wherein obtaining the plurality of sets of samples of the signal ƒ(x) comprises: selecting a plurality of sample indices of ƒ(x); and obtaining samples of y(x) at a plurality of sample indices of y(x), each one being determined via the inverse mapping and a respective selected sample index of ƒ(x).

47. The method of claim 46, wherein the non-zero frequency coefficients of the signal ƒ(x) represent non-zero frequency coefficients of the multi-dimensional signal y(x).

48. The method of claim 46, wherein the multi-dimensional signal y(x) comprises a multi-dimensional signal obtained from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

49. The method of claim 41, wherein: a non-zero frequency coefficient comprises a frequency coefficient having a magnitude greater than a specified threshold value; and a frequency coefficient designated as a zero coefficient has a magnitude less than or equal to the specified threshold value.

50. The method of claim 41, wherein: the step of determining indices of the non-zero frequency coefficients is performed using a first value (N.sub.1) of a number of samples of the signal ƒ(x) according to a Nyquist rate, to obtain a first support set .sup.1; the step of determining indices of the non-zero frequency coefficients is repeated using a second value (N.sub.2) of the number of samples, to obtain a second support set .sup.2, the method further comprising: 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.

51. A system for generating non-zero frequency coefficients of a signal ƒ(x), the system comprising: a first processor; and a first memory in electrical communication with the first processor, the first memory comprising instructions which, when executed by a processing unit comprising at least one of the first processor and a second processor, and in electronic communication with a memory module comprising at least one of the first memory and a second memory, program the processing unit to: obtain at a receiver a plurality of sets of samples of the signal ƒ(x) received by an antenna or a sensor and comprising a first set corresponding to a first sampling interval and a different second set, wherein obtaining each of the sets of samples comprises subsampling by the receiver the signal at a rate less than a Nyquist rate corresponding to the signal; determine indices of the non-zero frequency coefficients of the signal ƒ(x) using at least the first set, the first sampling interval corresponding to the indices of non-zero frequency coefficients; and determine values of the non-zero frequency coefficients of the signal ƒ(x) using the indices thereof and at least the second set.

52. The system of claim 51, wherein to determine the indices of the non-zero frequency coefficients the processing unit is programmed to: obtain: (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); select a sample size K, a number of iterations P, a plurality of support set growth factors ρ.sub.k, and a plurality of index limits M.sub.k, each index limit representing an index limit of a respective set of candidate support coefficients .sub.k, wherein M.sub.k=ρ.sub.kK and $N=K\underset{k=1}{\overset{P}{.Math.}}{\rho }_k;$ and determine a current set of aliased support coefficients .sub.k using the first set of samples of the signal ƒ(x), the first set comprising less than N samples.

53. The system of claim 51, wherein to determine the indices of the non-zero frequency coefficients the processing unit is programmed to: perform at least one iteration, comprising: obtaining a current set of candidate support coefficients .sub.k by dealiasing the current set of aliased support coefficients .sub.k; determining a next set of aliased support coefficients using the current set of candidate support coefficients .sub.k and at least one set of samples from the plurality of sets of samples of the signal ƒ(x), each one of the at least one set of samples being obtained using a sampling interval based on M.sub.k; and designating the next set of aliased support coefficients as the current set of aliased support coefficients .sub.k.

54. The system of claim 51, wherein to determine the indices of the non-zero frequency coefficients the processing unit is programmed to: perform at least one iteration, comprising: obtaining a current set of candidate support coefficients .sub.k by dealiasing the current set of aliased support coefficients .sub.k; determining a next set of aliased support coefficients using the current set of candidate support coefficients .sub.k and at least one set of samples from the plurality of sets of samples of the signal ƒ(x), each one of the at least one set of samples being obtained using a sampling interval based on M.sub.k; and designating the next set of aliased support coefficients as final set of aliased support coefficients , the final set of aliased support coefficients representing the indices of the non-zero frequency coefficients.

55. The system of claim 51, wherein the signal ƒ(x) comprises a one-dimensional signal obtained from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

56. The system of claim 51, wherein the processing unit is further programmed to: map a multi-dimensional signal y(x) into a one-dimensional signal ƒ(x); and store an inverse mapping identifying a sample index for y(x) corresponding to a sample index for ƒ(x), wherein to obtain the plurality of sets of samples of the signal ƒ(x) the processing unit is programmed to: select a plurality of sample indices of ƒ(x); and obtain samples of y(x) at a plurality of sample indices of y(x), each one being determined via the inverse mapping and a respective selected sample index of ƒ(x).

57. The system of claim 56, wherein the non-zero frequency coefficients of the signal ƒ(x) represent non-zero frequency coefficients of the multi-dimensional signal y(x).

58. The system of claim 56, wherein the multi-dimensional signal y(x) comprises a multi-dimensional signal obtained from a signal source comprising at least one of a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, and an acoustic sensor.

59. The system of claim 51, wherein: a non-zero frequency coefficient comprises a frequency coefficient having a magnitude greater than a specified threshold value; and a frequency coefficient designated as a zero coefficient has a magnitude less than or equal to the specified threshold value.

60. The system of claim 51, wherein the processing unit is programmed to: determine the indices of the non-zero frequency coefficients using a first value (N.sub.1) of a number of samples of the signal ƒ(x) according to a Nyquist rate, to obtain a first support set .sup.1; determine the indices of the non-zero frequency coefficients again, using a second value (N.sub.2) of the number of samples, to obtain a second support set .sup.2; generate a final support set as a union of the first and second support sets; and designate coefficients of the final support set as the indices of the non-zero frequency coefficients.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) In the following description, various embodiments of the present invention are described with reference to the following drawings, in which:

(2) FIG. 1 is a table showing the runtime and sampling complexities of various known sFFT techniques and the last row shows the runtime and sampling complexities of a system according to one embodiment of our technique;

(3) FIG. 2 depicts an overall process of determining the locations of the non-zero frequency components of a signal, and determining their values, according to one embodiment;

(4) FIG. 3 shows an example of determining the locations of the non-zero frequency components of a signal, according to one embodiment;

(5) FIG. 4 depicts a process of iteratively determining the locations of the non-zero frequency components of a signal, according to one embodiment;

(6) FIG. 5 depicts a subprocess used by the process of FIG. 4, according to one embodiment;

(7) FIG. 6 shows an example of determining the aliased locations of the non-zero frequency components of a signal, according to one embodiment of the subprocess shown in FIG. 5;

(8) FIG. 7 depicts a process of determining the values of the non-zero frequency components of a signal when their locations are known, according to one embodiment;

(9) FIGS. 8A and 8B illustrate a mapping of a multi-dimensional signal into a one-dimensional signal, according to one embodiment;

(10) FIG. 9A schematically depicts a system for determining the locations of the non-zero frequency components of a signal, determining their values, and collecting signal samples that are required, according to one embodiment;

(11) FIG. 9B schematically depicts a subsystem for collecting samples of a signal, where the samples are used for determining the locations and/or values of the non-zero frequency components of the signal, according to one embodiment;

(12) FIG. 10 is a table of various parameters used for experimental validation of a system according to one embodiment; and

(13) FIGS. 11A and 11B show comparison results for the performance of an embodiment of a system described herein with that of a known system.

DETAILED DESCRIPTION

(14) Statement of the Problem and Preliminaries We introduce the notation used throughout the remainder of the discussion. Unless otherwise stated, we consider a one-dimension (1D) function or signal ƒ(x) of the form:

(15) $\begin{array}{cc}f\left(x\right)=\underset{j=0}{\overset{N-1}{.Math.}}{e}^{-2\pi \mathrm{ixj}}\left({\hat{f}}_j+\eta {\hat{v}}_j\right),& \left(3\right)\end{array}$
for some finite 0<N∈ and noise level 0≤√{square root over (N)}∥{circumflex over (ν)}∥.sub.2≤η. It is further assumed that the vector {circumflex over (ƒ)} has real nonnegative elements, i.e., {circumflex over (ƒ)}.sub.j≥0, ∀j, and that its support:
:={j∈{0,1, . . . ,N−1}:|{circumflex over (ƒ)}.sub.j|≠0}
satisfies 0≤#≤R<N<∞, where # indicates cardinality. In particular, we are interested in the case where R«N. Given some accuracy parameter ϵ above the noise level η, the problem involves computing an R-sparse vector ƒ* such that:
∥ƒ*−{circumflex over (ƒ)}∥.sub.2≤ε∥{circumflex over (ƒ)}∥.sub.2
We shall denote by the Fourier transform (and * its inverse/adjoint), i.e.,
[{circumflex over (ƒ)}(ξ)](x)=∫.sub.de.sup.−2πix.Math.ξ{circumflex over (ƒ)}(ξ)dξ
where d represents the ambient dimension. The size-N Discrete Fourier Transform (DFT) is defined as:

(16) $\begin{array}{cc}{f}_{n_{;}N}={\left[F_N\hat{f}\right]}_n=\underset{j=0}{\overset{N-1}{.Math.}}{e}^{-2\pi i\frac{nj}{N}}{\hat{f}}_j,n=0,1,\mathrm{.Math.},N-1.& \left(4\right)\end{array}$
A Sparse FFT in 1D

(17) We describe a fast way to compute the one-dimensional DFT of a bandlimited and periodic function ƒ(x) of the form of Eq. (3). Our approach to solving this problem can be broken into two separate steps: in the first step, the support of the vector {circumflex over (ƒ)} is recovered, and in the second step, the nonzero values of {circumflex over (ƒ)} are computed using the knowledge of the recovered support. We describe the algorithm in the noiseless case, followed by a discussion of its stability to noise. Steps of the signal sampling and frequency-coefficient computation processes are provided in Algorithms 1-4.

(18) Finding the Support

(19) Referring to the example in FIG. 3, from a high-level perspective, our support-finding scheme uses three major ingredients: 1) sampling or sub-sampling, 2) shuffling, and 3) low-pass filtering. Sampling or Sub-sampling reduces the size of the problem to a manageable level, but leads to aliasing. Nonetheless, when the nonnegativity assumption is satisfied, an aliased Fourier coefficient is nonzero if and only if its corresponding aliased lattice contains an element of the true support (the positivity condition avoids cancellation). This provides a useful criterion to discriminate between frequency elements that belong to the support and frequency elements that do not.

(20) To illustrate, we proceed through an example and refer to FIG. 3. To begin with, consider k, N, M.sub.k∈, 0<α<1 and .sub.k, .sub.k, .sub.k⊂{0, 1, . . . , N−1}. We define the following: the aliased support .sub.k at iteration/step k corresponds to the indices of the elements of the true support modulo M.sub.k; the working support at iteration/step k corresponds to the set .sub.k:={0, 1, . . . , M.sub.k−1}; a candidate support at iteration/step k is any set satisfying .sub.k⊂⊂.sub.k of size (ρR log(R)). Line 0 (FIG. 3) represents a lattice (thin tickmarks) of size,

(21) $N=40=5\underset{i=1}{\overset{3}{.Math.}}2=K\underset{i=1}{\overset{P}{.Math.}}{\rho }_i$
which contains only 3 positive frequencies (black dots; ={1, 23, 35}). In the beginning, (step k=0) only the fact that ⊂{0, 1, . . . , N−1} is known.

(22) The first step (k=1) is performed as follows: Letting:

(23) $M_1=\frac{N}{\underset{i=2}{\overset{P}{.Math.}}{\rho }_i}={\rho }_{1}K=\mathrm{��}\left(R\log \left(R\right)\right)$
sample the function ƒ(x) at,

(24) $x_{n_1\underset{i=2}{\overset{P}{.Math.}}{\rho }_i;N}=\frac{n_1\underset{i=2}{\overset{P}{.Math.}}{\rho }_i}{N}=\frac{n_1}{M_1}=x_{n_1;M_1}$
to obtain,

(25) $\begin{array}{cc}{f}_{n_1{\Pi }_{i=2}^P\rho i;N}=\underset{j=0}{\overset{N-1}{.Math.}}{e}^{-2\pi i\frac{n_1{\Pi }_{i=2}^P{\rho }_{i}j}{N}}\widehat{f_j}=\underset{l=0}{\overset{M_1-1}{.Math.}}{e}^{-2\pi i\frac{n_{1}l}{M_1}}\left(\underset{j:{\mathrm{jmodM}}_1=l}{.Math.}{\hat{f}}_j\right)=\underset{l=0}{\overset{M_1-1}{.Math.}}{e}^{-2\pi i\frac{n_{1}l}{M_1}}{\hat{f}}_j^{\left(1\right)}=f_{n_1;M_1}& \left(5\right)\end{array}$
for n.sub.1∈.sub.1:={0, 1, . . . , M.sub.1−1} defined as the candidate support in the first step.

(26) The samples correspond to a DFT of size M.sub.1 of the vector {circumflex over (ƒ)}.sup.(1) with entries that are an aliased version of those of the original vector {circumflex over (ƒ)}. These can be computed through the FFT in order (M.sub.1 log(M.sub.1))=(R log.sup.2(R)). In this first step, it is further possible to rapidly identify the aliased support .sub.1 from the knowledge of {circumflex over (ƒ)}.sup.(1) since the former correspond to the set,
{l∈{0,1, . . . ,M.sub.1−1}:{circumflex over (ƒ)}.sub.l.sup.(1)≠0}
due to the fact that
{circumflex over (ƒ)}.sub.l.sup.(1):=Σ.sub.j:jmodM.sub.1.sub.=l{circumflex over (ƒ)}.sub.j>0⇔l∈.sub.1
following the nonnegativity assumption. In our example, M.sub.1=ρ.sub.1K=2.Math.5=10 which leads to
.sub.1={1 mod 10,23 mod 10,35 mod 10}={1,3,5}={l∈{0,1, . . . ,9}:{circumflex over (ƒ)}.sub.l.sup.(1)≠0}
.sub.1=.sub.1={0,1, . . . ,9}.
This is shown on line 1 of FIG. 3. For this first step, the working support .sub.1 is equal to the candidate support .sub.1.

(27) Then, proceed to the next step (k=2) as follows: let,

(28) 0$M_2={\rho }_2{M}_1=K\underset{i=1}{\overset{2}{.Math.}}{\rho }_i=5.Math.2^2=20$
and consider the samples,

(29) $f_{n_2{\Pi }_{i=3}^P{\rho }_i;N}=\underset{l=0}{\overset{M_2-1}{.Math.}}{e}^{-2\pi i\frac{n_{2}l}{M_2}}\left(\underset{j:\mathrm{jmo}d{M}_2=l}{.Math.}\widehat{f_j}\right)=\underset{l=0}{\overset{M_2-1}{.Math.}}{e}^{-2\pi i\frac{n_{2}l}{M_2}}{\widehat{f_l}}^{\left(2\right)}=f_{n_2;M_2}$
for n.sub.2=0, 1, . . . , M.sub.2−1 as before. Here however, knowledge of .sub.1 is incorporated. Indeed, since M.sub.2 is a multiple of M.sub.1, it follows upon close examination that,
.sub.2⊂∪.sub.k=0.sup.ρ.sup.1.sup.−1(.sub.1+kM.sub.1):=.sub.2.

(30) That is, the set .sub.2, defined as the union of ρ.sub.1=O(1) translated copies of .sub.1, must itself contain .sub.2. Furthermore, it is of size (ρ.sub.1#.sub.1)=(ρR log(R)) by construction. It is thus a proper candidate support (by definition). In our example, one obtains
∪.sub.k=0.sup.1(.sub.1+kM.sub.1)={1,3,5}∪{1+10,3+10,5+10}={1,3,5,11,13,15}=.sub.2,
which contains the aliased support,
.sub.2={1 mod 20,23 mod 20,35 mod 20}={1,3,15}
as shown on line 2 of FIG. 3.

(31) The working support becomes .sub.2:={0, 1, . . . , 19}. Once again, it is possible to recover .sub.2 by leveraging the fact that {l∈{0, 1, . . . , M.sub.2−1}:{circumflex over (ƒ)}.sub.l.sup.(2)≠0}=.sub.2. Here however, the cost is higher since computing {circumflex over (ƒ)}.sup.(2) involves performing an FFT of size M.sub.2=20. Continuing in the fashion of the first step, the cost would increase exponentially with k, so different steps are required to contain the cost. Such steps involve a special kind of shuffling and filtering of the samples followed by an FFT, and we describe this in detail below. Altogether, it is shown that .sub.k can now be recovered from the knowledge of .sub.k at any step k using merely (R log(R)) samples and (R log.sup.2(R)) computational steps.

(32) Following the rapid recovery of .sub.2, we proceed in a similar fashion until .sub.k:={0, 1, . . . , N−1} at which point .sub.k=. Throughout this process, the size of the aliased support .sub.k and candidate support .sub.k remain of order (R log(R)) while the size of the working support increases exponentially fast; i.e.,

(33) $\#k_{}=\mathrm{��}\left(K\underset{i=1}{\overset{k}{.Math.}}{\rho }_i\right)\ge 2^k.Math.R.$
This therefore implies

(34) $P=\mathrm{��}\left(\log \left(\frac{N}{R}\right)\right)$
“dealiasing” steps, and thus a total cost of

(35) $\widetilde{\mathrm{��}}\left(R\log \left(R\right)\log \left(\frac{N}{R}\right)\right)$
samples and

(36) $\widetilde{\mathrm{��}}\left(R{\log }^2\left(R\right)\log \left(\frac{N}{R}\right)\right)$
computational steps to identify .

(37) Thus, with reference to FIG. 3, the support is computed as follows. Line 0: Initialization; (unknown) elements of correspond to black dots and lie in the grid {0, 1, . . . , N−1}. Line 1: First step; elements of the candidate support .sub.1 are represented by thin tickmarks and those of the aliased support .sub.1 by thick tickmarks. .sub.1 is a subset of .sub.1 and both lie in the working support {0, 1, . . . , M.sub.1−1}. Line 2: Second step; elements of the candidate support .sub.2 correspond to thin tickmarks and are obtained through de-aliasing of .sub.1. Elements of the aliased support .sub.2 correspond to thick tickmarks. Both lie in the working support {0, 1, . . . , M.sub.2−1}. M.sub.2 is a constant factor/multiple of M.sub.1. Line 3: The final step correspond to the step when the working set is equal to {0, 1, . . . , N−1}. While FIG. 3 shows three iterations/steps, in general

(38) $P=\mathrm{��}\left(\log \left(\frac{N}{R}\right)\right)$
iterations may be performed, where P can be any number such as 1, 5, 8, 10, 30, etc.

(39) The steps of this support-recovery algorithm are described in Algorithm 2, the correctness of which is guaranteed by the following proposition.

(40) Proposition 1. In the noiseless case, Algorithm 2 outputs , the support of the nonnegative R-sparse vector {circumflex over (ƒ)} with probability at least (1−p) using

(41) $\widetilde{\mathrm{��}}\left(R\log \left(R\right)\log \left(\frac{N}{R}\right)\right)$
samples and

(42) $\widetilde{\mathrm{��}}\left(R{\log }^2\left(R\right)\log \left(\frac{N}{R}\right)\right)$
computational steps.

(43) Proof: Refer to Algorithm 2 and Proposition 5 as well as the above discussion.

(44) From the knowledge of , it is possible to recover the actual values of {circumflex over (ƒ)} rapidly and with few samples. This is the second major step of the sMFFT which we describe below.

(45) Rapid Recovery of .sub.k from Knowledge of .sub.k

(46) Details are given here as to how to solve the problem of rapidly recovering the aliased support .sub.k from the knowledge of a candidate support .sub.k. Before proceeding, a few definitions are introduced.

(47) Definition 1 Let 1≤K≤M∈. Then, define the set (K; M) as,

(48) $\left(K;M\right):=\left\{m\in \left\{0,1,\mathrm{.Math.},M-1\right\}:m\le \frac{K}{2}\mathrm{or}.Math.m-M.Math.<\frac{K}{2}\right\}$
Definition 2 Let 0<M∈. Then, we define the set Q(M) as,
Q(M):={q∈[0,M)∩:q⊥M},
where the symbol ⊥ between two integers indicates they are coprime.

(49) Algorithm 3 (FIG. 5) shows how to solve the aliased support recovery problem rapidly; correctness is guaranteed by Proposition 2, which relies on Proposition 5. Proposition 5 states that if the elements of an aliased vector of size M.sub.k with aliased support .sub.k containing at most R nonzeros are shuffled (according to appropriate random permutation) and subsequently convoluted with a (low-frequency) Gaussian (or another low pass filter function), then the probability that the resulting value at a location m∈.sub.k.sup.c is of order (1) is small. If m∈.sub.k the value at m is of order O(1) with probability 1.

(50) This realization allows us to develop an efficient statistical test to identify .sub.k from the knowledge of .sub.k. The process is shown schematically in FIG. 6. Specifically, the following four steps are performed: 1) permute samples randomly, 2) apply a diagonal Gaussian filter, 3) compute a small FFT, 4) eliminate elements that do not belong to the aliased support. Instead of a Gaussian filter, any other low-pass filter can be applied.

(51) To illustrate, we once again proceed through an example shown in FIG. 6. To begin with, assume M.sub.k=40, and
.sub.k={1,23,35},.sub.k={1,3,15,21,23,35},.sub.k={1,2, . . . ,39}
as in step/iteration k=3 shown in Algorithm 2. This alias support is shown at line A of FIG. 6.

(52) The first step is to randomly shuffle the elements of .sub.k within .sub.k by applying a permutation operator Π.sub.Q(⋅) in the sample space:

(53) 0$\begin{array}{cc}\underset{Q}{.Math.}\left(f_{n;M_k}\right)=f_{\left({n\left[Q\right]}_{M_k}^{-1}\right)mod{M}_k;M_k}& \left(6\right)\end{array}$
for some integer Q∈Q(M.sub.k) ([Q].sub.M.sub.k.sup.−1 being the unique inverse of Q modulo M.sub.k) by Lemma 3. This is equivalent to shuffling in frequency space as: j.fwdarw.(jQ) mod M.sub.k. Indeed, after shuffling, evaluating at (jQ) mod M.sub.k gives,

(54) ${\hat{f}}_{{\left(\left(\mathrm{jQ}\right)\mathrm{mo}d{M}_k\right)\left[Q\right]}_{M_k}^{-1})mod{M}_k}^{\left(k\right)}={\widehat{f_j}}^{\left(k\right)}.$
Furthermore, Lemma 4 shows that if Q is chosen uniformly at random within Q(M.sub.k), the mapped elements of the candidate support .sub.k will be more or less uniformly distributed within the working support .sub.k,

(55) $(.Math.\left(\mathrm{iQ}\right)\mathrm{mod}{M}_k-\left(\mathrm{jQ}\right)\mathrm{mod}{M}_k.Math.\le C.Math.i\ne j)\le \mathrm{��}\left(\frac{C}{M_k}\right).$

(56) This step is followed by the application of a diagonal Gaussian filtering operator Ψ.sub.σ(⋅) having elements:

(57) $\begin{array}{cc}{g}_{\sigma }\left(\frac{m}{M_k}\right)=\sqrt{\pi }\sigma \underset{h\in \mathbb{Z}}{.Math.}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{m+hM}{M}\right)}^2}& \left(7\right)\end{array}$
in the sample space (step 2). Instead of a Gaussian filter, any other low-pass filter can be applied. By the properties of the Fourier transform, this is equivalent to a convolution in the frequency space (line C, FIG. 6), implying the equality:

(58) $\begin{array}{cc}\left[{\Psi }_{\sigma }\left({\Pi }_Q\left(f_{n;M_k}\right)\right)\right]\left(\xi \right)=ℱ*\left[\underset{j\in k_{}}{.Math.}{\widehat{f_j}}^{\left(k\right)}{e}^{-\frac{{.Math.x-\left(\mathrm{jQ}\right){\mathrm{modM}}_k.Math.}^2}{{\sigma }^2}}\right]\left(\xi \right).& \left(8\right)\end{array}$

(59) The function is now bandlimited (with bandwidth of order O(K) due to our choice for a, the standard deviation of the Gaussian filter (or a low-pass filter, in general) that is applied, as shown in Algorithm 3, so this expression can be discretized (samples x, FIG. 6) to produce our main expression,

(60) $\begin{array}{cc}{\varphi }_n^{\left(k\right)}\left(Q\right)={F_{\mathrm{��}\left(K,M_k\right)}\left[{\Psi }_{\sigma }\left({\Pi }_Q\left(f_{n;M_k}\right)\right)\right]}_n=\frac{1}{M_k}\underset{m\in \left(K,M_k\right)}{.Math.}{e}^{2\pi i\frac{nm}{K}}{g}_{\sigma }\left(\frac{m}{M_k}\right)f_{\left({m\left[Q\right]}_{M_k}^{-1}\right){\mathrm{modM}}_k;M_k}& \left(9\right)\end{array}$
In particular, we note that if n if of the form

(61) $j\frac{M_k}{K}$
for j=0, . . . , K−1, the last step can be performed through a small size-K FFT. This corresponds to step 3 of the aliased support recovery algorithm. The knowledge of,

(62) ${\left\{{\varphi }_{j\frac{M_k}{K}}^{\left(k\right)}\right\}}_{j=0}^{K-1}$
can be used to recover .sub.k from .sub.k rapidly, which can be seen as follows.

(63) By construction, ϕ.sub.n.sup.(k) can be “large” only if the distance between n and some element belonging to the shuffled and aliased support, i.e., some element of {(lQ) mod M.sub.k}.sub.l∈.sub.k, is smaller than O(σ), which in turn occurs only if the distance between

(64) $\left[n\frac{K}{M_k}\right]\frac{M_k}{K}$
and the shuffled elements of the aliased support is smaller than O(σ), as well (by the triangle inequality). However, because of the randomness introduced by the shuffling, and because of the particular choice of σ, Proposition 5 shows that for any fixed n∈.sub.k, the probability that a computed element

(65) ${\varphi }_{\left[\left({\mathrm{nQmodM}}_k\right)\frac{K}{M_k}\right]\frac{M_k}{K}}^{\left(k\right)}$
is “large” for multiple independent trials is small if n∈.sub.k∩.sub.k.sup.c and equal to 1 if n∈.sub.k∩.sub.k. This fact allows for the construction of an efficient statistical test based on the knowledge of the quantities found in Eq. (9) to discriminate between the points of .sub.k∩.sub.k.sup.c and those of .sub.k∩.sub.k (step 4). Such a test forms the core of Algorithm 3, and its correctness follows from the Proposition 2.

(66) Proposition 2. In the noiseless case, Algorithm 3 outputs .sub.k, the aliased support of the vector {circumflex over (ƒ)} at step k, with probability at least (1−p) using (log(p)R log(R)) samples and (log(p)R log.sup.2(R)) computational steps.

(67) Proof. Refer to Algorithm 3 and Proposition 5 as well as the above discussion.

(68) Thus, according to FIG. 6, finding the aliased support .sub.k from knowledge of .sub.k involves receiving .sub.k (line A)). First, indices are shuffled in sample space leading to a shuffling in frequency space (line B)). A Gaussian or another low pass filter is applied followed by a small FFT (line C)) on a grid G(x). The points of .sub.k for which the value of the result of the last step at their closest neighbor in G is small are discarded, leaving only the aliased support .sub.k.

(69) As for the computational cost, the permutation and filtering (multiplication) steps (1 and 2) both incur a cost of (R log(R)) since only the samples for which the filter is of order O(1) are considered (and there are (K)=(R log(R)) of them following our choice of σ and K). These are followed by an FFT (step 3) of size O(K) which carries a cost of order (R log.sup.2(R)). Finally, step 4 involves checking a simple property on each of the M.sub.k elements of .sub.k incurring a cost of (R log(R)). This is repeated (log(p)) times for a probability (1−p) of success. Thus, extracting .sub.k from .sub.k requires merely (log(p)R log(R)) samples and (log(p)R log.sup.2(R)) computational time for fixed p.

(70) Recovering Values from Knowledge of the Support

(71) Assume a set size O(K) containing the support S has been recovered. We now show how the values of the nonzero Fourier coefficients of {circumflex over (ƒ)} in Eq. (3) can be rapidly computed using this information. For this purpose, assume ƒ(x) can be sampled at locations:

(72) 0${\left\{\frac{q\mathrm{mod}{P}^{\left(t\right)}}{P^{\left(t\right)}}\right\}}_{q=0}^{P^{\left(t\right)}-1}$
for t=0, 1, . . . , T, and {P.sup.(t)}.sub.t=1.sup.T are some random prime or co-prime numbers on the order of (R log.sub.R(N)) (see Algorithm 4, shown in FIG. 7). It follows that

(73) $\begin{array}{cc}{f}_{q\mathrm{mod}{P}^{\left(t\right)};p\left(t\right)}^{\left(t\right)}=\underset{j\in }{\overset{}{.Math.}}{e}^{-2\pi i\frac{q\left(j{\mathrm{modP}}^{\left(t\right)}\right)}{P^{\left(t\right)}}}\widehat{f_j}=\underset{l=0}{\overset{P^{\left(t\right)}-1}{.Math.}}{e}^{-2\pi i\frac{ql}{P^{\left(t\right)}}}\left(\underset{j\in :j\mathrm{mod}{P}^{\left(t\right)}=l}{.Math.}\widehat{f_j}\right)& \left(10\right)\end{array}$
for t=0, 1, . . . , T. The outer sum is seen to be a DFT of size P.sup.(t) of a shuffled and aliased vector, whereas the inner sum can be expressed as the application of a binary matrix B.sub.q,j.sup.(t) with entries

(74) $B_{q,j}^{\left(t\right)}=\{\begin{array}{ccc}1& \mathrm{if}& j\mathrm{mod}{P}^{\left(t\right)}=q\\ 0& \mathrm{else}& \end{array}$
to the vector with entries' index corresponding to those of the support of {circumflex over (ƒ)}. In particular, each such matrix is sparse with exactly #S=O(R) nonzero entries.

(75) Eq. (10) can further be written in matrix form as:

(76) $\begin{array}{cc}\left[FB\right]\hat{f}=\left[\begin{array}{cccc}{F}^{\left(1\right)}& 0& \mathrm{.Math.}& 0\\ 0& F^{\left(2\right)}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& F^{\left(T\right)}\end{array}\right]\left[\begin{array}{c}{B}^{\left(1\right)}\\ B^{\left(2\right)}\\ \mathrm{.Math.}\\ B^{\left(T\right)}\end{array}\right]\hat{f}=\left[\begin{array}{c}{f}^{\left(1\right)}\\ f^{\left(2\right)}\\ \mathrm{.Math.}\\ f^{\left(T\right)}\end{array}\right]=f_0,& \left(11\right)\end{array}$
where F.sup.(t) is a standard DFT matrix of size P.sup.(t). Proposition 6 states that if T=O(1) is sufficiently large, then with nonzero probability

(77) $\frac{1}{T}\left(\mathrm{FB}\right)*\left(\mathrm{FB}\right)=I+,$
where I is the identity matrix and is a perturbation with 2-norm smaller than ½. When this occurs, one can solve a linear system through the Neumann series, i.e.,

(78) $\hat{f}=\underset{n=0}{\overset{\infty }{.Math.}}{\left(I-B*B\right)}^n\left(\mathrm{FB}\right)*f_0$
This forms the core of Algorithm 4. Proposition 3 provides the correctness of this algorithm.

(79) Since each matrix B.sup.(t) contains exactly R nonzero entries, both B and B*B can be applied in order RT=(R log.sub.R (N)) steps. In addition, since F is a block diagonal matrix with T=O(1) blocks consisting of DFT matrices of size (R log.sub.R(N)), it can be applied in order (R log(N)) using the FFT computation. Finally, for an accuracy η the Neumann series can be truncated after (log( )) terms, and the process may be repeated at most log(p) times for a probability p of success. Therefore, the cost of computing the nonzero values of {circumflex over (ƒ)} is bounded by ((log(p)+log(η))R log(N)) and uses at most ((log(p)+log(η))R log.sub.R(N)) samples.

(80) Proposition 3. Given the support of {circumflex over (ƒ)} is known, Algorithm 4 outputs an approximation to the nonzero elements of {circumflex over (ƒ)} with error bounded by η in the .sup.2-norm, with probability greater than or equal to 1−p using ((log(p)+log(η))R log.sub.R(N)) samples and ((log(p)+log(η))R log(N)) computational steps.

(81) Stability to Low-Level Noise

(82) As discussed previously, the theory underlying Algorithms 1 through 4 has been designed for vectors that are exactly sparse. Here, we discuss the effect of low-level noise. In fact, we show that if the sparse vector of Fourier coefficients takes the form {circumflex over (ƒ)}+{circumflex over (ν)}, where √{square root over (N)}∥{circumflex over (ν)}∥.sub.2<η for some “small” η, an embodiment of the sMFFT technique can recover the support and values of {circumflex over (ƒ)} with the same guarantees as described earlier.

(83) Support Recovery

(84) An important quantity for the fast recovery of the support is the set of frequency coefficients of shuffled and filtered samples, as described in Eq. (9), so in the presence of noise,

(85) $\begin{array}{cc}{\varphi }_n^{\left(k\right)}\left(Q\right)=\frac{1}{M_k}\underset{m\in \left(K,M_k\right)}{.Math.}{e}^{2\pi i\frac{nm}{K}}{g}_{\sigma }\left(\frac{m}{M_k}\right)\left(f_{\left({m\left[Q\right]}_{M_k}^{-1}\right)\mathrm{mod}{M}_k;M_k}+v_{\left({m\left[Q\right]}_{M_k}^{-1}\right){\mathrm{modM}}_k;M_k}\right).& \left(12\right)\end{array}$
The second term in this expression is the error term and can be uniformly bounded by the following lemma:

(86) Lemma 1. Assuming the noise term {circumflex over (ν)} is such that

(87) ${.Math.\hat{v}.Math.}_2<\frac{\eta }{\sqrt{N}},$
the error term of the computed value in Eq. (12) is uniformly bounded by the error term of the computed value in

(88) ${.Math.{\psi }_n^{\left(k\right)}\left(Q\right).Math.}_{\infty }={.Math.\frac{1}{M_k}\underset{m\in \left(K,M_k\right)}{.Math.}{e}^{2\pi i\frac{nm}{K}}{g}_{\sigma }\left(\frac{m}{M_k}\right)v_{\left({m\left[Q\right]}_{M_k}^{-1}\right)\mathrm{mod}{M}_k;M_k}.Math.}_{\infty }<\mathrm{��}\left(\eta \right).$
Algorithm 2 tests whether

(89) $.Math.{\varphi }_{\left[\left(i{Q}^{\left(l\right)}mod{M}_k\right)\frac{K}{M_k}\right]\frac{M_k}{K}}^{\left(k\right)}\left(Q^{\left(l\right)}\right).Math.>\mathrm{\delta \mu }$
in order to discriminate between elements of the candidate and aliased supports. The presence of noise can skew this test in two ways: 1) by bringing the computed value below the threshold when i∈.sub.k, or 2) by bringing the value above the threshold multiple times when i.Math..sub.k. Either way, if η is small enough, i.e., such that

(90) 0${.Math.{\psi }_n^{\left(k\right)}\left(Q^{\left(l\right)}\right).Math.}_{\infty }\le \frac{\delta \mu }{2},$
the conclusion of Proposition 5 follows through with similar estimate, by simply replacing δ with

(91) $\frac{\delta }{2}$
in the proof.
Recovering Values from Knowledge of the Support

(92) It is quickly observed that the recovery of the values is a well-conditioned problem. Indeed, since

(93) $\frac{1}{T}\left(\mathrm{FB}\right)*\left(\mathrm{FB}\right)=I-,$
and

(94) ${.Math..Math.}_2\le \frac{1}{2}$
with high probability by Proposition 6, a simple argument based on the singular value decomposition produces the following corollary,

(95) Corollary 1. Under the hypothesis of Proposition 6,

(96) ${\left(\frac{1}{T}\left(\mathrm{FB}\right)*\left(\mathrm{FB}\right)\right)}^{-1}$
exists, and

(97) ${.Math.{\left(\frac{1}{T}\left(\mathrm{FB}\right)*\left(\mathrm{FB}\right)\right)}^{-1}.Math.}_2\le 2$
with probability greater than or equal to ½.
Therefore, the output of Algorithm 4 is such that

(98) ${.Math.{\hat{f}}^{\mathrm{sMFFT}}-\hat{f}.Math.}_2\le {.Math.{\left(\frac{1}{T}\left(\mathrm{FB}\right)*\left(\mathrm{FB}\right)\right)}^{-1}.Math.}_2{.Math.\left(\mathrm{FB}\right)*v.Math.}_2\le 2{.Math.B.Math.}_2{.Math.v.Math.}_2=\mathrm{��}\left(\eta \right).$
This, together with Proposition 3, demonstrates the stability of Algorithm 4 in the noisy case.
The Multi-Dimensional Sparse FFT

(99) Whenever dealing with the multidimensional DFT/FFT, it is assumed that the function of interest is both periodic and bandlimited with fundamental period [0,1).sup.d i.e.,

(100) $f\left(x\right)=\underset{j\in {\left(\left[0,M\right].Math.\right)}^d}{.Math.}{e}^{-2\pi \mathrm{ix}.Math.j}\widehat{f_j}$
for some finite M∈ and j∈.sup.d, up to some rescaling. Computing the Fourier coefficients is then equivalent to computing the d-dimensional integrals,
{circumflex over (ƒ)}.sub.n=∫.sub.[0,1].sub.de.sup.−2πin.Math.xƒ(x)dx,
and this can be achieved through a “dimension-by-dimension” trapezoid rule

(101) $\begin{array}{cc}{\hat{f}}_{\left(j_1,j_2,\mathrm{.Math.},j_d\right)}=\underset{n_1=0}{\overset{M-1}{.Math.}}\frac{e^{2\pi i\frac{j_1{n}_1}{M}}}{M}\left(\mathrm{.Math.}\left(\underset{n_{d-1}=0}{\overset{M-1}{.Math.}}\frac{e^{2\pi i\frac{j_{d-1}{n}_{d-1}}{M}}}{M}\left(\underset{n_d=0}{\overset{M-1}{.Math.}}\frac{e^{2\pi i\frac{j_d{n}_d}{M}}}{M}{f}_{\left(n_1,n_2,\mathrm{.Math.},n_d\right)}\right)\right)\right).& \left(13\right)\end{array}$

(102) However, Proposition 4 below shows that it is also possible to re-write the d-dimensional DFT as that of a 1D function with Fourier coefficients equal to those of the original function, but with different ordering.

(103) Proposition 4. (Rank-1 d-dimensional DFT) Assume the function ƒ:[0,1).sup.d.fwdarw. has form (13). Then,

(104) $\begin{array}{cc}{\int }_{{\left[0,1\right]}^d}{e}^{-2\pi \mathrm{ij}.Math.x}f\left(x\right)dx=\frac{1}{N}\underset{n=0}{\overset{N-1}{.Math.}}{e}^{-2\pi \mathrm{ij}.Math.x_n}f\left(x_n\right)& \left(14\right)\end{array}$
for all j∈[0, M).sup.d ∩.sup.d, where

(105) 0$x_n=\frac{\mathrm{ng}\mathrm{mod}N}{N},g=\left(1,M,M^2,\mathrm{.Math.},M^{d-1}\right)$
and N=M.sup.d.
Now, the right-hand side of Eq. (14) can be written in two different ways (due to periodicity); namely,

(106) $\frac{1}{N}\underset{n=0}{\overset{N-1}{.Math.}}{e}^{-2\pi ij\frac{\mathrm{ng}\mathrm{modN}}{N}}f\left(\frac{\mathrm{ng}\mathrm{mod}N}{N}\right)=\frac{1}{N}\underset{n}{.Math.}{e}^{-2\pi i\frac{\left(j.Math.g\right)n}{N}}f\left(\frac{ng}{N}\right)$

(107) Geometrically, the left-hand side represents a quadrature rule with points

(108) $x_n=\frac{\mathrm{ng}\mathrm{mod}N}{N}$
distributed (more-or-less uniformly) in [0, 1).sup.d, as shown at points 802 in FIG. 8A. The right-hand side represents an equivalent quadrature where the points

(109) $x_n=\frac{ng}{N}$
now lie on a line embedded in R.sup.d, as shown at points 806 in FIG. 8B. The location at which the lattice (thin black lines) intersects represents the standard multidimensional DFT samples. Thus, Proposition 4 allows one to write any d-dimensional DFT as a one-dimensional DFT by picking the appropriate sample points (Proposition 4) and proceeding to a re-ordering of the Fourier coefficients through the isomorphism
ń:{n∈([0,M)∩).sup.D}.fwdarw.n.Math.g=n.sub.0+n.sub.1M+ . . . +n.sub.D−1M.sup.D−1∈[0,M.sup.D)∩.sup.D.

(110) We use this sort of sampling to treat of the multidimensional problem, which we first convert to a 1D problem. Then, we address the 1D problem through Algorithms 1 through 4, and then map the result back to its original d-dimensional space. Specifically, the sampling points corresponding to a 1-D signal that is equivalent to the multi-dimensional signal are obtained using Algorithms 1 through 4. Using inverse mapping between the multi-dimensional signal and the 1-D signal, the corresponding points of the multi-dimensional signal are identified and sampled.

(111) FIGS. 8A and 8B generally describe geometric interpretation of rank-1 d-dimensional DFT in 2D. The thick black boxes represent fundamental periodic domain. The dots 802, 804 represent rank-1 discretization points. The 2D grid represents standard discretization points. In FIG. 8A, the rank-1 d-dimensional quadrature interpreted as a 2D discretization over the fundamental periodic region. In FIG. 8B, the rank-1 d-dimensional quadrature interpreted as a uniform discretization over a line in R.sup.2.

(112) In general, there is no need for the signal's spectrum to reside on a regular grid. Indeed if the signal is sparse in the continuous domain, one can approximate such function on the hypercube by a function which spectrum does lie on a regular grid, and which spectral support correspond to the closest neighbor of the original spectral support within the regular grid. For some error ϵ, a grid with spacing

(113) $\mathrm{��}\left(\frac{1}{\mathrm{.Math.}}\right)$
containing

(114) $N=\mathrm{��}\left(\frac{1}{{\mathrm{.Math.}}^d}\right)$
unknowns should suffice.
An Example System for Generating Frequency Components of a Signal

(115) With reference to FIG. 9A, an sMFFT generation system 900 can generate the non-zero frequency components of a received signal ƒ(x). The signal ƒ(x) can be received from one or more signal sources such as a radar, a lidar, an antenna, a sonar, a camera, an infra-red sensor, an electro-magnetic radiation sensor, or an acoustic sensor. In some implementations, a signal received from a source may be preprocessed, e.g., filtered to limit the bandwidth of the signal. The bandwidth of the signal ƒ(x) and/or a maximum frequency corresponding to the signal ƒ(x) may be known or estimated. Using the bandwidth and/or the maximum frequency, the number of samples N that must be collected during a window of length τ seconds so as to facilitate the conventional DFT/FFT transformation of the signal ƒ(x) is determined. In general, if the maximum frequency is ƒ.sup.max according to the Nyquist sampling rate N≥2*ƒ.sup.max*τ. N can also be the total number of frequency coefficients of the signal ƒ(x). In addition, an upper bound R on the number of non-zero frequency components of the signal ƒ(x) is also provided to the system 900 or is estimated.

(116) Using the parameters N and R, one or more of other parameters of the system may be computed or selected. These parameters include an FFT size K, a number of iterations P, an index limit M.sub.k for the candidate support set .sub.k used in the k-th iteration, and a support set growth factor ρ.sub.k. In various implementations, these parameters may be related according to the expressions

(117) $N=K\underset{k=1}{\overset{P}{.Math.}}{\rho }_k;$
and M.sub.k=Kρ.sub.k=O(R log R). The system 900 may also determine various shuffling parameters Q.sub.k.sup.l, where l ranges from 1 through a parameter L, using the values M.sub.k. In addition, the system 900 may determine, using the values M.sub.k, various value recovery parameters P.sup.(l,t), where l ranges from 1 through the parameter L and t ranges from 1 through a parameter T, discussed above. In some implementations, L=1. The shuffling and value recovery parameters can be used for sampling the signal ƒ(x).

(118) In particular, the system 900 includes a sampling/sub-sampling module 902 and a sampling and shuffling module 904. In some implementations, the shuffling module can be different module. The sampling/sub-sampling module 902 samples the signal ƒ(x) using the parameters M.sub.1 and N, and generates a set of samples ƒ.sup.1(x.sup.1). A Fourier transform of these samples (e.g., an M.sub.1-point FFT) is performed in the module 906 to obtain M.sub.1 frequency coefficients, which are designated as the candidate support set .sub.1. The non-zero coefficients from the set .sub.1 are selected in the module 908 and designated as the support set .sub.1.

(119) In a de-aliasing module 910, the value of M.sub.k is used to de-alias a previously generated support set .sub.k−1, as discussed above. For example, the support set .sub.1 can be de-aliased to obtain the candidate support set .sub.2. De-aliasing can be performed in the frequency domain, i.e., using indices of the aliased frequency components, and without requiring values of those components or samples of the signal ƒ(x). Thereafter, the significant non-zero coefficients of the candidate support set .sub.k are selected in the module 912, to obtain the support set S.sub.k. This step requires one or more sets of signal samples, which may have been obtained a priori or may be obtained as needed by the sampling/shuffling module 904. During this process, the index limit M.sub.k of the candidate support set .sub.k grows according to the growth factor ρ.sub.k.

(120) The operations in the modules 910, 912 are repeated using the new support set S.sub.k if it is determined that M.sub.k<N. These iterations are generally described above with reference to Algorithms 2 and 3. After P iterations, i.e., when M.sub.k is no longer less than N or is equal to N, the support set S.sub.p generated in the last iteration is considered to be the final support set . The values of the frequency components corresponding to the indices in the final support set are determined at the recovery module 914 according to Algorithm 4 discussed above. The determination of the values also requires one or more sets of samples of the signal ƒ(x). These sample set(s) can be obtained a priori or as needed by the sampling/shuffling module 904.

(121) FIG. 9B shows a sampler 952 that uses the parameters M.sub.1 and N to obtain samples of ƒ(x) that can be used to generate the candidate support set .sub.1, as discussed above. The sampling/sub-sampling module 902 (FIG. 9A) can be or may include the sampler 952. The sampling/shuffling module 904 (FIG. 9A) may include one or more of the other modules depicted in FIG. 9B.

(122) Specifically, the modules 954 use the shuffling parameters Q.sub.k.sup.l, where l ranges from 1 through the parameter L (which can have the value 1), and k ranges from 1 or 2 through the number of iterations P. The values of Q.sub.k.sup.l depend on the values of the index limits M.sub.k, as described above. A distinct set of samples of the signal ƒ(x) may be obtained for each value of Q.sub.k.sup.l. In some implementations, instead of using a different sampler 954 for each Q.sub.k.sup.l, fewer samplers, and as few as a single sampler 954 may be used. One or more samplers 954 may be tunable to perform sampling according to one or more values Q.sub.k.sup.l. In some embodiments, the samplers 954 also include the corresponding shufflers. In other implementations, shuffling, as described above, may be performed in different module(s).

(123) The modules 956 use the recovery parameters P.sup.(l,t), where l ranges from 1 through the parameter L (which can have the value 1), and t ranges from 1 through T, as described above. The values of P.sup.(l,t) depend on the parameters N and/or R, as described above. Here again, a distinct set of samples of the signal ƒ(x) may be obtained for each value of P.sup.(l,t). In some implementations, instead of using a different sampler 956 for each P.sup.(l,t), fewer samplers, and as few as a single sampler 956 may be used. One or more samplers 956 may be tunable to perform sampling according to one or more values P.sup.(l,t). In some embodiments, the samplers 956 also include the corresponding shufflers. In other implementations, shuffling, as described above, may be performed in different module(s).

(124) In some implementations, all the sampling is performed a priori, i.e., before the determination of any or at least some of the support sets .sub.k. The different sample sets are stored in memory and are accessed as needed. In other implementations, the sampling and the computation of the support sets and/or the recovery of the frequency coefficients may be interleaved. In some implementations, a set of samples is obtained by performing sampling only when it is needed for the determination of the support set or for recovery of the frequency coefficient values.

(125) Generalization to General Complex Sparse Vectors

(126) We now describe certain additional steps for transforming the sMFFT for real positive vectors into a reliable process for general complex vectors. To achieve this task, two major hurdles, both associated with the support-recovery portion of the scheme, are overcome; the first one is associated with the initial aliasing of the signal described, as discussed above. As shown in Eq.(5), at each step aliasing implies Fourier coefficients of the form,
{circumflex over (ƒ)}.sub.l.sup.(k)=Σ.sub.j:j mod M.sub.k.sub.=l{circumflex over (ƒ)}.sub.j,l=0,1, . . . ,M.sub.k.
When the original nonzero coefficients are all strictly positive, this expression is positive if and only if the lattice

(127) $l+i{M}_k,i=0,1,\mathrm{.Math.}\frac{N}{M_k}-1$
contains one of the original non-zero coefficients. When the non-zero coefficients are complex, however, this is no longer true.

(128) The second potential issue pertains to the resulting filtering step found in Algorithm 3. As described by Eq. (8), the result takes the form,

(129) $\left[{\Psi }_{\sigma }\left({\Pi }_Q\left(f_{n;M_k}\right)\right)\right]\left(\xi \right)=F^{*}\left[\underset{j\in k_{}}{.Math.}{\widehat{f_j}}^{\left(k\right)}{e}^{-\frac{{.Math.x-\left({j\left[Q\right]}_{M_k}^{-1}mod{M}_k\right).Math.}^2}{{\sigma }^2}}\right]\left(\xi \right).$
which corresponds to the Fourier transform of the aliased signal convoluted with a Gaussian. Once again, the statistical test used in Algorithm 3 relies on this quantity being positive if and only if a point lies in the vicinity of an element of the (shuffled and aliased) support .sub.k. Such statement does not hold true if we allow the coefficients to be general complex numbers (as some elements might cancels out).

(130) It follows from these observations that as a consequence of the lack of positivity, it is possible that elements belonging to .sub.k∩.sub.k might be wrongfully eliminated in Algorithm 3, i.e., the false negative identification rate is nontrivial. To alleviate these issues, we describe a modified scheme, where we accommodate for the possibility that the output of Algorithm 3 may be missing one or more elements of .sub.k, by launching multiple independent iterations of the FIND_SUPPORT(⋅) process in Algorithm 1, and by taking the union of the outputs. In this sense, although it is possible to miss an element with a single run, we expect that the probability of a miss over multiple independent run is very small. In one implementation, up to log p independent iterations, using different values of N, are launched.

(131) This modification does not have any effect on the fundamental computational complexity; indeed, close examination shows that these additional steps only increase the algorithmic constant by some small quantity independent of N and/or R. We have implemented this modification and it provided good numerical results in line with our expectation based on the previous discussion, and very similar to those obtained in the real-positive case.

(132) Numerical Results

(133) We have implemented our sMFFT algorithms in MATLAB® and present a few numerical results which exhibit the efficient scaling. All simulations were carried out on a small cluster possessing 4 Intel Xeon E7-4860 v2 processors and 256 GB of RAM, with the MATLAB flag—singleCompThread to ensure fairness through the use of a single computational thread. The numerical experiments presented here fall in two categories: 1) dependence of running time as a function of the total number of unknowns N for a fixed number of nonzero frequencies R, and 2) dependence of running time as a function of the number of nonzero frequencies R for a fixed total number of unknowns N. All experiments were carried out in three dimensions (3D) with additive Gaussian noise with variance r. The nonzero values of {circumflex over (ƒ)} were picked randomly and uniformly at random in [0.5, 1.5], and the remaining parameters were set according to Table 2, shown in FIG. 10. All comparisons are perdormed with the MATLAB fftn(⋅) function, which uses a dimension-wise decomposition of the DFT and a 1D FFT routine along each dimension.

(134) For case 1), we picked R=50 nonzero frequencies distributed uniformly at random on a 3D lattice having N.sup.1/3 elements in each dimension for different values of N∈[10.sup.3, 10.sup.10]. The results are shown in FIG. 11A, which depicts the running time vs. the number of unknowns (N) for the MATLAB fftn(−) and the sMFFT in three dimensions (3D), with R=50 nonzeros and noise n=10.sup.−3. It can be observed that the cost of computing the DFT through the sMFFT remains more or less constant with N, whereas that the cost of the MATLAB fftn(⋅) function increases linearly. This is the expected behavior and demonstrates the advantages of the sMFFT technique over the FFT. Also note that the largest relative .sup.2-error observed was 9. 3 10.sup.−3 which is on the order of the noise level, as predicted by the theory.

(135) For case 2), we fixed N=O(10.sup.8) and proceeded to compare the sMFFT algorithm with the MATLAB fftn(−) function as before (w/ parameters found in Table 2). The results are shown in FIG. 11B, which depicts the running time vs number of nonzero frequencies (R) for the MATLAB fftn(−) and the sMFFT in three dimensions (3D) and for N=10.sup.8 and noise n=10.sup.−3. In this case, the theory states that the sMFFT algorithm should scale quasi-linearly with the number of nonzero frequencies R. A close look shows that it is indeed the case. For this case, the largest relative .sup.2-error observed was 1.1.Math.10.sup.−2, again on the order of the noise level and in agreement the theory. Finally, the cost of the MATLAB fftn(−) function remains constant as the FFT scales like (N log(N)) and is oblivious to R.

(136) We have introduced a sparse multidimensional FFT (sMFFT) for computing the DFT of a N×1 sparse, real-positive or complex vector (having R non-zeros) that is stable to low-level noise, that exhibits a sampling complexity of ((R log(R)log(N))) and a computational complexity ((R log.sup.2(R)log(N))). This technique also has a scaling of (log(p)) with respect to probability of failure, (log(η)) with respect to accuracy and O(1) with respect to dimension. We have provided a rigorous theoretical analysis of our approach demonstrating this efficient and accurate computation of the non-zero frequency coefficients of a signal. We have implemented our algorithm and provided numerical examples in 3D successfully demonstrating the benefits of the scaling that can be achieved using implementations of our technique.

(137) Proofs

(138) We present all proofs and accompanying results related to the statements presented in the discussion above.

(139) Lemma 2. Let 0<Q≤N∈, Q⊥N. Then the map,
n∈{0,1, . . . ,N−1}.fwdarw.nQ mod N⊂{0,1, . . . ,N−1}

(140) is an isomorphism.

(141) Proof. Since the range is discrete and a subset of the domain, it suffices to show that the map is injective. Surjectivity will then follow from the pigeon hole principle. To show injectivity, consider i, j∈{0, 1, . . . , N−1}, and assume,
iQ mod N=jQ mod N

(142) This implies (by definition) that there exists some integer p such that,
(i−j)Q=pN

(143) so that N divides (i−j)Q. However, N⊥Q so N must be a factor (i−j). Now, i, j are restricted to {0, 1, . . . , N−1} so,
|i−j|<N,

(144) and the only integer divisible by N that satisfies this equation is 0. Thus,
i−j=0⇔i=j

(145) which demonstrates injectivity.

(146) Lemma 3. Let 0<Q<M be an integer coprime to M and,

(147) $f_n=\underset{l=0}{\overset{M-1}{.Math.}}{e}^{-2\pi i\frac{nl}{M}}\widehat{f_l}$

(148) Then,

(149) 0$\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{mn}{M}}{f}_{\left(nQ\right)\mathrm{mod}M}={\hat{f}}_{\left({m\left[Q\right]}_M^{-1}\right)modM}$

(150) where 0<[Q].sub.M.sup.−1<M is the unique integer such that [Q].sub.M.sup.−1 Q mod M=1 mod M.

(151) Proof. Consider

(152) $\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{mn}{M}}{f}_{\left(nQ\right)\mathrm{mod}M}\underset{l=0}{\overset{M-1}{.Math.}}\left(\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{n}{M}\left(m-\mathrm{Ql}\mathrm{mod}M\right)}\right){\hat{f}}_l=\underset{l=0}{\overset{M-1}{.Math.}}\left(\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{\mathrm{nQ}}{M}\left({m\left[Q\right]}_M^{-1}\mathrm{mod}M-l\right)}\right){\hat{f}}_l$

(153) However,

(154) $\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{\mathrm{nQ}}{M}\left({m\left[Q\right]}_M^{-1}\mathrm{mod}M-l\right)}=\frac{1}{M}\underset{j=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{j}{M}\left({m\left[Q\right]}_M^{-1}\mathrm{mod}M-l\right)}={\delta }_{{m\left[Q\right]}_M^{-1}\mathrm{mod}M,l},$

(155) where the second equality follows from the fact that m.fwdarw.m[Q].sub.M.sup.−1 mod M is an isomorphism (Lemma 2).

(156) This implies that

(157) $\frac{1}{M}\underset{n=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{mn}{M}}{f}_{\left(nQ\right)modM}={\hat{f}}_{\left({m\left[Q\right]}_M^{-1}\right)modM}$

(158) as claimed.

(159) Lemma 4. Let M∈/{0} and let Q be a uniform random variable over Q(M) (Definition 2). Then,

(160) $\left(.Math.\mathrm{jQ}\mathrm{mod}M.Math.\le C\right)\le \mathrm{��}\left(\frac{C}{M}\right)$

(161) for all 0<j<M (up to a log(log(M)) factor).

(162) Proof. Fix 0<j, k<M and let γ=gcd(j, M). Consider,

(163) $\left(\mathrm{jQ}\mathrm{mod}M=k\right)=\underset{q\in Q\left(M\right)}{.Math.}\left(\mathrm{jq}\mathrm{mod}M=k.Math.Q=q\right)\left(Q=q\right)=\underset{q\in Q\left(M\right)}{.Math.}{\mathrm{jq}\mathrm{modM}=k}_{}\left(q\right)\left(Q=q\right)$

(164) and note that,

(165) $\left(Q=q\right)=\frac{1}{\#Q\left(M\right)}=\frac{1}{\varphi \left(M\right)}\le \frac{e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}}{M}$

(166) following bounds on the Euler totient function ϕ(⋅) ([33]), where ζ is the Euler-Mascheroni constant, and since Q is uniformly distributed in Q(M). Therefore,

(167) $\left(\mathrm{jQ}\mathrm{mod}M=k\right)\le \frac{e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}}{M}\underset{q}{.Math.}{\mathrm{jq}\mathrm{modM}=k}_{}\left(q\right)$

(168) We now show that the quantity Σ.sub.q.sub.jq mod M=k(q) is bounded above and below by,

(169) $\gamma -1\le \underset{q}{.Math.}$

(170) To see this, first note that this quantity corresponds to the number of integers q which hash to the integer k through the map q.fwdarw.(jq) mod M. Now, assume there exists some q such that
jq mod M≡k,  (15)
which implies that
jq+iM=k  (16)

(171) for some integer i∈. This is a Diophantine equation which has infinitely many solutions if and only if gcd(j, M)=y divides k. Otherwise, it has no solution. Assuming it does and (q.sub.0, i.sub.0) is a particular solution, all remaining solutions must take the form

(172) $q=q_0+u\frac{M}{\gamma },i_0-u\frac{j}{y}$

(173) where u∈. However, since 0≤q<M the number of possible solutions must be such that,

(174) 0$\gamma -1\le \#\left\{q\in \left[0,M\right):q=q_0+u\frac{M}{\gamma }\right\}\le \gamma .$
which proves the claim. Thus,

(175) $\left(\mathrm{jQ}\mathrm{mod}M=k\right)\le \left(e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}\right)\frac{\gamma }{M}$

(176) We can now treat (|jQ mod M|≤C). Before we proceed however, recall that Eq.(15) has a solution if and only if γ|k. We then write,

(177) $\left(.Math.\mathrm{jQ}\mathrm{mod}M.Math.\le C\right)=\underset{0\le k\le C}{.Math.}{\gamma |k}_{}\left(k\right)\left(\mathrm{nQ}\mathrm{mod}M=k\right)$

(178) from which it follows that,

(179) $\left(.Math.\mathrm{jQ}\mathrm{mod}M.Math.\le C\right)\le \left(e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}\right)\frac{\gamma }{M}\underset{0\le k\le C}{.Math.}{\gamma |k}_{}\left(k\right)\le \left(e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}\right)\frac{\gamma }{M}\frac{C}{\gamma }\le \left(e^{\zeta }\log \left(\log \left(M\right)\right)+\frac{3}{\log \left(\log \left(M\right)\right)}\right)\frac{C}{M}$

(180) since the number of integers in 0≤k≤C that are divisible by γ is bounded above by

(181) $\frac{C}{\gamma }.$
Finally, since this holds regardless of our choice of j, this proves the desired result.

(182) Lemma 5. Consider a function ƒ(x) of the form of Eq.(3) and satisfying a certain constraint. Let

(183) $0<\sigma =\mathrm{��}\left(\frac{M}{R\sqrt{\log \left(R\right)}}\right),0<\delta <1\mathrm{and},\mu =\underset{j\in }{\min }.Math.\widehat{f_j}.Math.$$\Delta =\frac{\underset{j\in }{\max }.Math.\widehat{f_j}.Math.}{\underset{j\in }{\min }.Math.\widehat{f_j}.Math.}=\frac{{.Math.\hat{f}.Math.}_{\infty }}{\mu }$

(184) Finally, let (⋅) and Ψ.sub.σ(⋅) be the operators found in Eq.(9). Then, there exists a constant 1<C<∞ such that if

(185) $\begin{array}{cc}K\ge C\frac{M\sqrt{\log \left(\frac{2\Delta }{\delta }\right)}}{\pi \sigma }& \left(17\right)\end{array}$

(186) the inequality,

(187) ${\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left(f_{k;M}\right)\right)\right]}_{n\frac{M}{K}}\ge \delta \mu$

(188) implies that

(189) $\begin{array}{cc}\underset{j\in }{\inf }.Math.n\frac{M}{K}-j.Math.\le \sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}& \left(18\right)\end{array}$

(190) and,

(191) $\underset{j\in }{\inf }.Math.n\frac{M}{K}-j.Math.\le \sigma \sqrt{\log \left(\frac{1}{\delta }\right)},$

(192) implies that

(193) 0${\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left(f_{k;M}\right)\right)\right]}_{n\frac{M}{K}}\ge \delta \mu$

(194) for all n∈{0, 1, . . . , K−1}.

(195) Proof. Consider the quantity

(196) $\begin{array}{cc}\frac{1}{M}\underset{m\in \left(K;M\right)}{.Math.}{e}^{2\mathrm{\pi i}\frac{n\frac{M}{K}m}{M}}{\widehat{ℊ}}_{\sigma }\left(\frac{m}{M}\right)f_{m;M}=\underset{j\in }{.Math.}\left(\frac{1}{M}\underset{m\in \left(K;M\right)}{.Math.}{e}^{2\mathrm{\pi i}\frac{m}{K}\left(n-j\right)}{\widehat{ℊ}}_{\sigma }\left(\frac{m}{M}\right)\right){\hat{f}}_j& \left(19\right)\end{array}$

(197) an recall that,

(198) $\frac{1}{M}{\hat{g}}_{\sigma }\left(\frac{m}{M}\right)=\frac{\sqrt{\pi }\sigma }{M}\underset{h\in }{.Math.}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{m+hM}{M}\right)}^2}$

(199) where

(200) $\left(K;M\right):=\left\{m\in \left\{0,1,\mathrm{.Math.},M-1\right\}:m\le \frac{K}{2}\mathrm{or}.Math.m-M.Math.<\frac{K}{2}\right\}$
(Definition 1).
From this expression, it is apparent that there exists some constant 1<C<∞ such that by choosing

(201) $K\ge C\frac{M\sqrt{\log \left(\frac{2\Delta }{\delta }\right)}}{\pi \sigma },$
one has,

(202) $.Math.\underset{j\in }{.Math.}\left(\frac{1}{M}\underset{m\in c^{}\left(K;M\right)}{.Math.}{e}^{2\pi i\frac{m}{K}\left(n-j\right)}{\hat{g}}_{\sigma }\left(\frac{m}{M}\right)\right)\widehat{f_j}.Math.\le \underset{m\in c^{}\left(K;m\right)}{.Math.}\left(\frac{\sqrt{\pi }\sigma }{M}\underset{h\in }{.Math.}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{m+hM}{M}\right)}^2}\right){.Math.\hat{f}.Math.}_{\infty }\le \frac{\delta \mu }{2}$

(203) Indeed, by the integral test,

(204) $\underset{m\in c^{}\left(K;m\right)}{.Math.}\left(\frac{\sqrt{\pi }\sigma }{M}\underset{h\in }{.Math.}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{m+hM}{M}\right)}^2}\right)\le A\frac{\sqrt{\pi }\sigma }{M}\underset{m\ge \frac{K}{2}}{.Math.}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{m}{M}\right)}^2}\le A\frac{\sqrt{\pi }\sigma }{M}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{K}{2M}\right)}^2}+\frac{\sqrt{\pi }\sigma }{M}{\int }_{\frac{K}{2}}^{\infty }{e}^{-{x^2\left(\frac{\pi \sigma }{M}\right)}^2}\mathrm{dx}\le A\frac{\sqrt{\pi }\sigma }{M}{e}^{-{\pi }^2{{\sigma }^2\left(\frac{K}{2M}\right)}^2}+\frac{1}{\sqrt{\pi }}\mathrm{erfc}\left(\frac{\pi K\sigma }{2M}\right)\le B{e}^{-{\pi }^2{{\sigma }^2\left(\frac{K}{2M}\right)}^2}$

(205) for some positive constants A, B, and where the last inequality follows from estimates on the complementary error function and the fact that

(206) $\frac{\sqrt{\pi }\sigma }{M}=\mathrm{��}\left(\frac{1}{R\sqrt{\log \left(R\right)}}\right)$
by assumption. Therefore,

(207) ${\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left(f_{m;M}\right)\right)\right]}_{n\frac{M}{K}}=\frac{1}{M}\underset{m\in \left(K;M\right)}{.Math.}{e}^{2\pi i\frac{nm}{K}}{\hat{g}}_{\sigma }\left(\frac{m}{M}\right)f_{m;M}=\frac{1}{M}\underset{m=0}{\overset{M-1}{.Math.}}{e}^{2\pi i\frac{n\frac{M}{K}m}{M}}{\hat{g}}_{\sigma }\left(\frac{m}{M}\right)f_{m;M}+{\mathrm{.Math.}}_n=\underset{j\in }{.Math.}{e}^{-\frac{\left(n\frac{M}{K}-j\right)}{{\sigma }^2}}\widehat{f_j}+{\mathrm{.Math.}}_n$

(208) where

(209) ${\max }_n.Math.{\mathrm{.Math.}}_n.Math.\le \frac{\delta \mu }{2}.$
Now assume:

(210) 0$.Math.{\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left(f_{k;M}\right)\right)\right]}_{n\frac{M}{K}}.Math.\ge \mathrm{\delta \mu }.$
Then, the triangle inequality and the previous equation imply that,

(211) $\begin{array}{cc}\underset{j\in }{.Math.}{e}^{-\frac{\left(n\frac{M}{K}-j\right)}{{\sigma }^2}}\widehat{f_j}\ge .Math.{\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left(f_{m;M}\right)\right)\right]}_{n\frac{M}{K}}.Math.-\frac{\delta \mu }{2}\ge \left(\delta -\frac{\delta }{2}\right)\mu =\frac{\delta }{2}\mu .& \left(20\right)\end{array}$

(212) We claim that this cannot occur unless,

(213) $\begin{array}{cc}\underset{j\in }{\inf }.Math.n\frac{M}{K}-j.Math.\le \sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}.& \left(21\right)\end{array}$

(214) We proceed by contradiction. Assume the opposite holds. Then,

(215) $\underset{j\in }{.Math.}{e}^{-\frac{\left(n\frac{M}{K}-j\right)}{{\sigma }^2}}\widehat{f_j}\le {.Math.\hat{f}.Math.}_{\infty }\underset{j\in }{.Math.}{e}^{-\frac{{\left(n\frac{M}{K}-j\right)}^2}{{\sigma }^2}}<{.Math.\hat{f}.Math.}_{\infty }\frac{\delta }{2\Delta }=\frac{\delta \mu }{2}$

(216) by assumption. This is a contradiction. Thus, Eq.(21) must indeed hold. This proves the first part of the proposition. For the second part, assume

(217) $\begin{array}{cc}\underset{j\in }{\inf }.Math.n\frac{M}{K}-j.Math.\le \sigma \sqrt{\log \left(\frac{1}{\delta }\right)}& \left(22\right)\end{array}$

(218) holds. Letting j* be such that

(219) $.Math.n\frac{M}{K}-j^{*}.Math.={\inf }_{j\in }.Math.n\frac{M}{K}-j.Math.,$
we note that,

(220) $\underset{j\in }{.Math.}{e}^{-\frac{{\left(n\frac{M}{K}-j\right)}^2}{{\sigma }^2}}\widehat{f_j}\ge e^{-\frac{{\left(n\frac{M}{K}-j^{*}\right)}^2}{{\sigma }^2}}{\hat{f}}_{j^{*}}\ge \mathrm{\delta \mu }$ since ƒ and the Gaussian are all positive by assumption. This shows the second part. We are now ready to prove the validity of the Algorithm 1.

(221) Proposition 5. (Correctness of Algorithm 3) Consider a function ƒ(x) of the form of Eq.(3) and satisfying the nonnegativity hypothesis, and let Π.sub.Q(⋅), Ψ.sub.σ(⋅) and F.sub.K(⋅) be the operators found in Eq.(9) where δ, μ, Δ and, K are as in Lemma 5 and, K satisfies the additional constraint

(222) $K>\frac{R}{\alpha }\sqrt{\frac{\log \left(\frac{2R\Delta }{\delta }\right)}{\log \left(\frac{1}{\delta }\right)}},$
and

(223) $\begin{array}{cc}\sigma =\frac{\alpha \frac{M}{2R}}{\sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}}& \left(23\right)\end{array}$

(224) for some 0<α<1. Assume further that the integers {Q.sup.(l)}.sub.l=1.sup.L are chosen independently and uniformly at random within Q(M), for some 1≤L∈. Consider

(225) $\begin{array}{cc}{\varphi }_{\left[i\frac{K}{M}\right]\frac{M}{K}}\left(Q^{\left(l\right)}\right):={\left[F_{\left(K;M\right)}\left({\Psi }_{\sigma }\left({\Pi }_{Q^{\left(l\right)}}\left(f_{k;M}\right)\right)\right)\right]}_{\left[i\frac{K}{M}\right]\frac{M}{K}}& \left(24\right)\end{array}$
Then,

(226) 00$\begin{array}{cc}\left({.Math.}_{l=1}^L\left\{.Math.{\varphi }_{\left[i\frac{K}{M}\right]\frac{M}{K}}\left(Q^{\left(l\right)}\right).Math.\ge \mathrm{\delta \mu }\right\}\right)\le {\alpha }^L& \left(25\right)\end{array}$
for every i such that (i[Q].sub.M.sup.−1) mod M.Math..sup.c, and

(227) 01$.Math.{\varphi }_{\left[i\frac{K}{M}\right]\frac{M}{K}}\left(Q^{\left(l\right)}\right).Math.\ge \mathrm{\delta \mu }$
almost surely for all Q.sup.(l) land every i such that (i[Q].sub.M.sup.−1) mod M∈.

(228) Proof. From independence, the probability in Eq. (25) is equal to,

(229) 02$\underset{l=1}{\overset{L}{.Math.}}\left(.Math.{\varphi }_{\left[\left(\mathrm{iQ}\mathrm{mod}M\right)\frac{K}{M}\right]\frac{M}{K}}\left(Q\right).Math.\ge \mathrm{\delta \mu }\right).$

(230) So it is sufficient to consider a fixed l. Assume first that i[Q].sub.M.sup.−1 mod M.Math.. As a consequence of Lemma 5 and Lemma 3 we have the inclusion,

(231) 03$\left\{.Math.{\varphi }_{\left[i\frac{K}{M}\right]\frac{M}{K}}\left(Q\right).Math.\ge \mathrm{\delta \mu }\right\}.Math.\left\{\underset{j\in }{\inf }.Math.\left[i\frac{K}{M}\right]\frac{M}{K}-\left(\mathrm{jQ}\right)\mathrm{mod}M.Math.\le \sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}\right\}.Math.{.Math.}_{j\in }\left\{.Math.\left(\left(\left({i\left[Q\right]}_M^{-1}\mathrm{mod}M\right)-j\right)Q\right)\mathrm{mod}M.Math.\le \sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}+\frac{M}{2K}\right\},$

(232) which implies that the probability for each fixed l is bounded by,

(233) 04$\left(.Math.{\varphi }_{\left[i\frac{K}{M}\right]\frac{M}{K}}\left(Q^{\left(l\right)}\right).Math.\ge \mathrm{\delta \mu }\right)\le \underset{j\in }{.Math.}\left(.Math.\left(\left({i\left[Q\right]}_M^{-1}\mathrm{mod}M-j\right)Q^{\left(l\right)}\right)\mathrm{mod}M.Math.\le \sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}+\frac{M}{2K}\right)\le \mathrm{��}\left(R\left(\frac{\sigma \sqrt{\log \left(\frac{2R\Delta }{\delta }\right)}+\frac{M}{2K}}{M}\right)\right)=\mathrm{��}\left(\alpha \right),$