SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

This article is part of the series Advanced Statistical Tools for Enhanced Quality Digital Imaging with Realistic Capture Models.

Open Access Research

Universal and efficient compressed sensing by spread spectrum and application to realistic Fourier imaging techniques

Gilles Puy12*, Pierre Vandergheynst1, Rémi Gribonval3 and Yves Wiaux145

Author Affiliations

1 Institute of Electrical Engineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

2 Institute of the Physics of Biological Systems, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

3 Centre de Recherche INRIA Rennes-Bretagne Atlantique, F-35042 Rennes cedex, France

4 Institute of Bioengineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

5 Department of Radiology and Medical Informatics, University of Geneva (UniGE), CH-1211 Geneva, Switzerland

For all author emails, please log on.

EURASIP Journal on Advances in Signal Processing 2012, 2012:6  doi:10.1186/1687-6180-2012-6


The electronic version of this article is the complete one and can be found online at: http://asp.eurasipjournals.com/content/2012/1/6


Received:6 July 2011
Accepted:12 January 2012
Published:12 January 2012

© 2012 Puy et al; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We advocate a compressed sensing strategy that consists of multiplying the signal of interest by a wide bandwidth modulation before projection onto randomly selected vectors of an orthonormal basis. First, in a digital setting with random modulation, considering a whole class of sensing bases including the Fourier basis, we prove that the technique is universal in the sense that the required number of measurements for accurate recovery is optimal and independent of the sparsity basis. This universality stems from a drastic decrease of coherence between the sparsity and the sensing bases, which for a Fourier sensing basis relates to a spread of the original signal spectrum by the modulation (hence the name "spread spectrum"). The approach is also efficient as sensing matrices with fast matrix multiplication algorithms can be used, in particular in the case of Fourier measurements. Second, these results are confirmed by a numerical analysis of the phase transition of the ℓ1-minimization problem. Finally, we show that the spread spectrum technique remains effective in an analog setting with chirp modulation for application to realistic Fourier imaging. We illustrate these findings in the context of radio interferometry and magnetic resonance imaging.

1 Introduction

In this section we concisely recall some basics of compressed sensing, emphasizing on the role of mutual coherence between the sparsity and sensing bases. We discuss the interest of improving the standard acquisition strategy in the context of Fourier imaging techniques such as radio interferometry and magnetic resonance imaging (MRI). Finally we highlight the main contributions of our study advocating a universal and efficient compressed sensing strategy coined spread spectrum, and describe the organization of this article.

1.1 Compressed sensing basics

Compressed sensing is a recent theory aiming at merging data acquisition and compression [1-7]. It predicts that sparse or compressible signals can be recovered from a small number of linear and non-adaptative measurements. In this context, Gaussian and Bernouilli random matrices, respectively with independent standard normal and ± 1 entries, have encountered a particular interest as they provide optimal conditions in terms of the number of measurements needed to recover sparse signals [3-5]. However, the use of these matrices for real-world applications is limited for several reasons: no fast matrix multiplication algorithm is available, huge memory requirements for large scale problems, difficult implementation on hardware, etc. Let us consider s-sparse digital signals x ∈ ℂN in an orthonormal basis ψ = (ψ1,...,ψN) ∈ ℂN×N. The decomposition of x in this basis is denoted <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M1">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M76">View MathML</a> (·* denotes the conjugate transpose), and contains s non-zero entries. The original signal x is then probed by projection onto m randomly selected vectors of another orthonormal basis ϕ = (ϕ1,...,ϕN) ∈ ℂN×N. The indices Ω = {l1,...,lm} of the selected vectors are chosen independently and uniformly at random from {1,...,N}. We denote <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M2">View MathML</a> the m × N matrix made of the selected rows of ϕ*. The measurement vector y ∈ ℂm thus reads as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M3">View MathML</a>

(1)

We also denote A = ϕ* ψ ∈ ℂN×N. Finally we aim at recovering α by solving the ℓ1-minimization problem

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M4">View MathML</a>

(2)

where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M5">View MathML</a> (|·| denotes the complex magnitude). The reconstructed signal x* satisfies x* = ψα*.

The theory of compressed sensing already demonstrates that a small number m N of random measurements are sufficient for an accurate and stable reconstruction of x [6,7]. However, the recovery conditions depend on the mutual coherence μ between ϕ and ψ. This value is a similarity measure between the sensing and sparsity bases. It is defined as μ = max1 ≤ i,j N |〈ϕi,ψj〉| and satisfies N-1/2 μ ≤ 1. The performance is optimal when the bases are perfectly incoherent, i.e., μ = N-1/2, and unavoidably decreases when μ increases.

1.2 Fourier imaging applications and mutual coherence

The dependence of performance on the mutual coherence μ is a key concept in compressed sensing. It has significant implications for Fourier imaging applications, in particular radio interferometry or MRI, where signals are probed in the orthonormal Fourier basis. In radio interferometry, one of the main challenges is to reconstruct accurately the original signal from a limited number of accessible measurements [8-12]. In MRI, accelerating the acquisition process by reducing the number of measurements is of huge interest in, for example, static and dynamic imaging [13-17], parallel MRI [18-20], or MR spectroscopic imaging [21-23]. The theory of compressed sensing shows that Fourier acquisition is the best sampling strategy when signals are sparse in the Dirac basis. The sensing system is indeed optimally incoherent. Unfortunately, natural signals are usually rather sparse in multi-scale bases, e.g., wavelet bases, which are coherent with the Fourier basis. Many measurements are thus needed to reconstruct accurately the original signal. In the perspective of accessing better performance, sampling strategies that improve the incoherence of the sensing scheme should be considered.

1.3 Main contributions and organization

In the present study, we advocate a compressed sensing strategy coined spread spectrum that consists of a wide bandwidth pre-modulation of the signal x before projection onto randomly selected vectors of an orthonormal basis. In the particular case of Fourier measurements, the pre-modulation amounts to a convolution in the Fourier domain which spreads the power spectrum of the original signal x (hence the name "spread spectrum"), while preserving its norm. Equivalently, this spread spectrum phenomenon acts on each sparsity basis vector describing x so that information of each of them is accessible whatever the Fourier coefficient selected. This effect implies a decrease of coherence between the sparsity and sensing bases and enables an enhancement of the reconstruction quality.

In Section 2, we study the spread spectrum technique in a digital setting for arbitrary pairs of sensing and sparsity bases (ϕ, ψ). We consider a digital pre-modulation <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M6">View MathML</a> with |cl| = 1 and random phases identifying a random Rademacher or Steinhaus sequence. We show that the recovery conditions do not depend anymore on the coherence of the system but on a new parameter β (ϕ, ψ) called modulus-coherence and defined as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M7">View MathML</a>

(3)

where ϕki and ψkj are respectively the kth entries of the vectors ϕi and ψj. We then show that this parameter reaches its optimal value β (ϕ, ψ) = N-1/2 whatever the sparsity basis Ψ, for particular sensing matrices ϕ including the Fourier matrix, thus providing universal recovery performances. It is also efficient as sensing matrices with fast matrix multiplication algorithms can be used, thus reducing the need in memory requirement and computational power. In Section 3, these theoretical results are confirmed numerically through an analysis of the empirical phase transition of the ℓ1-minimization problem for different pairs of sensing and sparsity bases. In Section 4, we show that the spread spectrum technique remains effective in an analog setting with chirp modulation for application to realistic Fourier imaging, and illustrate these findings in the context of radio interferometry and MRI. Finally, we conclude in Section 5.

In the context of compressed sensing, the spread spectrum technique was already briefly introduced by the authors for compressive sampling of pulse trains in [24], applied to radio interferometry in [25,26] and to MRI in [27-30]. This article provides theoretical foundations for this technique, both in the digital and analog settings. Note that other acquisition strategies can be related to the spread spectrum technique as discussed in Section 2.5.

Let us also acknowledge that spread spectrum techniques are very popular in telecommunications. For example, one can cite the direct sequence spread spectrum (DSSS) and the frequency hopping spread spectrum (FHSS) techniques. The former is sometimes used over wireless local area networks, the latter is used in Bluetooth systems [31]. In general, spread spectrum techniques are used for their robustness to narrowband interference and also to establish secure communications.

2 Compressed sensing by spread spectrum

In this section, we first recall the standard recovery conditions of sparse signals randomly sampled in a bounded orthonormal system. These recovery results depend on the mutual coherence μ of the system. Hence, we study the effect of a random pre-modulation on this value and deduce recovery conditions for the spread spectrum technique. We finally show that the number of measurements needed to recover sparse signals becomes universal for a family of sensing matrices ϕ which includes the Fourier basis.

2.1 Recovery results in a bounded orthonormal system

For the setting presented in Section 1, the theory of compressed sensing already provides sufficient conditions on the number of measurements needed to recover the vector α from the measurements y by solving the ℓ1-minimization problem (2) [6,7].

Theorem 1 ([7], Theorem 4.4). Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M8">View MathML</a>be an s-sparse vector, Ω = {l1,...,lm} be a set of m indices chosen independently and uniformly at random from {1,...,N}, and y = Aα ∈ ℂm. For some universal constants C > 0 and γ > 1, if

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M9">View MathML</a>

(4)

then α is the unique minimizer of the 1-minimization problem (2) with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10">View MathML</a>.

Let us acknowledge that even if the measurements are corrupted by noise or if α is non-exactly sparse, the theory of compressed sensing also shows that the reconstruction obtained by solving the ℓ1-minimization problem remains accurate and stable:

Theorem 2 ([7], Theorem 4.4). Let A = ϕ*ψ, Ω = {l1,...,lm} be a set of m indices chosen independently and uniformly at random from {1,...,N}, and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M11">View MathML</a>be the best s-sparse approximation of the (possibly non-sparse) vector α ∈ ℂN. Let the noisy measurements y = AΩ α + n ∈ ℂm be given with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M12">View MathML</a>. For some universal constants D, E > 0 and γ > 1, if relation (4) holds, then the solution α* of the 1-minimization problem

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M13">View MathML</a>

(5)

satisfies

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M14">View MathML</a>

(6)

with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10">View MathML</a>.

In the above theorems, the role of the mutual coherence μ is crucial as the number of measurements needed to reconstruct x scales quadratically with its value. In the worst case where ϕ and ψ are identical, μ = 1 and the signal x is probed in a domain where it is also sparse. According to relation (4), the number of measurements necessary to recover x is of order N. This result is actually very intuitive. For an accurate reconstruction of signals sampled in their sparsity domain, all the non-zero entries need to be probed. It becomes highly probable when m N. On the contrary, when ϕ and ψ are as incoherent as possible, i.e., μ = N-1/2, the energy of the sparsity basis vectors spreads equally over the sensing basis vectors. Consequently, whatever the sensing basis vector selected, one always gets information of all the sparsity basis vectors describing the signal x, therefore reducing the need in the number of measurements. This is confirmed by relation (4) which shows that the number of measurements is of the order of s when μ = N-1/2. To achieve much better performance when the mutual coherence is not optimal, one would naturally try to modify the measurement process to achieve a better global incoherence. We will see in the following section that a simple random pre-modulation is an efficient way to achieve this goal whatever the sparsity matrix ψ.

2.2 Pre-modulation effect on the mutual coherence

The spread spectrum technique consists of pre-modulating the signal x by a wide-band signal <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M15">View MathML</a>, with |cl| = 1 and random phases, before projecting the resulting signal onto m vectors of the basis ϕ. The measurement vector y thus satisfies

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M16">View MathML</a>

(7)

where the additional matrix C ∈ ℝN×N stands for the diagonal matrix associated to the sequence c.

In this setting, the matrix Ac is orthonormal. Therefore, the recovery condition of sparse signals sampled with this matrix depends on the mutual coherence μ = max1 ≤ i,j N |〈ϕi, C ψj〉|. With a pre-modulation by a random Rademacher or Steinhaus sequence, Lemma 1 shows that the mutual coherence μ is essentially bounded by the modulus-coherence β (ϕ, ψ) defined in Equation (3).

Lemma 1. Let ϵ ∈ (0, 1), c ∈ ℂN be a random Rademacher or Steinhaus sequence and C ∈ ℂN×N be the associated diagonal matrix. Then, the mutual coherence μ = max1 ≤ i,j N |〈 ϕi, C ψj〉| satisfies

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M17">View MathML</a>

(8)

with probabilty at least 1 - ϵ.

The proof of Lemma 1 relies on a simple application of the Hoeffding's inequality and the union bound.

Proof. We have <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M18">View MathML</a>, where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M19">View MathML</a>. An application of the Hoeffding's inequality shows that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M20">View MathML</a>

for all u > 0 and 1 ≤ i, j N, with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M21">View MathML</a>. The union bound then yields

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M22">View MathML</a>

for all u > 0. As <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M23">View MathML</a> then <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M24">View MathML</a> for all 1 ≤ i, j N, and the previous relation becomes

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M25">View MathML</a>

for all u > 0.Taking <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M26">View MathML</a> terminates the proof.

2.3 Sparse recovery with the spread spectrum technique

Combining Theorem 1 with the previous estimate on the mutual coherence, we can state the following theorem:

Theorem 3. Let c ∈ ℂN, with N ≥ 2, be a random Rademacher or Steinhaus sequence, C ∈ ℂN×N be the diagonal matrix associated to c, α ∈ ℂN be an s-sparse vector, Ω = {l1,...,lm} be a set of m indices chosen independently and uniformly at random from {1,...,N}, and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M27">View MathML</a>, with Ac = ϕ*. For some constants 0 < ρ < log3(N) and Cρ > 0, if

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M28">View MathML</a>

(9)

then α is the unique minimizer of the 1-minimization problem (2) with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M29">View MathML</a>.

Proof. It is straightforward to check that C*C = CC* = I, where I is the identity matrix. The matrix Ac = ϕ*is thus orthonormal and Theorem 1 applies. To keep the notations simple, let F denotes the event of failure of the ℓ1-minimization problem (2), X be the event of m CNμ2 s log4 (N), and Y be the event of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M30">View MathML</a>. According to Theorem 1 and Lemma 1, the probability of F given X satisfies <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M31">View MathML</a> and the probability of Y satisfies <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M32">View MathML</a>.

We will see, at the end of this proof, that for a proper choice of ϵ, when condition (9) holds, we have

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M33">View MathML</a>

(10)

Using this fact, we compute the probability of failure <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M34">View MathML</a> of the ℓ1 minimization problem. We start by noticing that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M35">View MathML</a>

where Xc denotes the complement of event X. In the first inequality, the probability <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M75">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M36">View MathML</a> are saturated to 1. One can also note that if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M37">View MathML</a>, i.e., Y occurs, condition (10) implies that m CNμ2s log4 (N), i.e., X occurs. Therefore <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M38">View MathML</a> and

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M39">View MathML</a>

The probability of failure is thus bounded above by <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M40">View MathML</a>. Consequently if condition (10) holds with ϵ = Nand 0 < ρ < log3(N), α is the unique minimizer of the ℓ1-minimization problem (2) with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M41">View MathML</a>.

Finally, noticing that for ϵ = Nwith N ≥ 2, condition (10) always holds when condition (9), with Cρ = 2(3 + ρ)C, is satisfied, terminates the proof.

Note that relation (9) also ensures the stability of the spread spectrum technique relative to noise and compressibility by combination of Theorem 2 and Lemma 1.

2.4 Universal sensing bases with ideal modulus-coherence

Theorem 3 shows that the performance of the spread spectrum technique is driven by the modulus-coherence β(ϕ, ψ). In general the spread spectrum technique is not universal and the number of measurements required for accurate reconstructions depends on the value of this parameter.

Definition 1. (Universal sensing basis) An orthonormal basis ϕ ∈ ℂN×N is called a universal sensing basis if all its entries ϕki, 1 ≤ k,i N, are of equal complex magnitude.

For universal sensing bases, e.g., the Fourier transform or the Hadamard transform, we have |ϕki| = N-1/2 for all 1 ≤ k, i N. It follows that β (ϕ, ψ) = N-1/2 and μ N-1/2, i.e., its optimal value up to a logarithmic factor, whatever the sparsity matrix considered! For such sensing matrices, the spread spectrum technique is thus a simple and efficient way to render a system incoherent independently of the sparsity matrix.

Corollary 1. (Spread spectrum universality) Let c ∈ ℂN, with N ≥ 2, be a random Rademacher or Steinhaus sequence, C ∈ ℂN×N be the diagonal matrix associated to c, α ∈ ℂN be an s-sparse vector, Ω = {l1,...,lm} be a set of m indices chosen independently and uniformly at random from <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M42">View MathML</a>, with Ac = ϕ*. For some constants 0 < ρ < log3(N), Cρ > 0, and universal sensing bases ϕ ∈ ℂN×N, if

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M43">View MathML</a>

(11)

then α is the unique minimizer of the 1-minimization problem (2) with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M41">View MathML</a>.

For universal sensing bases, the spread spectrum technique is thus universal: the recovery condition does not depend on the sparsity basis and the number of measurements needed to reconstruct sparse signals is optimal in the sense that it is reduced to the sparsity level s. The technique is also efficient as the pre-modulation only requires a sample-by-sample multiplication between x and c. Furthermore, fast multiplication matrix algorithms are available for several universal sensing bases such as the Fourier or Hadamard bases.

In light of Corollary 1, one can notice that sampling sparse signals in the Fourier basis is a universal encoding strategy whatever the sparsity basis ψ - even if the original signal is itself sparse in the Fourier basis! We will confirm these results experimentally in Section 3.

2.5 Related work

Let us acknowledge that the techniques proposed in [32-37] can be related to the spread spectrum technique. The benefit of a random pre-modulation in the measurement system is already briefly suggested in [32]. The proofs of the claims presented in that conference paper have very recently been accepted for publication in [33] during the review process of this article. The authors obtain similar recovery results as those presented here. In [34], the author proposes to convolve the signal x with a random waveform and randomly under-sample the result in time-domain. The random convolution is performed through a random pre-modulation in the Fourier domain and the signal thus spreads in time-domain. In our setting, this method actually corresponds to taking ϕ as the Fourier matrix and ψ as the composition of the Fourier matrix and the initial sparsity matrix. In [35], the authors propose a technique to sample signals sparse in the Fourier domain. They first pre-modulate the signal by a random sequence, then apply a low-pass antialiasing filter, and finally sample it at low rate. Finally, random pre-modulation is also used in [36,37] but for dimension reduction and low dimensional embedding.

We recover similar results, albeit in a different way. We also have a more general interpretation. In particular, we proved that changing the sensing matrix from the Fourier basis to the Hadamard does not change the recovery condition (11).

3 Numerical simulations

In this section, we confirm our theoretical predictions by showing, through a numerical analysis of the phase transition of the ℓ1-minimization problem, that the spread spectrum technique is universal for the Fourier and Hadamard sensing bases.

3.1 Settings

For the first set of simulations, we consider the Dirac, Fourier, and Haar wavelet bases as sparsity basis ψ and choose the Fourier basis as the sensing matrix ϕ. We generate complex s-sparse signals of size N = 1,024 with s ∈ {1,...,N}. The positions of the non-zero coefficients are chosen uniformly at random in {1,...,N}, their phases are set by generating a Steinhaus sequence, and their amplitudes follow a uniform distribution over [0, 1]. The signals are then probed according to relation (1) or (7) and reconstructed from different number of measurements m ∈ {s,...,10s} by solving the ℓl-minimization problem (2) with the SPGL1 toolbox [38,39]. For each pair (m, s), we compute the probability of recoverya over 100 simulations.

For the second set of simulations, the same protocol is applied with the same sparsity basis but with the Hadamard basis as the sensing matrix ϕ.

3.2 Results

Figure 1 shows the phase transitions of the ℓ1-minimization problem obtained for sparse signals in the Dirac, Haar, and Fourier sparsity bases and probed in the Fourier basis with and without random pre-modulation. Figure 2 shows the same graphs but with measurements performed in the Hadamard basis. In the absence of pre-modulation, one can note that the phase transitions depend on the mutual coherence of the system as predicted by Theorem 1. For the pairs Fourier-Dirac and Hadamard-Dirac, the mutual coherence is optimal and the experimental phase transitions match the one of Donoho-Tanner (dashed green line) [5]. For all the other cases, the coherence is not optimal and the region where the signals are recovered is much smaller. The worst case is obtained for the pair Fourier-Fourier for which μ = 1. In the presence of pre-modulation, Corollary 1 predicts that the performance should not depend on the sparsity basis and should become optimal. It is confirmed by the phases transition showed on Figures 1 and 2 as they all match the phase transition of Donoho-Tanner, even for the pair Fourier-Fourier!

thumbnailFigure 1. Phase transition of the ℓ1-minimization problem for different sparsity bases and random selection of Fourier measurements without (left panels) and with (right panels) random modulation. The sparsity bases considered are the Dirac basis (top), the Haar wavelet basis (center), and the Fourier basis (bottom). The dashed green line indicates the phase transition of Donoho-Tanner [5]. The color bar goes from white to black indicating a probability of recovery from 0 to 1.

thumbnailFigure 2. Phase transition of the ℓ1-minimization problem for different sparsity bases and random selection of Hadamard measurements without (left panels) and with (right panels) random modulation. The sparsity bases considered are the Dirac basis (top), the Haar wavelet basis (center), and the Fourier basis (bottom). The dashed green line indicates the phase transition of Donoho-Tanner [5]. The color bar goes from white to black indicating a probability of recovery from 0 to 1.

4 Application to realistic Fourier imaging

In this section, we discuss the application of the spread spectrum technique to realistic analog Fourier imaging such as radio interferometric imaging or MRI. Firstly, we introduce the exact sensing matrix needed to account for the analog nature of the imaging problem. Secondly, while our original theoretical results strictly hold only in a digital setting, we derive explicit performance guarantees for the analog version of the spread spectrum technique. We also confirm on the basis of simulations that the spread spectrum technique drastically enhances the quality of reconstructed signals.

4.1 Sensing model

Radio interferometry dates back to more than 60 years ago [40-43]. It allows observations of the sky with angular resolutions and sensitivities inaccessible with a single telescope. In a few words, radio telescope arrays synthesize the aperture of a unique telescope whose size would be the maximum projected distance between two telescopes of the array on the plane perpendicular to line of sight. Considering small field of views, the signal probed can be considered as a planar image on the plane perpendicular to the pointing direction of the instrument. Measurements are obtained through correlation of the incoming electric fields between each pair of telescopes. As showed by the van Cittert-Zernike theorem [43], these measurements correspond to the Fourier transform of the image multiplied by an illumination function. In general, the number of spatial frequencies probed are much smaller than the number of frequencies required by the Nyquist-Shannon theorem, so that the Fourier coverage is incomplete. An ill-posed inverse problem is thus defined for reconstruction of the original image. To address this problem, approaches based on compressed sensing have recently been developed [10-12].

Magnetic resonance images are created by nuclear magnetic resonance in the tissues to be imaged. Standard MR measurements take the form of Fourier (also called k-space) coefficients of the image of interest. These measurements are obtained by application of linear gradient magnetic fields that provides the Fourier coefficient of the signal at a spatial frequency proportional the gradient strength and its duration of application. Accelerating the acquisition process, or equivalently increasing the achievable resolution for a fixed acquisition time, is of major interest for MRI applications. To address this problem, recent approaches based on compressed sensing seek to reconstruct the signal from incomplete information. In this context, several approaches have been designed [13,28-30,44-48].

In light of the results of Section 2, Fourier imaging is a perfect framework for the spread spectrum technique, apart from the analog nature of the corresponding imaging problems. In the quoted applications, the random pre-modulation is replaced by a linear chirp pre-modulation [25-30]. In radio interferometry, this modulation is inherently part of the acquisition process [25,26]. In MRI, it is easily implemented through the use of dedicated coils or RF pulses [29,30]. For two-dimensional signals, the linear chirp with chirp rate w ∈ ℝ reads as a complex-valued function <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M44">View MathML</a> of the spatial variable τ ∈ ℝ2. Note that for high chirp rates w, i.e., for chirp whose band-limit is of the same order of the band-limit of the signal, this chirp shares the following important properties with the random modulation: it is a wide-band signal which does not change the norm of the signal x, as |c(τ)| = 1 whatever τ ∈ ℝ2.

In this setting, the complete linear relationship between the signal and the measurements is given by

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M45">View MathML</a>

(12)

In the above equation, the matrix U represents an up-sampling operator needed to avoid any aliasing of the modulated signal due to a lack of sampling resolution in a digital description of the originally analog problem. The convolution in Fourier space induced by the analog modulation implies, in contrast with the digital setting studied before, that the band limit of the modulated signal is the sum of the individual band limits of the original signal and of the chirp c. We assume here that, on its finite field of view L, the signal x is approximately band-limited with a cut-off frequency at B, i.e., its energy beyond the frequency B is negligible. The signal x is thus discretized on a grid of N = 2LB points. On this field of view L, the linear chirp c may be approximated by a band limited function of band limit identified by its maximum instantaneous frequency |w|L/2. This band limit can also be parametrized in terms of a discrete chirp rate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M46">View MathML</a> and thus <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M47">View MathML</a>. Therefore, an up-sampled grid with at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M48">View MathML</a> points needs to be considered and the modulated signal is correctly obtained by applying the chirp modulation on the signal after up-sampling on the Nw points grid.b The up-sampling operator U, implemented in Fourier space by zero padding, is of size Nw × N and satisfies U*U = I ∈ ℂN×N. Finally, the matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M49">View MathML</a> is the diagonal matrix implementing the chirp modulation on this up-sampled grid and the matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M50">View MathML</a> stands for the discrete Fourier basis on the same grid. The indices Ω = {l,...,lm} of the Fourier vectors selected to probe the signal are chosen independently and uniformly at random from {1,...,Nw}.

4.2 Illustration

Up to the introduction of the matrix U and the substitution of the linear chirp modulation for the random modulation, we are in the same setting as the one studied in Section 2. To illustrate the effectiveness of the spread spectrum technique, we consider two images of size N = 256 × 256 showed in Figure 3. The first image shows the radio emission associated with the encounter of a galaxy with its northern neighbor. It was acquired with the very large array in New Mexico [49]. The second image shows a brain acquired in an MRI scanner. This image is part of the BRAINIX database [50]. These images are probed according to relation (12) in the absence (<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M51">View MathML</a>) and presence (<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M52">View MathML</a>) of a linear chirp modulation. Independent and identically distributed Gaussian noise with zero-mean is also added to the measurements. The variance σ2 of the noise is defined such that the input <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M53">View MathML</a> is 30 dB (Σx stands for the sampled standard deviation of x). The images are reconstructed from m = 0.4 N complex Fourier measurements by solving the ℓ1-minimization problem (5). Note that in order to stay in the setting of the theorems presented so far, no reality constraint is enforced in the reconstructions, so that the reconstructed images are complex valued. The sparsity bases ψ used are the Daubechies-6 and Haar wavelet bases for the galaxy and the brain, respectively. In each case, 20 reconstructions are performed for different noise and mask realizations. The complex magnitudes of reconstructed images with median mean squared errors are presented in Figure 3.

thumbnailFigure 3. Top panels: Image of the giant elliptical galaxy NGC1316 (center of the image) devouring its small northern neighbor. The image shows the radio emission associated with this encounter superimposed on an optical image. The radio emission was imaged using the very large array in New Mexico (Image courtesy of NRAO/AUI and Uson). The image size is N = 256 × 256. Bottom panels: MRI image of a brain from the BRAINIX database. From left to right: original image; complex magnitudes of the reconstructed images from m = 0.4N measurements without chirp modulation; complex magnitudes of the reconstructed images from m = 0.4N measurements with chirp modulation.

In the absence of linear chirp modulation, the quality of the reconstructed image is very low. However, one can already note that the fine scale structures are much better reconstructed than those at large scales. The fine details live at the small scales of the wavelet decomposition whereas the large structures live at larger scales. The small scale wavelets being more incoherent with the Fourier basis than the larger wavelets, the high frequency details are naturally better recovered.

In the presence of the linear chirp modulation, all the wavelets in ψ become optimally incoherent with the Fourier basis thanks to the universality of the spread spectrum technique. Consequently, as one can observe on Figure 3, the low and high frequency details are better reconstructed and the image quality is drastically enhanced. Note that much better reconstructions can be obtained for the brain image by substituting the total variation normc for the ℓ1 norm in (5) [13,30]. However, Theorems 1 and 3 do not hold for such a norm.

Let us acknowledge that these simulations are not fully realistic. For example, in radio-interferometry the spatial frequency cannot be chosen at random. To simulate realistic acquisitions, one would have to consider non-random measurements in the continuous Fourier plane. Such a study is beyond the scope of this study. However, in the context of MRI, part of the authors implemented and tested this technique on a real scanner with in vivo acquisitions [29,30].

4.3 Modified recovery condition

Because of the modifications introduced in (12) to account for the analog nature of the problem, the digital theory associated with the measurement matrix (7) does not explicitly apply. Nevertheless, the previous illustration shows that the spread spectrum technique is indeed still very effective in this analog setting. Actually, performance guarantees similar to Theorem 1 may be obtained in this setting.

Theorem 4. Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M54">View MathML</a>be an s-sparse vector, Ω = {l1,...,lm} be a set of m indices chosen independently and uniformly at random from {1,...,Nw}, and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M55">View MathML</a>. For some universal constants C > 0 and γ > 1, if the number of measurements m satisfies

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M56">View MathML</a>

(13)

then α is the unique minimizer of the 1-minimization problem (2) with probability at least <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M10">View MathML</a>.

Proof. The proof follows directly from Theorem 4.4 in [7]. Indeed, Theorem 4.4 applies to any matrices AΩ associated to an orthonormal system (with respect to the probability measure used to draw Ω) that satisfies the so-called boundedness condition (see Section 4.1 of [7] for more details).

Let us denote <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M57">View MathML</a> the normalized measurement matrix. It is easy to check that this new matrix is orthonormal relative to the discrete uniform probability measure on {1,...,Nw}. Indeed

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M58">View MathML</a>

as U*U = I ∈ ℂN×N and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M59">View MathML</a>. Furthermore, the boundedness condition is satisfied as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M60">View MathML</a>. Applying Theorem 4.4 in [7] to the matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M61">View MathML</a> terminates the proof.

Note that Theorem 4.4 in [7] also ensures that our analog sensing scheme is stable relative to noise and non exact sparsity if condition (13) is satisfied. Also note that one can obtain a similar results using Theorems 1.1 and 1.2 of the very recently accepted article [51].

In view of this theorem, one can notice that the number of measurements needed for accurate reconstructions of sparse signals is proportional to the sparsity s times the product <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a>, with factors depending on the chirp rate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a>. In the analog framework, two effects are actually competing. On the one hand, the mutual coherence <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M64">View MathML</a> of the system is decreasing with the spread spectrum phenomenon, but, on the other hand, the number of accessible frequencies Nw that bear information is increasing linearly with the chirp rate. The optimal recovery conditions are reached for a chirp rate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a> that ensures that the product <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a> is at its minimum. This minimum might change depending on the sparsity matrix, so the universality of the recovery is formally lost.

To illustrate this effect, Table 1 shows values of the product <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a> for different chirp rates <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a> and two different sparsity matrices: the Fourier and the Dirac bases. The values are computed numerically for a size of signal N = 1,024. In the case of the Dirac basis, one can notice that the product <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a> slightly increases with the chirp rate, thus predicting that the performance should even slightly diminish in the presence of a chirp. On the contrary, the product is drastically reduced for the Fourier basis as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a> increases, predicting a significant performance improvement in the presence of chirp modulation.

Table 1. Influence of a chirp modulation on <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a>

4.4 Experiments

To confirm the theoretical predictions of the previous section, we consider the Dirac and Fourier bases as sparsity matrices ψ. We then generate complex s-sparse signals of size N = 1, 024 with s = 10. The positions of the non-zero coefficients are chosen uniformly at random in {1,..., N}, their signs are set by generating a Steinhaus sequence, and their amplitudes follow a uniform distribution in [0, 1]. The signals are then probed according to relation (12) and reconstructed from different number of measurements m ∈ {s,..., N} by solving the ℓ1-minimization problem (2) with the SPGL1 toolbox. For each pair (m, s), we compute the probability of recovery over 100 simulations for different chirp rate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M69">View MathML</a>.

Figure 4 shows the probability of recovery ϵ as a function of the number of measurements.

thumbnailFigure 4. Probability of recovery ϵ of 10-sparse signals as a function of the number measurement m obtained with the measurement matrix (12) for two different sparsity basis: the Dirac basis (left) and the Fourier basis (right). The continuous black curve corresponds to the probability of recovery for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M70">View MathML</a>. The dot-dashed blue curve corresponds to the probability of recovery for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M71">View MathML</a>. The dashed black curve corresponds to the probability of recovery for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M72">View MathML</a>. The continuous red curve corresponds to the probability of recovery for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M73">View MathML</a>.

First, in the case where ψ is the Dirac basis, one can notice that the number of measurements needed to reach a probability of recovery of 1 slightly increases with the chirp rate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a>. This is in line with the value in Table 1 and Theorem 4.

Second, in the case where ψ is the Fourier basis, the performance becomes much better in the presence of a chirp. As predicted by the value in Table 1, the improvement is drastic when <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M63">View MathML</a> goes from 0 to 0.1 and then starts to saturate between 0.1 and 0.5.

Third, according to Table 1, the product <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M62">View MathML</a> is equal to 4.13 for the Fourier basis when <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M73">View MathML</a>. This is nearly the value obtained with the Dirac basis for the same chirp rate, suggesting the same probability of recovery for the same number of measurements. Indeed, one can notice on Figure 4 that the number of measurement needed to reach a probability of recovery of 1 is around 100 in both cases.

Finally, these results also suggest that the spread spectrum technique in the modified setting is almost universal in practice. Indeed, for the perfectly incoherent pair Fourier-Dirac of sensing-sparsity bases, the number of measurements needed for perfect recovery is around 100 and this number remains almost unchanged in presence of the linear chirp modulation. Furthermore, for the pair Fourier-Fourier, the spread spectrum technique allows to reduce the number of measurements for perfect recovery close to this optimal value.

5 Conclusion

We have presented a compressed sensing strategy that consists of a wide bandwidth pre-modulation of the signal of interest before projection onto randomly selected vectors of an orthonormal basis. In a digital setting with a random pre-modulation, the technique was proved to be universal for sensing bases such as the Fourier or Hadamard bases, where it may be implemented efficiently. Our results were confirmed through a numerical analysis of the phase transition of the ℓ1-minimization problem for different pairs of sensing and sparsity bases.

The spread spectrum technique was also shown to be of great interest for realistic analog Fourier imaging. In applications such as radio interferometry and MRI, the originally digital random pre-modulation may be mimicked by an analog linear chirp. Explicit performance guarantees for the analog version of the technique with a chirp modulation were derived. It shows that recovery results are still enhanced in this setting, though universality does not strictly hold anymore. Numerical simulations have shown that the quality of reconstructed signals is drastically enhanced in this more realistic setting, also for pairs of sensing-sparsity bases initially highly coherent, such as the Fourier-Fourier pair.

Competing interests

Part of this study was funded by Merck Serono S.A.

Authors' contributions

GP carried out the theoretical study, the numerical experiments, and wrote most of the manuscript. PV and RG conceived of the study, participated in its design and coordination, and guided GP in the theoretical study. YW contributed to the setup of the technique for realistic Fourier imaging, participated in the design of the associated numerical experiments and in the writing of the manuscript. All authors discussed the results, read and approved the final manuscript.

Endnotes

aperfect recovery is considered if the ℓ2 norm between the original signal x and the reconstructed signal x* satisfy: ||x - x*||2 ≤ 10-3||x||2.

bIn a full generality, natural signals are not necessarily band-limited. The spread spectrum technique can easily be adapted to this case. The sensing model should simply be modified to account for the fact that, if measurements are performed at frequencies up to a band limit B, they unavoidably contain energy of the signal up to band limit <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/6/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/6/mathml/M74">View MathML</a>.

c1norm of the magnitude of the gradient.

Acknowledgements

This study was supported in part by the Center for Biomedical Imaging (CIBM) of the Geneva and Lausanne Universities, EPFL, and the Leenaards and Louis-Jeantet foundations, in part by the Swiss National Science Foundation (SNSF) under grant PP00P2-123438, also by the EU FET-Open project FP7-ICT-225913-SMALL: Sparse Models, Algorithms and Learning for Large-Scale data, and by the EPFL-Merck Serono Alliance award.

References

  1. E Candes, J Romberg, T Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 52, 489–509 (2006)

  2. DL Donoho, Compressed sensing. IEEE Trans Inf Theory 52, 1289–1306 (2006)

  3. R Baraniuk, Compressive sensing. IEEE Signal Process. Mag 24, 118–121 (2007)

  4. E Candes, T Tao, Near optimal signal recovery from random projections: universal encoding strategies. IEEE Trans Inf Theory 52, 5406–5425 (2006)

  5. DL Donoho, J Tanner, Counting faces of randomly-projected polytopes when the projection radically lowers dimension. J Am Math Soc 22, 1–53 (2009)

  6. E Candes, J Romberg, Sparsity and incoherence in compressive sampling. Inverse Probl 23, 969–985 (2007). Publisher Full Text OpenURL

  7. H Rauhut, Compressive Sensing and Structured Random Matrices. Theoretical Foundations and Numerical Methods for Sparse Recovery, [Radon Series on Computational and Applied Mathematics] (De Gruyter, Berlin, 2010) 9, pp. 1–92

  8. JA Högbom, Aperture synthesis with a non-regular distribution of interferometer baselines. Astron Astrophys 15, 417–426 (1974)

  9. TJ Cornwell, KF Evans, A simple maximum entropy deconvolution algorithm. Astron Astrophys 143, 77–83 (1985)

  10. Y Wiaux, L Jacques, G Puy, AMM Scaife, P Vandergheynst, Compressed sensing imaging techniques for radio interferometry. Mon Not R Astron Soc 395, 1733–1742 (2009). Publisher Full Text OpenURL

  11. S Wenger, M Magnor, Y Pihlström, S Bhatnagar, U Rau, SparseRI: A compressed sensing framework for aperture synthesis imaging in radio astronomy. Publ Astron Soc Pac 122, 1397–1374 (2010). Publisher Full Text OpenURL

  12. F Li, TJ Cornwell, F de Hoog, The application of compressive sampling to radio astronomy I: Deconvolution. Astron Astrophys 528, A31 (2011)

  13. M Lustig, D Donoho, JM Pauly, Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 58, 1182–1195 (2007). PubMed Abstract | Publisher Full Text OpenURL

  14. H Jung, JC Ye, EY Kim, Improved k-t BLAST and k-t SENSE using FOCUSS. Phys Med Biol 52, 3201–3226 (2007). PubMed Abstract | Publisher Full Text OpenURL

  15. U Gamper, P Boesiger, S Kozerke, Compressed sensing in dynamic MRI. Magn Reson Med 59, 365–373 (2008). PubMed Abstract | Publisher Full Text OpenURL

  16. H Jung, K Sung, K Nayak, EY Kim, JC Ye, k-t FOCUSS: A general compressed sensing framework for high resolution dynamic MRI. Magn Reson Med 61, 103–116 (2009). PubMed Abstract | Publisher Full Text OpenURL

  17. M Usman, C Prieto, T Schaeffter, P Batchelor, k-t group sparse: A method for accelerating dynamic MRI. Magn Reson Med 66, 1163–1176 (2011). PubMed Abstract | Publisher Full Text OpenURL

  18. D Liang, B Liu, J Wang, L Ying, Accelerating sense using compressed sensing. Magn Reson Med 62, 1574–1584 (2009). PubMed Abstract | Publisher Full Text OpenURL

  19. M Lustig, J Pauly, Spirit: Iterative self-consistent parallel imaging reconstruction from arbitrary k-space. Magn Reson Med 64, 457–471 (2010). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. R Otazo, D Kim, L Axel, D Sodickson, Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI. Magn Reson Med 64, 767–776 (2010). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. S Hu, M Lustig, AP Chen, J Crane, AK A, DA Kelley, R Hurd, J Kurhanewicz, SJ Nelson, JM Pauly, D Vigneron, Compressed sensing for resolution enhancement of hyperpolarized 13C flyback 3D-MRSI. J Magn Reson 192, 258–264 (2008). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. S Hu, M Lustig, AP Chen, A Balakrishnan, PEZ Larson, R Bok, J Kurhanewicz, SJ Nelson, A Goga, JM Pauly, DB Vigneron, 3D compressed sensing for highly accelerated hyperpolarized 13C MRSI with in vivo applications to transgenic mouse models of cancer. Magn Reson Med 63, 312–321 (2010). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  23. PEZ Larson, S Hu, M Lustig, AB Kerr, SJ Nelson, J Kurhanewicz, JM Pauly, DB Vigneron, Fast dynamic 3D MR spectroscopic imaging with compressed sensing and multiband excitation pulses for hyperpolarized 13C studies. Magn Reson Med 65, 610–619 (2011). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. FM Naini, R Gribonval, L Jacques, P Vandergheynst, Compressive sampling of pulse trains: spread the spectrum! In. Proc IEEE Int Conf Acoustics Speech and Signal Process (Taipei, Taiwan, 2009), pp. 2877–2880

  25. Y Wiaux, G Puy, Y Boursier, P Vandergheynst, Compressed sensing for radio interferometry: spread spectrum imaging techniques. Proc SPIE Conf WAVELET XIII (San Diego, CA, USA, 2009) 7446, p. 74460J

  26. Y Wiaux, G Puy, Y Boursier, P Vandergheynst, Spread spectrum for imaging techniques in radio interferometry. Mon Not R Astron Soc 400, 1029–1038 (2009). Publisher Full Text OpenURL

  27. G Puy, Y Wiaux, R Gruetter, JP Thiran, DV de Ville, P Vandergheynst, Spread spectrum for interferometric and magnetic resonance imaging. Proc IEEE Int Conf on Acoustics Speech and Signal Process (Dallas, TX, USA, 2010), pp. 2802–2805

  28. Y Wiaux, G Puy, R Gruetter, JP Thiran, DV de Ville, P Vandergheynst, Spread spectrum for compressed sensing techniques in magnetic resonance imaging. Proc IEEE Int Sym on Biomed Imaging (Rotterdam, Netherlands, 2010), pp. 756–759

  29. G Puy, J Marques, R Gruetter, JP Thiran, DV de Ville, P Vandergheynst, Y Wiaux, Accelerated MR imaging with spread spectrum encoding. Proc Intl Soc Mag Reson Med (Montreal, Canada, 2011) 19, p. 2808

  30. G Puy, J Marques, R Gruetter, JP Thiran, DV de Ville, P Vandergheynst, Y Wiaux, Spread spectrum magnetic resonance imaging. IEEE Trans Med Image (2012)

  31. J Schiller, Mobile Communications, 2nd edn. (Addison-Wesley, London, 2003)

  32. T Do, T Tran, L Gan, Fast compressive sampling with structurally random matrices. Proc IEEE Int Conf on Acoustics Speech and Signal Process (Las Vegas, NV, USA, 2008), pp. 3369–3372

  33. T Do, L Gan, N Nguyen, T Tran, Fast and efficient compressive sensing using structurally random matrices. IEEE Trans Signal Process 60, 139–154 (2012)

  34. J Romberg, Compressive sensing by random convolution. SIAM J Imag Sci 2, 1098–1128 (2009). Publisher Full Text OpenURL

  35. JA Tropp, JN Laska, MF Duarte, JK Romberg, RG Barianuk, Beyond nyquist: efficient sampling of sparse bandlimited signals. IEEE Trans Inf Theory 56, 520–544 (2010)

  36. F Krahmer, R Ward, New and improved Johnson-Lindenstrauss embeddings via the restricted isometry property. SIAM J Math Anal 43, 1269–1281 (2011). Publisher Full Text OpenURL

  37. JA Tropp, Improved analysis of the subsampled randomized Hadamard transform. Adv Adapt Data Anal 3(1-2), 115–126 (2012)

  38. E van den Berg, MP Friedlander, Probing the Pareto frontier for basis pursuit solutions. SIAM J Sci Comput 31, 890–912 (2008)

  39. SGPL1:, A solver for large-scale sparse reconstruction [http://www.cs.ubc.ca/labs/scl/spgl1/] webcite

  40. M Ryle, DD Vonberg, Solar Radiation on 175 Mc./s. Nature 1958, 339–340

  41. M Ryle, A Hewish, J Shakeshaft, The synthesis of large radio telescopes by the use of radio interferometers. IRE Trans Antennas Propag 7, 120–124 (1959)

  42. M Ryle, A Hewish, The Synthesis of large radio telescopes. Mon Not R Astron Soc 120, 220–230 (1960)

  43. AR Thompson, JM Moran, GW Swenson, Interferometry and Synthesis in Radio Astronomy, 2nd edn. (Wiley-Interscience, New York, 2001)

  44. F Sebert, YM Zou, L Ying, Compressed sensing MRI with random B1 field. Proc ISMRM (Toronto, Canada, 2008), p. 3151

  45. D Liang, G Xu, H Wang, KF King, D Xu, L Ying, Toeplitz random encoding MR imaging using compressed sensing. Proc IEEE Int Symp Biomed Imag From Nano to Macro (ISBI) (Boston, MA, USA, 2009), pp. 270–273

  46. H Wang, D Liang, K King, L Ying, Toeplitz random encoding for reduced acquisition using compressed sensing. Proc ISMRM (Honolulu, Hawai, 2009), p. 2669

  47. M Seeger, H Nickisch, R Pohmann, B Scholkopf, Optimization of k-space trajectories for compressed sensing by Bayesian experimental design. Magn Reson Med 63, 116–126 (2010). PubMed Abstract | Publisher Full Text OpenURL

  48. J Haldar, D Hernando, ZP Liang, Compressed-sensing in MRI with random encoding. IEEE Trans Med Imag 30, 898–903 (2011)

  49. National Radio Astronomy Observatory, [http://www.vla.nrao.edu/] webcite

  50. DICOM sample image sets, [http://pubimage.hcuge.ch:8080/] webcite

  51. E Candes, Y Plan, A Probabilistic and RIPless Theory of Compressed Sensing. IEEE Trans Inf Theory 52, 1289–1306 (2006)