[lecture NOTES]
Compressive Sensing
T he Shannon/Nyquist sam-
pling theorem specifies that
to avoid losing information
when capturing a signal, one
must sample at least two
times faster than the signal bandwidth. In
many applications, including digital
image and video cameras, the Nyquist rate
is so high that too many samples result,
making compression a necessity prior to
storage or transmission. In other applica-
tions, including imaging systems (medical
scanners and radars) and high-speed ana-
log-to-digital converters, increasing the
sampling rate is very expensive.
This lecture note presents a new
method to capture and represent com-
pressible signals at a rate significantly
below the Nyquist rate. This method,
called compressive sensing, employs
nonadaptive linear projections that pre-
serve the structure of the signal; the sig-
nal is then reconstructed from these
projections using an optimization
process [1], [2].
RELEVANCE
The ideas presented here can be used to
illustrate the links between data acquisi-
tion, compression, dimensionality reduc-
tion, and optimization in undergraduate
and graduate digital signal processing,
statistics, and applied mathematics
courses.
PREREQUISITES
The prerequisites for understanding this
lecture note material are linear algebra,
basic optimization, and basic probability.
PROBLEM STATEMENT
COMPRESSIBLE SIGNALS
Consider a real-valued, finite-length,
one-dimensional, discrete-time signal x,
Richard G. Baraniuk
which can be viewed as an N × 1 column
vector in RN with elements x[n],
n = 1, 2, . . . , N. (We treat an image or
higher-dimensional data by vectorizing it
into a long one-dimensional vector.) Any
signal in RN can be represented in terms
of a basis of N × 1 vectors {ψi}N
i=1. For
simplicity, assume that the basis is
orthonormal. Using the N × N basis
matrix = [ψ1|ψ2| . . .|ψN] with the
vectors {ψi} as columns, a signal x can
be expressed as
x = N
i=1
si ψi
or
x = ψs
(1)
where s is the N × 1 column vector of
weighting coefficients si = x, ψi = ψ T
i x
and ·T denotes transposition. Clearly, x
and s are equivalent representations of
the signal, with x in the time or space
domain and s in the domain.
The signal x is K-sparse if it is a linear
combination of only K basis vectors; that
is, only K of the si coefficients in (1) are
nonzero and (N − K) are zero. The case
of interest is when K N. The signal x
is compressible if the representation (1)
has just a few large coefficients and many
small coefficients.
TRANSFORM CODING
AND ITS INEFFICIENCIES
The fact that compressible signals are
well approximated by K-sparse represen-
tations forms the foundation of trans-
form coding [3]. In data acquisition
systems (for example, digital cameras)
transform coding plays a central role: the
full N-sample signal x is acquired; the
complete set of transform coefficients
{si} is computed via s = Tx; the K
largest coefficients are located and the
(N − K) smallest coefficients are dis-
carded; and the K values and locations of
the largest coefficients are encoded.
Unfortunately, this sample-then-com-
press framework suffers from three
inherent inefficiencies. First, the initial
number of samples N may be large even
if the desired K is small. Second, the set
of all N transform coefficients {si} must
be computed even though all but K of
them will be discarded. Third, the loca-
tions of the large coefficients must be
encoded, thus introducing an overhead.
THE COMPRESSIVE
SENSING PROBLEM
Compressive sensing address these ineffi-
ciencies by directly acquiring a com-
pressed signal representation without
going through the intermediate stage of
acquiring N samples [1], [2]. Consider a
general linear measurement process that
computes M < N inner products
between x and a collection of vectors
{φj}M
j=1 as in yj = x, φj. Arrange the
measurements yj in an M × 1 vector y
and the measurement vectors φ T
j as rows
in an M × N matrix . Then, by substi-
tuting from (1), y can be written as
y = x = s = s
(2)
where = is an M × N matrix. The
measurement process is not adaptive,
meaning that is fixed and does not
depend on the signal x. The problem
consists of designing a) a stable meas-
urement matrix such that the salient
information in any K-sparse or com-
pressible signal is not damaged by the
dimensionality reduction from x ∈ RN
to y ∈ RM and b) a reconstruction algo-
rithm to recover x from only M ≈ K
measurements y (or about as many
measurements as the number of coeffi-
cients recorded by a traditional trans-
form coder).
IEEE SIGNAL PROCESSING MAGAZINE [118]
IEEE SIGNAL PROCESSING MAGAZINE [118]
JULY 2007
JULY 2007
1053-5888/07/$25.00©2007IEEE
y
Φ
Ψ
S
M
=
Θ
S
y
=
N
K-sparse
x
(a)
(b)
[FIG1] (a) Compressive sensing measurement process with a random Gaussian
measurement matrix and discrete cosine transform (DCT) matrix . The vector of
coefficients s is sparse with K = 4. (b) Measurement process with = . There are
four columns that correspond to nonzero si coefficients; the measurement vector y is a
linear combination of these columns.
SOLUTION
DESIGNING A STABLE
MEASUREMENT MATRIX
The measurement matrix must allow
the reconstruction of the length-N signal
x from M < N measurements (the vector
y). Since M < N, this problem appears
ill-conditioned. If, however, x is K-sparse
and the K locations of the nonzero coef-
ficients in s are known, then the problem
can be solved provided M ≥ K. A neces-
sary and sufficient condition for this sim-
plified problem to be well conditioned is
that, for any vector v sharing the same K
nonzero entries as s and for some > 0
1 − ≤ v2
v2
≤ 1 + .
(3)
That is, the matrix must preserve the
lengths of these particular K-sparse vec-
tors. Of course, in general the locations
of the K nonzero entries in s are not
known. However, a sufficient condition
for a stable solution for both K-sparse
and compressible signals is that satis-
fies (3) for an arbitrary 3K-sparse vector
v. This condition is referred to as the
restricted isometry property (RIP) [1]. A
related condition, referred to as incoher-
ence, requires that the rows {φj} of
cannot sparsely represent the columns
{ψi} of (and vice versa).
Direct construction of a measure-
ment matrix such that = has
the RIP requires verifying (3) for each
of the
possible combinations of K
nonzero entries in the vector v of
N
K
length N. However, both the RIP and
incoherence can be achieved with high
probability simply by selecting as a
random matrix. For instance, let the
matrix elements φj,i be independent
and identically distributed (iid) random
variables from a Gaussian probability
density function with mean zero and
variance 1/N [1], [2], [4]. Then the
measurements y are merely M different
randomly weighted linear combinations
of the elements of x, as illustrated in
Figure 1(a). The Gaussian measure-
ment matrix has two interesting and
useful properties:
■ The matrix is incoherent with
the basis = I of delta spikes with
high probability. More specifically, an
M × N
iid Gaussian matrix
= I = can be shown to have
the RIP with high probability if
M ≥ cK log(N/K), with c a small
constant [1], [2], [4]. Therefore, K-
sparse and compressible signals of
length N can be recovered from
only M ≥ cK log(N/K) N random
Gaussian measurements.
■ The matrix is universal in the
sense that = will be iid
Gaussian and thus have the RIP with
high probability regardless of the
choice of orthonormal basis .
DESIGNING A SIGNAL
RECONSTRUCTION ALGORITHM
The signal reconstruction algorithm
must take the M measurements in the
vector y, the random measurement
matrix (or the random seed that gen-
IEEE SIGNAL PROCESSING MAGAZINE [119]
JULY 2007
erated it), and the basis and recon-
struct the length-N signal x or, equiva-
lently, its sparse coefficient vector s. For
K-sparse signals, since M < N in (2)
there are infinitely many s that satisfy
s = y. This is because if s = y then
(s + r) = y for any vector r in the null
space N () of . Therefore, the signal
reconstruction algorithm aims to find
the signal’s sparse coefficient vector in
the (N − M)-dimensional translated null
space H = N () + s.
(sp)p =
■ Minimum 2 norm reconstruction:
Define the p norm of the vector s as
i=1 |si|p. The classical
approach to inverse problems of this
type is to find the vector in the trans-
lated null space with the smallest 2
norm (energy) by solving
s = argmins
2 such that s
N
= y.
(4)
This optimization has the convenient
closed-form solution s = T(T)−1 y.
Unfortunately, 2 minimization will
almost never find a K-sparse solution,
returning instead a nonsparse s with
many nonzero elements.
■ Minimum 0 norm reconstruction:
Since the 2 norm measures signal
energy and not signal sparsity, con-
sider the 0 norm that counts the
number of non-zero entries in s.
(Hence a K-sparse vector has 0
norm equal to K.) The modified opti-
mization
s = argmins
0 such that s
= y
(5)
can recover a K-sparse signal exactly
with high probability using only
M = K + 1 iid Gaussian measure-
ments [5]. Unfortunately, solving (5)
is both numerically unstable and NP-
complete, requiring an exhaustive
enumeration of all
possible loca-
tions of the nonzero entries in s.
■ Minimum 1 norm reconstruction:
Surprisingly, optimization based on
the 1 norm
s = argmins
1 such that s
N
K
= y
(6)
[lecture NOTES] continued
can exactly recover K-sparse signals and
closely approximate compressible signals
with high probability using only
M ≥ cK log(N/K) iid Gaussian meas-
urements [1], [2]. This is a convex opti-
mization problem that conveniently
reduces to a linear program known as
basis pursuit [1], [2] whose computation-
al complexity is about O(N 3). Other,
related reconstruction algorithms are
proposed in [6] and [7].
DISCUSSION
The geometry of the compressive sensing
problem in RN helps visualize why 2
reconstruction fails to find the sparse
solution that can be identified by 1
reconstruction. The set of all K-sparse
vectors s in RN is a highly nonlinear
space consisting of all K-dimensional
hyperplanes that are aligned with the
coordinate axes as shown in Figure 2(a).
The translated null space H = N () + s
is oriented at a random angle due to the
randomness in the matrix as shown in
Figure 2(b). (In practice N, M, K 3, so
any intuition based on three dimensions
may be misleading.) The 2 minimizer s
from (4) is the point on H closest to the
origin. This point can be found by blow-
ing up a hypersphere (the 2 ball) until it
contacts H. Due to the random orienta-
tion of H, the closest point s will live
away from the coordinate axes with high
probability and hence will be neither
sparse nor close to the correct answer s.
In contrast, the 1 ball in Figure 2(c) has
points aligned with the coordinate axes.
Therefore, when the 1 ball is blown up,
it will first contact the translated null
space H at a point near the coordinate
axes, which is precisely where the sparse
vector s is located.
While the focus here has been on dis-
crete-time signals x, compressive sensing
also applies to sparse or compressible
analog signals x(t) that can be represent-
ed or approximated using only K out of
N possible elements from a continuous
basis or dictionary {ψi(t)}N
i=1 . While
each ψi(t ) may have large bandwidth
(and thus a high Nyquist rate), the signal
x(t ) has only K degrees of freedom and
thus can be measured at a much lower
rate [8], [9].
H
S
S
(b)
H
S
S
(c)
S
(a)
[FIG2] (a) The subspaces containing two sparse vectors in R3 lie close to the
coordinate axes. (b) Visualization of the 2 minimization (5) that finds the non-
sparse point-of-contact s between the 2 ball (hypersphere, in red) and the
translated measurement matrix null space (in green). (c) Visualization of the 1
minimization solution that finds the sparse point-of-contact s with high probability
thanks to the pointiness of the 1 ball.
Scene
Photodiode
Bitstream
A/D
Reconstruction
Image
DMD
Array
RNG
(a)
(b)
(c)
[FIG3] (a) Single-pixel, compressive sensing camera. (b) Conventional digital camera
image of a soccer ball. (c) 64 × 64 black-and-white image x of the same ball (N = 4,096
pixels) recovered from M = 1,600 random measurements taken by the camera in (a).
The images in (b) and (c) are not meant to be aligned.
PRACTICAL EXAMPLE
As a practical example, consider a sin-
gle-pixel, compressive digital camera
that directly acquires M random linear
measurements without first collecting
the N pixel values [10]. As illustrated in
Figure 3(a), the incident light-field cor-
responding to the desired image x is
reflected off a digital micromirror device
(DMD) consisting of an array of N tiny
mirrors. (DMDs are present in many
computer projectors and projection tele-
visions.) The reflected light is then col-
lected by a second lens and focused onto
a single photodiode (the single pixel).
Each mirror can be independently ori-
ented either towards the photodiode
(corresponding to a 1) or away from the
photodiode (corresponding to a 0). To
collect measurements, a random number
generator (RNG) sets the mirror orienta-
tions in a pseudorandom 1/0 pattern to
create the measurement vector φj. The
voltage at the photodiode then equals yj,
which is the inner product between φj
and the desired image x. The process is
repeated M times to obtain all of the
entries in y.
(continued on page 124)
IEEE SIGNAL PROCESSING MAGAZINE [120]
JULY 2007
[dsp TIPS&TRICKS] continued
line structure in the spectrum of the
recursive CORDIC, and the phase error
correction is not applied to suppress phase
error artifacts but rather to complete the
phase rotation left incomplete due to the
residual phase term in the angle accumu-
lator. This is a very different DDS!
IMPLEMENTATION
As a practical note, there are truncating
quantizers between the AGC multipliers
and the feedback delay element regis-
ters. As such, the truncation error circu-
lates in the registers and contributes an
undesired dc component to the complex
sinusoid output. This dc component can
(and should) be suppressed by using a
sigma delta-based dc cancellation loop
between the AGC multipliers and the
feedback delay elements [6].
CONCLUSIONS
We modified the traditional recursive
DDS complex oscillator structure to a
[lecture NOTES] continued from page 120
An image acquired with the single-
pixel camera using about 60% fewer ran-
dom measurements than reconstructed
pixels is illustrated in Figure 3(c); com-
pare to the target image in Figure 3(b).
The reconstruction was performed via a
total variation optimization [1], which is
closely related to the 1 reconstruction in
the wavelet domain. In addition to requir-
ing fewer measurements, this camera can
image at wavelengths where is difficult or
expensive to create a large array of sen-
sors. It can also acquire data over time to
enable video reconstruction [10].
CONCLUSIONS:
WHAT WE HAVE LEARNED
Signal acquisition based on compressive
sensing can be more efficient than tradi-
tional sampling for sparse or compressible
signals. In compressive sensing, the famil-
iar least squares optimization is inadequate
for signal reconstruction, and other types
of convex optimization must be invoked.
ACKNOWLEDGMENTS
This work was supported by grants from
NSF, DARPA, ONR, AFOSR, and the Texas
tangent/cosine configuration. The tan(θ)
computations were implemented by
CORDIC rotations avoiding the need for
multiply operations. To minimize output
phase angle error, we applied a post-
CORDIC clean-up angle rotation. Finally,
we stabilized the DDS output amplitude
by an AGC loop. The phase-noise per-
formance of the DDS is quite remarkable
and we invite you, the reader, to take a
careful look at its structure. A MATLAB-
code implementation of the DDS is avail-
able at http://apollo.ee.columbia.edu/
spm/?i=external/tipsandtricks.
ACKNOWLEDGMENT
Thanks to Rick Lyons for patience and
constructive criticism above and beyond
the call of duty.
AUTHOR
Fred Harris (fred.harris@sdsu.edu)
teaches DSP and modem design at San
Diego State University. He holds 12
patents on digital receivers and DSP
technology. He has written over 140
journal and conference papers and is the
author of the book Multirate Signal
Processing for Communication Systems
(Prentice Hall Publishing).
REFERENCES
[1] C. Dick, F. Harris, and M. Rice, “Synchronization
in software defined radios—Carrier and timing
recovery using FPGAs,” in Proc. IEEE Symp. Field-
Programmable Custom Computing Machines, Napa
Valley, CA, pp. 195–204, Apr. 2000.
[2] J. Valls, T. Sansaloni, A. Perez-Pascual, V. Torres,
and V. Almenar, “The use of CORDIC in software
defined radios: A tutorial,” IEEE Commun. Mag.,
vol. 44, no. 9, pp. 46–50, Sept. 2006.
[3] F. Harris, C. Dick, and R. Jekel, “An ultra low
phase noise DDS,” presented at Software Defined
Radio Forum Tech. Conf. (SDR-2006), Orlando FL,
Nov. 2006.
[4] R. Lyons, Understanding Digital Signal
Processing, 2nd ed. Upper Saddle River, NJ: Prentice
Hall, pp. 576–578, 2004.
[5] C. Turner, “Recursive discrete-time sinusoidal
oscillators,” IEEE Signal Processing Mag., vol. 20,
no. 3, pp. 103–111, May 2003.
[6] C. Dick and F. Harris, “FPGA signal processing using
sigma-delta modulation,” IEEE Signal Processing Mag.,
[SP]
vol. 17, no. 1, pp. 20–35, Jan. 2000.
Instruments (TI) Leadership University
Program. Special thanks are due to TI for
the DMD array used in the single-pixel
camera. Thanks also to the Rice DSP
group and Ron DeVore for many enlight-
ening discussions and Justin Romberg
for help with the reconstruction in
Figure 3.
AUTHOR
Richard G. Baraniuk (richb@rice.edu)
is the Victor E. Cameron Professor of
Electrical and Computer Engineering
at Rice University. His research inter-
ests include multiscale analysis,
inverse problems, distributed signal
processing, and sensor networks. He is
a Fellow of the IEEE.
REFERENCES
Additional compressive sensing resources are avail-
able at dsp.rice.edu/cs.
[1] E. Candès, J. Romberg, and T. Tao, “Robust
uncertainty principles: Exact signal reconstruction
from highly incomplete frequency information,”
IEEE Trans. Inform. Theory, vol. 52, no. 2,
pp. 489–509, Feb. 2006.
[2] D. Donoho, “Compressed sensing,” IEEE Trans.
Inform. Theory, vol. 52, no. 4, pp. 1289–1306, Apr.
2006.
[3] S. Mallat, A Wavelet Tour of Signal Processing.
New York: Academic, 1999.
[4] R.G. Baraniuk, M. Davenport, R. DeVore, and
M.B. Wakin, “A simple proof of the restricted isome-
try principle for random matrices (aka the
Johnson-Lindenstrauss lemma meets compressed
sensing),” Constructive Approximation, 2007
[Online].
http://dsp.rice.edu/
cs/jlcs-v03.pdf
Available:
[5] D. Baron, M.B. Wakin, M. Duarte, S. Sarvotham,
and R.G. Baraniuk, “Distributed compressed
sensing,” 2005 [Online]. Available: http://dsp.
rice.edu/cs/DCS112005.pdf
[6] J. Tropp and A.C. Gilbert, “Signal recovery
from partial information via orthogonal matching
pursuit,” Apr. 2005 [Online]. Available: http://www-
personal.umich.edu/~jtropp/papers/TG06-Signal-
Recovery.pdf
[7] J. Haupt and R. Nowak, “Signal reconstruction
from noisy random projections,” IEEE Trans. Inform.
Theory, vol. 52, no. 9, pp. 4036–4048, Sept. 2006.
[8] S. Kirolos, J. Laska, M. Wakin, M. Duarte, D.
Baron, T. Ragheb, Y. Massoud, and R.G. Baraniuk,
“Analog-to-information conversion via random
demodulation,” in Proc. IEEE Dallas Circuits
Systems Workshop, Oct. 2006, pp. 71-74.
[9] M. Vetterli, P. Marziliano, and T. Blu, “Sampling
signals with finite rate of innovation,” IEEE Trans.
Signal Processing, vol. 50, no. 6, pp. 1417–1428,
June 2002.
[10] D. Takhar, V. Bansal, M. Wakin, M. Duarte, D.
Baron, J. Laska, K.F. Kelly, and R.G. Baraniuk, “A
compressed sensing camera: New theory and an
implementation using digital micromirrors,” in
Proc. Comput. Imaging IV SPIE Electronic Imaging,
[SP]
San Jose, Jan. 2006.
IEEE SIGNAL PROCESSING MAGAZINE [124]
JULY 2007