Sparsity Order Estimation and its Application in ... - IEEE Xplore
2116
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
Sparsity Order Estimation and its Application in Compressive Spectrum Sensing for Cognitive Radios Yue Wang, Member, IEEE, Zhi Tian, Senior Member, IEEE, and Chunyan Feng
Abstract—Compressive sampling techniques can effectively reduce the acquisition costs of high-dimensional signals by utilizing the fact that typical signals of interest are often sparse in a certain domain. For compressive samplers, the number of samples Mr needed to reconstruct a sparse signal is determined by the actual sparsity order Snz of the signal, which can be much smaller than the signal dimension N . However, Snz is often unknown or dynamically varying in practice, and the practical sampling rate has to be chosen conservatively according to an upper bound Smax of the actual sparsity order in lieu of Snz , which can be unnecessarily high. To circumvent such wastage of the sampling resources, this paper introduces the concept of sparsity order estimation, which aims to accurately acquire Snz prior to sparse signal recovery, by using a very small number of samples Me less than Mr . A statistical learning methodology is used to quantify the gap between Mr and Me in a closed form via data fitting, which offers useful design guideline for compressive samplers. It is shown that Me ≥ 1.2Snz log(N/Snz + 2) + 3 for a broad range of sampling matrices. Capitalizing on this gap, this paper also develops a two-step compressive spectrum sensing algorithm for wideband cognitive radios as an illustrative application. The first step quickly estimates the actual sparsity order of the wide spectrum of interest using a small number of samples, and the second step adjusts the total number of collected samples according to the estimated signal sparsity order. By doing so, the overall sampling cost can be minimized adaptively, without degrading the sensing performance. Index Terms—Sparsity order estimation, curve fitting, compressive sampling, cognitive radio, wideband spectrum sensing.
Manuscript received March 22, 2011; revised August 14, 2011; accepted September 23, 2011. The associate editor coordinating the review of this paper and approving it for publication was M. Morelli. This work was supported in part by the US National Science Foundation (NSF) under Grant ECS-0925881, the National Science Foundation of China (NSFC) under Grant 60902047, and the National High Technology Research and Development Program of China (“863” Program) under Grant 2009AA011802. Part of this work was presented at the IEEE GLOBECOM Conference, Miami, FL, USA, Dec. 2010. Y. Wang was with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing 100876, China. He is now with the Research Department of Hisilicon, Huawei Technologies Co., Ltd., Beijing 100095, China (e-mail: [email protected]). This work was performed when Y. Wang was a visiting Ph.D. student at Michigan Technological University. Z. Tian is with the Department of Electrical and Computer Engineering, Michigan Technological University, Houghton, MI 49931, USA (e-mail: [email protected]). C. Feng is with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing 100876, China (e-mail: [email protected]). Digital Object Identifier 10.1109/TWC.2012.050112.110505
I. I NTRODUCTION OMPRESSIVE sampling (CS) is an emerging framework for efficient acquisition and processing of sparse signals. It deals with the reconstruction of sparse signals from a small number of compressive samples collected from linear projections of the signals [1], [2]. Compared with the traditional Nyquist sampling theory, the number of samples required in CS techniques can be made much smaller than the signal dimension by exploiting the signal sparsity property exhibited in a certain domain. An Snz -sparse signal s ∈ RN can be expressed as a linear projection on some proper representation basis Ψ ∈ RN ×N in the form of s = Ψθ, where the vector of expansion coefficients θ ∈ RN only has Snz nonzero entries with Snz N . Essentially, the task of CS is to recover θ from a reduced set of M (< N ) samples x ∈ RM that are acquired by undersampling s with a linear compression matrix Φ ∈ RM×N , as follows:
C
x = Φs = ΦΨθ
(1)
where the number of rows in Φ decides the number of compressive samples collected. Theoretical work indicates that θ can be recovered exactly from x when the multiplication of Φ and Ψ satisfies the so-called restricted isometry property (RIP) [1]–[3]. The RIP offers a lower bound on the minimum number of samples needed for sparse signal reconstruction. For example, for Gaussian and Bernoulli random matrices, a log-like expression M ≥ αSnz log (N/Snz ) holds for some constant α, which means that the number of samples needed for recovering θ and then s can be much smaller than the number of unknowns N when Snz N [3]. As a result, the minimum sampling rate R = RN q M/N for signal reconstruction is decided by Snz , given N and the Nyquist rate RN q . The smaller the sparsity order Snz is, the smaller M and R need to be in order to recover the signal precisely. However, in practice, the actual sparsity order Snz is usually unknown a priori or even dynamically changing; instead, the maximum sparsity order Smax can be made available as a statistical upper bound of Snz . Thus, the practical sampling rate has to be set conservatively based on Smax in lieu of Snz , which may result in an unnecessarily high sampling rate as Smax deviates from Snz . A remedy to this wastage of the sampling resources is by adopting sequential sampling techniques [4], [5]. In [4], one sample is taken at a time, and the sparse signal is recovered at each time based on batch processing of all the past samples. The
c 2012 IEEE 1536-1276/12$31.00
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
sampling process stops once the signal estimate converges. This method entails high computational costs, because batchform sparse signal recovery is performed at every time step. In [5], a recursive sparse signal recovery algorithm is developed to considerably reduce the computational cost compared with batch processing; yet, it aims to reduce the computation per step and track slowly-varying signals, rather than minimizing the sampling rate and sensing costs. To alleviate the wasteful sampling costs at low reconstruction load, this paper aims to learn the signal sparsity order on the fly by introducing a novel concept of sparsity order estimation. The sparsity order estimation problem has not been studied in the CS literature. Indeed, most work on the minimum number of samples needed in CS can be divided into two categories: one for sparse signal estimation and the other for sparse support estimation. The former seeks to exactly reconstruct the sparse signal including both the positions and the amplitudes of all nonzero elements [1]– [3], while the latter aims to identify the positions of these nonzero elements [6], [7]. Evidently, both estimation tasks demand more sampling resources than the sparsity order estimation problem, because the latter only needs to identify the number of nonzero elements in the sparse signal, but not the position or the amplitude of each nonzero element. The tenet of CS suggests the possibility that signals be sampled at a rate proportional to the amount of information content to be extracted, which alludes to a lower sampling rate for sparsity order estimation. In this paper, a general methodology is provided to quantify the minimum number of samples required to estimate the sparsity order for different types of sampling matrices. To get around the difficulty in theoretical order statistics analysis for the sparsity order estimation problem, a simulation-based approach coupled with the curve fitting technique is employed to depict the fundamental limits of sampling rates in fitted closed-form expressions. The study confirms that the sparsity order estimation problem requires a smaller number of samples than that for sparse signal or support estimation problems. The nontrivial gap in terms of the sampling resources needed for sparsity order estimation versus sparse signal reconstruction can be useful in signal processing applications. This paper offers an illustrative example by developing a twostep compressive spectrum sensing (TS-CSS) algorithm for wideband cognitive radios (CRs). The TS-CSS first estimates the sparsity order Snz of the unknown wideband spectrum using a small number of samples collected at the first step, and then uses the estimated Sˆnz to decide the number of additional samples to be collected at the second step. Afterward, all samples from both steps are used to accurately reconstruct the wideband spectrum and eventually make the final decision on the spectrum occupancy. Overall, the TS-CSS adaptively adjusts its sampling rate according to Sˆnz estimated on the fly. As a result, the average sampling rate of the TS-CSS can be considerably lower than the traditional compressive sampling rate determined conservatively by Smax . Such a saving in acquisition costs is confirmed via simulations, given the desired sensing accuracy. Noticeably, reducing the sampling costs in CS applications can also lower the computational burden of the sparse signal recovery operations, which reduces the overall
2117
(hardware) implementation costs in practice. The rest of this paper is organized as follows. Section II describes the general principle to quantify the gap between the minimum number of samples needed for sparsity order estimation and that for sparse signal reconstruction. Section III illustrates the usefulness of such a gap by developing a twostep compressive spectrum sensing algorithm for wideband CRs. Performance evaluation results are presented in Section IV, followed by conclusions in Section V. II. F UNDAMENTAL L IMITS IN S PARSITY O RDER E STIMATION This section starts by introducing a novel concept of sparsity order estimation, which aims to identify the number of nonzero elements of a sparse signal without having to know their positions or amplitudes. To validate this concept, an empirical study is set up to evaluate the minimum number of samples needed for accurate sparsity order estimation. Then, a general methodology based on the curve fitting technique is developed to delineate the closed-form expressions for several widely-used random sampling matrices. This study provides quantitative results on the nontrivial gap between the minimum number of samples required for sparsity order estimation and that for sparse signal reconstruction. A. Concept of Sparsity Order Estimation The sparsity order of an N × 1 sparse signal vector θ is defined by its 0 -norm in the form Snz = θ0 , which refers to the number of nonzero elements of θ. The sparsity order estimation problem concerns estimating Snz from the compressive samples x in (1). The number of samples needed for sparsity order estimation depends on not only N and Snz , but also the matrices Φ and Ψ, which are discussed next. B. Options of Matrix Patterns The sample quality and estimation performance are determined by the matrices used in the sampling process. Here several widely-used matrix patterns are reviewed for later experimental use. 1) Sparse Representation Basis Ψ: In the context of CS, a signal may have a sparse representation when it is expressed on a proper basis characteristic of its specific application, e.g., Fourier basis for signals with underutilized frequency spectrum [8], wavelet basis and discrete curvelet transform basis for natural image signals [1], [9], etc. This paper will study frequency-sparse signals as an illustrative example in underutilized wideband CR networks; correspondingly the focus is on the Fourier basis in the form of an N × N discrete Fourier transform (DFT) matrix. For testing the universality of our methodology, the identity matrix will also be studied as a special case of the representation basis, which is employed when the signal sparsity appears in the same domain of the sampling process. 2) Random Sampling Matrix Φ: It has been shown that Gaussian and Bernoulli random sampling matrices have a high probability of satisfying the RIP property [3]. It means that samples acquired from such matrices have a better chance to
2118
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
perfectly recover sparse signals than other types of sampling matrices, given the same number of samples M . The entries {φi,j }i,j of these two random sampling matrices are independent and identically distributed according to Gaussian and Bernoulli probability distributions as follows (w.p. stands for with probability): Gaussian: Bernoulli:
φi,j ∼ N (0, 1/M ) ; √ +1/ M w.p. 1/2, √ φi,j = −1/ M w.p. 1/2.
(2) (3)
There are several options for hardware implementations of random sampling matrices [10]. In practical settings, structured random matrices are often employed for improved implementation affordability [11]–[13]. An exemplary sampler that will be tested in this paper is the so-called analogto-information converter (AIC) [11], which compresses and samples analog signals at uniform sub-Nyquist sampling rates. It has a serial concatenated structure with three main components: a wideband modulator made of a pseudo-random chipping sequence, a low-pass filter, and a uniform low-rate conventional analog-to-digital converter (ADC). The resulting sampling matrix Φ turns out to be a banded random matrix with a great deal of structure. Because of the special structure, the AIC alleviates the sampling burden on the ADC component, at the expense of slightly degraded recovery performance compared with those fully random and dense Gaussian or Bernoulli samplers. C. Experimental Design We now present a general methodology to determine the minimum number of samples Me that is required to reliably identify the sparsity order Snz of a sparse signal θ, while that needed to reconstruct the signal itself Mr is also provided as a benchmark. The developed approach is simulation-based, in which a large set of Monte Carlo simulation trials is performed following the steps below, and the trial results are used to fit into closed-form expressions quantifying Me and Mr as functions of the signal size N and the actual sparsity order Snz , for a given measurement matrix ΦΨ. Simulation-based quantification is a valid statistical learning tool that has been theoretically justified and fruitfully adopted in e.g., Markov Chain Monte Carlo (MCMC) methods for Bayesian networks [14]. In CS research, [11] adopts a simulation-based approach to quantify the performance bound of practical AIC samplers. Even when the RIP condition is derived analytically to come up the well-known expression Mr ≥ αSnz log(N/Snz ) for Gaussian and Bernoulli sampling matrices, the theoretical analysis is generally lacking concerning the precise value of α [1]–[3]. If a specific value for α is required, then one typically must determine it experimentally. Besides, theoretical results tend to work in an asymptotic regime where the probabilistic characteristics of the random matrices are valid, and hence are useful when the problem size is very large but may appear to be too conservative for practical problems of reasonably small to medium sizes. For these reasons, this paper develops experimental design along with the curve fitting technique to bypass the difficulty in theoretically analyzing the sparsity order estimation problem. The developed methodology results
in usable expressions that fit the simulation data into a closed form, and is general enough to be readily applied offline for a wide range of sampling matrices under practical operating regimes, without incurring high costs in computation and analysis. The following steps are performed in each Monte Carlo test. 1) Generate a Sparse Signal: An Snz -sparse signal vector θ of size N is generated randomly according to a generic stochastic model. That is, given a preset integer Snz ( N ), the sparse support of θ is a uniformly distributed random set of the Snz nonzero positions as follows: Supp (θ) = {{k1 , . . . , kSnz } |θki = 0, ki ∈ {1, . . . , N }} where
Pr(ki ∈ Supp(θ)) = Snz /N,
i = 1, . . . , N.
(4)
Meanwhile, the amplitude of each element of θ is set according to one of the following commonly-used distributions:1 = 1, ki ∈ Supp (θ) , Binary: θki (5) = 0, ki ∈ / Supp (θ) ; ∼ N (1, 0.01) , ki ∈ Supp (θ) , Gaussian: θki (6) = 0, ki ∈ / Supp (θ) ; ∼ U (0.8, 1.2) , ki ∈ Supp (θ) , or Uniform: θki (7) = 0, ki ∈ / Supp (θ) . 2) Collect Compressive Samples: With the randomly generated θ and the chosen Ψ and Φ, M measurements are generated from (1) to form a length-M sample vector x. ˆ of the 3) Reconstruct the Sparse Signal: An estimate θ original θ is computed using the basis pursuit (BP) algorithm for the noise-free case, which offers a tractable solution by minimizing the 1 -norm of θ subject to a linear projection constraint [1], as follows: ˆ = arg min θ s.t. x = ΦΨθ. θ 1
(8)
For the noisy case, the least absolute shrinkage and selection operator (LASSO) algorithm can be used to yield [15] ˆ = arg min θ s.t. x − ΦΨθ ≤ θ 1 2
(9)
where bounds the amount of noise energy in the sample vector. 4) Estimate the Sparsity Order: Having obtained the signal ˆ the sparsity order Snz can be estimated as estimate θ, Sˆnz =
N ˆ θ [i] ≥ λ
(10)
i=1
where λ > 0 is a preset decision threshold. Empirically, the threshold is set as λ = μ/2 for the noise-free case or λ = (μ + σn ) /2 for the noisy case, where μ is the average absolute value of all the Snz nonzero elements in θ, σn is the standard deviation of the noise component, and λ lies in the middle of μ and σn . As Fig. 1 indicates, the amplitudes of the ˆ are well separated from those of the nonzero elements in θ 1 This paper focuses on these simple and widely-used probabilistic distributions, while the amplitude setting can be extended to other distributions as well. The choice of the amplitude distribution has little impact on the simulation-based approach, because θ can be recovered accurately for any amplitudes, once the nonzero support is correctly identified.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
70
0.5
60
0
0
5
10
15 (a)
20
25
30
1 0.5 0
0
5
10
15 (b)
20
25
30
Minimum number of samples
1
50
40 min Me (Gaussian) min Me (Bernoulli) min Me (AIC) min Mr (Gaussian) min Mr (Bernoulli) min Mr (AIC) 1.2Slog(N/S+2)+3.2 1.32Slog(N/S+4)+2.6 1.71Slog(N/S+1)+0.9
30
20
1
=0.5
0.5 0
10
0 0
5
10
15 (c)
20
25
30
1
0
0
5
10
15 (d)
20
0
5
10 Sparsity order
15
20
Fig. 2. Minimum Me and Mr vs. Snz for sparsity order estimation and sparse signal reconstruction with different sampling matrices, using DFT representation basis and N = 128 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
=0.5
0.5
2119
25
30
D. Quantification Methodology Fig. 1. Illustration of a noise-free case when N = 30, Snz = 5, and λ = 0.5. (a) The sparse input signal generated from (4) and (5), (b) a successful reconstruction of the input signal, (c) a successful estimation of both the nonzero support and the sparsity order, although the amplitudes are not recovered accurately, and (d) a successful estimation of sparsity order, even though errors occur in both the amplitudes and the nonzero support.
zero elements, which makes the decision in (10) robust to the threshold value λ. This is because the 1 -norm penalty term on θ enforces a sparse solution by shrinking small estimates to zero. With (10), a successful identification of the sparsity order is declared when the estimated sparsity order is equal to the true value, i.e., Sˆnz = Snz . To identify the minimum number of samples Me needed for reliable sparsity order estimation, a range of values for M is tested, and Me is given by the smallest M that results in a success rate of 99% in correctly identifying Snz among 500 Monte Carlo trials. The simulation-based approach can also be used to empirically observe the minimum number of samples Mr needed to accurately reconstruct the sparse input signal. In that case, ˆ = θ to machine a successful recovery is declared when θ precision, and the smallest value Mr is reported at which the probability of empirical successes is greater than 99% in 500 Monte Carlo trials. Fig. 1 illustrates several simulation outcomes of the sparse signal recovery problem in Step 3) and the sparsity order estimation problem in Step 4). Apparently, given the same sampling resources, the probability of successful sparsity order estimation can be higher, because Sˆnz = Snz is possible even when the amplitudes or positions of those nonzero elements are recovered inaccurately.
According to the experimental design described above, we now analyze the empirical results on the minimum numbers of samples Me and Mr for the noise-free case. Specifically, we seek to delineate Me and Mr as functions of the sparsity order Snz and the signal size N , given the sampling matrix pattern. As suggested by both analysis [2], [3] and experiments [11], the desired function for Mr typically takes on a log-like form, which also conforms well with the empirical data we have observed for both Mr and Me . Hence, our quantification methodology is to adopt the general expression ASnz log(N/Snz + B) + C, and use experimental results and statistical curve fitting techniques to find the values of A, B and C that minimize the data fitting error. First, we fix N and find the expressions of Me and Mr as functions of Snz . Consider for example the case of N = 128 and using the DFT representation basis, while Snz varies from 1 to 20. For different sampling matrices including Gaussian, Bernoulli and AIC, the empirical data of Me and Mr as well as the fitted expressions using statistical regression are depicted in Fig. 2. The fitted expression of the minimum Me for sparsity order estimation is given by Me,min ≈ 1.2Snz log (N/Snz + 2) + 3.2. Meanwhile, the fitted expressions of the minimum Mr for sparse signal recovery are in the form of Mr,min ≈ 1.32Snz log (N/Snz + 4) + 2.6 for Gaussian and Bernoulli sampling matrices, and Mr,min ≈ 1.71Snz log (N/Snz + 1) + 0.9 for the AIC sampling matrix [11]. Apparently, in order to achieve the same reconstruction accuracy, a structured AIC random sampler requires a bit more samples than samplers using Gaussian and Bernoulli random matrices, in return for simplified hardware implementation. Alternatively, we fix Snz and quantify the two numbers Me and Mr as functions of N . Fig. 3 shows the empirical results and the fitted curves for Snz = 10, where N increases from 64 to 1024. The fitted expressions from both Fig. 2 and Fig. 3 are near identical.
2120
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
100
80 70 60
90 80 Minimum number of samples
Minimum number of samples
90
100 min Me (Gaussian) min Me (Bernoulli) min Me (AIC) min Mr (Gaussian) min Mr (Bernoulli) min Mr (AIC) 1.2Slog(N/S+2)+2.9 1.32Slog(N/S+4)+3.5 1.69Slog(N/S+1)+0.3
50 40
60 50 40 30
30 20
70
min Me (Gaussian) min Mr (Gaussian) 1.5Slog(N/S+2)+6.6 2.5Slog(N/S−1)+10.9
20
2
10
3
10
10 Signal dimension
Fig. 3. Minimum Me and Mr vs. N for sparsity order estimation and sparse signal reconstruction with different sampling matrices, using DFT representation basis and Snz = 10 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
E. Fundamental Limits in Sample Rates The following Results (1-3) summarize the curve fitting results, from the two sets of experiments depicted in Fig. 2 and Fig. 3. Result 1: Successful estimation of the sparsity order can be achieved with 99% probability when the number of samples Me satisfies the lower bound: Me ≥ 1.2Snz log (N/Snz + 2) + 3.
(11)
Result 2: To obtain successful sparse signal reconstruction with 99% probability, the lower bounds of Mr for different samplers are Mr ≥ 1.32Snz log (N/Snz + 4) + 3
(12)
for Gaussian/Bernoulli samplers, and Mr ≥ 1.7Snz log (N/Snz + 1)
(13)
for the AIC sampler [11]. Comparison of Result 1 and Result 2 leads to Result 3. Result 3: For any given Snz and N , it always holds that Mr ≥ Me . The expressions (11)-(13) are concluded when a size-N DFT matrix is used as the sparse representation basis Ψ. This basis matrix is relevant to the compressive spectrum sensing problem to be investigated in Section III. The general principle of statistical learning can be used to quantify the fundamental sampling costs for any matrix patterns. For example, when Ψ = I is an identity matrix, the empirical results and fitted curves are illustrated in Fig. 4. The required sampling costs appear to be higher than those in (11)-(13), even though the general trends are quite similar. In general, the proposed simulation-based approach can be conducted offline for any choices of matrix patterns, and offers simple expressions to guide practical sampling systems.
0
5
10 Sparsity order
15
20
Fig. 4. Minimum Me and Mr vs. Snz for sparsity order estimation and sparse signal reconstruction with Gaussian random matrix as sampler and identity matrix as the sparse representation basis, for N = 128 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
III. A PPLICATION IN C OMPRESSIVE S PECTRUM S ENSING FOR W IDEBAND CR S As a new paradigm for efficient acquisition and processing of sparse signals, CS has found attractive applications in a broad range of signal processing problems such as image processing [16], [17], medical imaging [18], wireless sensor networks [19], cognitive spectrum sensing [20]–[24], to name a few. As an illustrative example, this section presents the application of sparsity order estimation in compressive spectrum sensing for wideband CRs. In CR networks where users dynamically share the spectrum, CRs need to rapidly and accurately identify spectrum opportunities over a very wide band. Subject to analog-todigital hardware limitations in the over-GHz regime, CS techniques are well motivated for wideband spectrum detection exploiting sparsity, because there are plenty of temporarily unused spectrum opportunities in wideband CR networks [25]. Several CS-based approaches have been developed to detect the frequency occupancy of primary users using sub-Nyquistrate samples [20], [21]. To determine the suitable sampling rate, they all implicitly assume that (the upper bound of) the sparsity order of the underutilized spectrum is known exactly. However, in practical CR applications, the actual sparsity order Snz corresponds to the instantaneous spectrum occupancy of wireless users which is dynamically varying. Hence, Snz is often unknown; instead, what is typically available is its upper bound Smax , which can be measured from the maximum spectrum utilization observed statistically over a time period. In practice, conservative determination of the sampling rate based on Smax can cause unnecessarily high sensing costs. To alleviate such wasteful acquisition costs and keep low reconstruction load, we develop a two-step compressive spectrum sensing (TS-CSS) algorithm for wideband CRs [23], which capitalizes on our key observation of the gap between the number of samples needed for sparsity order estimation and that for sparse signal reconstruction.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
A. System Model and Problem Statement
70
(14)
where xt of length M consists of the collected compressive samples, rt of length N is the discrete-time representation of r(t) at Nyquist rate with elements rt [n] = r (t) | t=n/RN q , n = 0, . . . , N − 1, F−1 is the N × N inverse DFT matrix, and Φ is the M × N (M < N ) AIC sampling matrix. Evidently, the actual sampling rate R = RN q M/N in collecting the sample vector xt is lower than the Nyquist rate RN q , while the recovery of the unknown rt or rf in (14) will enable spectrum sensing at full frequency resolution. B. Two-Step Compressive Spectrum Sensing (TS-CSS) In the presence of an unknown Snz , a TS-CSS algorithm is now developed for wideband spectrum sensing at low sampling costs. The basic idea is to utilize the nontrivial gap between the minimum numbers of samples for sparsity order estimation and for sparse signal reconstruction (c.f. Fig. 5). To do so, each sensing block is divided into two time slots of length T1 and T2 , used for sparsity order estimation and signal detection, respectively. A typical choice is to set T1 = T2 = T /2, while the effect of non-equal time division will be discussed later. The overall procedure of the TS-CSS algorithm is shown in Fig. 6, and details of each step are elaborated next.
Mc
60
Minimum number of samples
Suppose that a wideband primary network spans over a total of B Hz, which consists of N non-overlapping subchannels of equal bandwidth. At a given time and location, the network is underutilized with only Snz ( N ) channels occupied by active primary users, whereas the rest N − Snz channels are idle and offer spectrum opportunities to secondary CR users. Clearly, the sparse nature of the signal of interest exhibits in the frequency domain. However, neither the locations of the occupied Snz channels, nor the sparsity order Snz , are known a priori. Instead, only an upper bound Smax of the sparsity order can be available in practice. Given the received signal r (t) of bandwidth B Hz and a coarsely known Smax , the wideband spectrum sensing task for a CR detector is to detect the occupancy of the N spectrum bands. Those idle frequency bands present potential transmission opportunities for CRs, if and only if they can be detected correctly. A CR detector adopting CS techniques starts with an analog-to-digital sampler that collects M discrete-time samples from r (t). Here we utilize the AIC as the random sampler and the DFT matrix as the sparse representation basis. During spectrum sensing, each time window for sensing shall be small enough compared with the channel coherence time and user dynamics, such that transmission opportunities can be utilized before the spectrum occupancy profile has changed. It is hence justified to assume that the unknown frequency spectrum of interest is block stationary for each sensing window of time length T = N/RN q , given N and RN q . Spectrum sensing is performed block by block, aiming to recover the N × 1 unknown frequency response vector rf of invariant strengths and Snz in each block. Hence, the discrete-time representation of the CS sampling operation on the received r (t) at a CR detector can be expressed in the following general form: xt = Φrt = ΦF−1 rf
2121
M c M
Saving
50
^
M
`M
M M1
2
M1
For sparse signal reconstruction
40
30
20
For sparsity order estimation
10
0
S nz 0
2
4
6
8 10 12 Sparsity order
S max 14
16
18
20
Fig. 5. Basic idea of the TS-CSS utilizing the gap between the numbers of samples required for sparsity order estimation at the first step and for sparse signal reconstruction at the final step.
r (t )
Estimate sparsity order
M1
Reconstruct sparse signal
S nz
rˆ f
Make sensing decision
d
M2
J
S max
Fig. 6.
Predict sampling rate
Flow chart of the TS-CSS algorithm.
1) Estimate the Sparsity Order at the First Step: With the available Smax and using the solid curve for Me in Fig. 5, the first step is to decide M1 which is the number of samples needed for estimating the sparsity order within the time slot T1 . From Result 1, and substituting the unknown Snz in (11) by its upper bound Smax , M1 is decided as follows: M1 = 1.2Smax log (N/Smax + 1/2) .
(15)
Here, M1 is chosen conservatively according to Smax , such that it suffices to accurately estimate the sparsity order of any signal with Snz ≤ Smax . Since the M1 samples are acquired during the first half period of r(t), the sampling process at the first step can be described as xt1 = [Φ1 0] rt = [Φ1 0] F−1 rf
(16)
where Φ1 is the M1 × N/2 sampling matrix. Given xt1 , the TS-CSS aims to estimate Sˆnz of the sparse spectrum rf . This can be done by first recovering ˆrf via (8), where θ is replaced by rf and (16) is employed as the linear
2122
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
constraint. Then, following (10), Sˆnz is estimated as Sˆnz =
N
(|ˆrf [i]| ≥ λ) .
Algorithm: Two-Step Compressive Spectrum Sensing (TS-CSS) (17)
i=1
2) Sample Adaptively at the Second Step: At the second step of the TS-CSS, additional M2 samples are collected into a data vector xt2 over the second time period of r(t), using an M2 × N/2 sampling matrix Φ2 as follows: xt2 = [0 Φ2 ] rt = [0 Φ2 ] F−1 rf .
(18)
The primary issue at the second step is to decide the number M2 such that rf can be accurately estimated from a total of M = M1 + M2 compressive samples. This can be done by using the dashed curve in Fig. 5 along with the estimated Sˆnz . Specifically, replacing Snz in (13) by Sˆnz , M can be decided as M = 1.7Sˆnz log(N/Sˆnz + 1).
(19)
If the upper bound Smax is loose such that M1 ≥ M , then the M1 samples collected at the first step are already adequate for reliably recovering the sparse spectrum, and hence there is no need to sample further at the second step. Thus, M2 is decided as M2 = max(M − M1 , 0)
(20)
where M1 and M are calculated from (15) and (19), respectively. 3) Reconstruct the Signal Spectrum from All the Collected Samples: After the two-step sampling process, the unknown wideband spectrum rf can be reconstructed from all available samples using the 1 -norm minimization formula in (8). Concatenating the sample vectors xt1 and xt2 in (16) and (18) from both steps, (8) becomes ˆrf = arg min rf 1 s.t.
xt1 xt2
=
Φ1 0
0 Φ2
F−1 rf .
(21) 4) Make Spectrum Occupancy Decisions: Having obtained ˆrf from (21), the CR detector can finally make decisions on the spectrum occupancy by comparing the signal energy within each frequency band with a decision threshold γ, as follows: ˆ = |ˆrf |2 ≥ γ . d
(22)
ˆ of length N indicates the sensing The binary decision vector d outcomes on the monitored wideband spectrum. If an entry of ˆ is 0, then the corresponding frequency band is decided to d be idle, presenting a potential spectrum opportunity for CRs. If the entry is 1, then the corresponding frequency band is decided to be occupied by primary users. The above steps are summarized in Algorithm TS-CSS.
• • • • • •
set T1 = T2 = T /2; collect M1 samples from r(t) over the first time period of T1 seconds, where M1 is decided by (15); estimate Sˆnz from (17); collect M2 samples from r(t) over the second time period of T2 seconds, where M2 is decided by (20); estimate the sparse signal spectrum ˆrf from (21); ˆ from (22). make the spectrum occupancy decision d
C. Discussions Several remarks on this TS-CSS algorithm are in order. 1) Saving in Sampling Costs: It is of interest to compare the sampling cost of the TS-CSS with that of a traditional onestep CS approach, which would directly use Smax to decide the total number of samples M used for recovering the sparse signal within one step. Essentially, M is given by Result 2 after replacing the unknown Snz by its upper bound Smax , which is indicated in Fig. 5. Apparently, M exceeds the total number M used in the TS-CSS, resulting in higher samplingrate requirements for the one-step approach. The difference (M − M ) reflects the saving in sampling costs. Notice that if the prior knowledge is already very good such that Smax = Snz , then the saving decreases to 0. On the other hand, given a certain Snz , the largest amount of saving happens when M1 ≥ M , in which case there is a big gap between the prior knowledge and the actual sparsity order such that Smax Snz . At least, the saving is positive as long as Smax > Snz , which takes place in most cases. The quantitative saving in sampling costs will be presented in next section via computer simulations. 2) Time Ratio T1 /T : The ratio T1 /T may impact the sensing efficiency. Intuitively, an assignment that consumes fewer samples and shorter length in the first time slot results in more efficient sensing, since Snz can be quickly estimated to adaptively decide the number of additional samples required in the second time slot. Fig. 7 presents the minimum numbers of samples required in the first time slot for different values of the time ratio T1 /T . The observed data and corresponding fitted curves therein are obtained by the same methodology developed in Section II. Fig. 7 indicates that the lower the ratio is chosen, the fewer samples are needed to successfully estimate the sparsity order. This result can be explained by the higher peak sampling rate resulted from a shorter length of the first time slot, as shown in Fig. 8. However, the drawback of a small ratio is the limited range of sparsity order which can be estimated; for example, the range of predictable Snz is only from 1 to 12 for the ratio T1 /T = 1/4 in Fig. 8, which also explains why the curve corresponding to the ratio T1 /T = 1/4 has a shorter range than its counterparts in Fig. 7. Therefore, there exists a tradeoff between the sampling efficiency and the applicable range of the TS-CSS algorithm. For simplicity and without loss of generality, T1 /T = 1/2 is adopted. 3) Implementation of Sampling Rate Switching: In TSCSS, there are two different sampling rates for the two successive steps, each of which is a uniform sub-Nyquist rate. The rate switching can be conveniently incorporated in practical
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
80
Minimum number of samples
70 60 50 40 T1/T=1 T1/T=3/4 T1/T=1/2 T1/T=1/4 1.2Slog(N/S+2)+3 1.2Slog(N/S+1.5)+2.9 1.2Slog(N/S+0.5)+0.1 1.2Slog(N/S−1.5)+0.2
30 20 10 0
0
5
10
15 20 Sparsity order
25
30
35
Fig. 7. Minimum number of samples vs. sparsity order for sparsity order estimation, for different values of the ratio T1 /T when N = 128 in noise-free case.
1 0.9
Peak sampling rate / RNq
0.8 0.7 0.6
T1/T=1 T1/T=3/4 T1/T=1/2 T1/T=1/4 (1.2Slog(N/S+2)+3)/N 4/3(1.2Slog(N/S+1.5)+2.9)/N 2(1.2Slog(N/S+0.5)+0.1)/N 4(1.2Slog(N/S−1.5)+0.2)/N
0.4 0.3 0.2 0.1 0
Consider a wideband CR system with N = 1024, Smax = 20, and Snz = 15. The sparse nonzero support is chosen as in (4), and the amplitudes of the nonzero elements are Gaussian distributed as in (6). The AIC sampler is used and the LASSO algorithm is employed for sparse signal recovery in the presence of ambient noise. The key design parameters are: M1 is computed from (15), M is computed from (19), and M2 = M − M1 . Sensing performance metrics of interest are the probability of miss-detection Pm and that of false-alarm Pf . These two detection criteria indicate the level of interference to primary users and the wastefulness of available spectrum opportunities for CRs, respectively. Given the true state vector d ∈ {0, 1}N of the target wideband spectrum, Pm and Pf can be expressed as follows: ˆ ˆ dT (d = d) (1 − d)T (d = d) , Pf = E Pm = E 1T d N − 1T d (23) where E denotes expectation, and 1 denotes the length-N allone vector. Practical CR applications, e.g., IEEE 802.22, concern the quality of service (QoS) of coexisting primary systems. To this end, it is required that Pm should be no higher than a predefined QoS-specific threshold Pm,th , given a tolerable Pf . Hence, the sensing performance can be evaluated by the following probability: Pr (Pm ≤ Pm,th ) , given a fixed Pf .
0.5
0
5
10
15
20 25 Sparsity order
30
35
40
45
Fig. 8. Peak sampling rate vs. sparsity order for sparsity order estimation, for different values of the ratio T1 /T when N = 128 in noise-free case.
hardware implementations, without causing much complication. For example, a serial random sampler, such as the AIC sampler [11], can flexibly adjust its sampling rate at the last ADC block to yield the two different effective sampling rates for the two steps, specifically, R1 = RN q M1 /(N/2) at the first step, and R2 = RN q M2 /(N/2) at the second step. Alternatively, in a parallel random sampler implementation for general sampling matrices [10], a punctured switch can be used to selectively admit different numbers of parallel branches in the two steps, specifically M1 and M2 branches, respectively. IV. P ERFORMANCE E VALUATIONS This section presents simulation results to verify the sensing performance and sampling costs of the proposed TS-CSS, using the traditional one-step CS approach as a benchmark.
2123
(24)
In (24), the two sensing performance metrics Pm and Pf are jointly considered to reflect the best-effort QoS satisfaction. In the simulations, both Pf and Pm,th are set to be 0.1, which means that the achieved probability of detection shall exceed 90%. Fig. 9 depicts the sensing performance versus signal-tonoise ratio (SNR) in terms of the best-effort QoS satisfaction. The SNR is defined to be the total signal power of the wideband signal over the entire spectrum, scaled by the power of the white noise. For a given SNR, the total number of samples used is averaged over 500 Monte Carlo trials and marked along the curves to show the sampling costs. In Fig. 9, it is clear that the best-effort QoS satisfaction of the TS-CSS is very close to that of the traditional one-step CS approach, when both of them use the same prior knowledge of Smax . As SNR increases, both of the probabilities of QoS satisfaction approach 1. Meanwhile, the total number of samples used in the TS-CSS reduces to M = 108 corresponding to an average sampling rate of R = (108/1024)RN q , whereas the traditional one-step CS approach samples at a higher rate of R = (135/1024)RN q throughout. The reduction in samples means lower sensing costs and less energy consumption, both of which are important merits for CR detectors in practice. Fig. 10 evaluates the saving rate of the acquisition costs offered by the TS-CSS over the one-step approach, for various Snz and Smax . The saving in sampling costs is measured by the following rate: η=
M − max (M1 , M ) × 100%. M
(25)
2124
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
V. C ONCLUSIONS
1
108 114.7 109.7 116.8 110.1 108 108 117.2
0.9 0.8
118.7 135
Pr(Pm Pm,th)
0.7 0.6
119.5
0.5 0.4 119.7 0.3 0.2 119.9 0.1
traditional one−step CS proposed TS−CSS
123.7121.9 126.1 120.7
0 −6
−4
−2
0
2 SNR
4
6
8
10
Fig. 9. Comparison of the proposed TS-CSS and the traditional one-step CS approach in terms of the probability of QoS satisfaction and the total number of samples used for different SNR cases when N = 1024, Smax = 20, and Snz = 15.
The cornerstone of this work is the fact that the number of samples required for estimating the signal sparsity order can be much smaller than that for reconstructing the unknown sparse signal itself. To reveal such a nontrivial gap, this paper first presents a general methodology to quantify the gap via simulation-based probabilistic evaluation and the curve fitting technique. The fundamental limits in sparsity order estimation are delineated by data-fitted closed-form expressions with loglike form, for some widely-used sampling matrices. Next, to illustrate the usefulness of this gap in practical applications, a two-step compressive spectrum sensing (TS-CSS) algorithm is developed for wideband CRs that only coarsely know the spectrum sparsity order. The TS-CSS algorithm adaptively switches the sampling rates for the two sequential steps of sparsity order estimation and sparse spectrum reconstruction, which effectively avoids collecting unnecessary samples caused by the discrepancy between actual and maximum sparsity orders. Compared with the traditional one-step CS-based spectrum sensing, the TS-CSS utilizes the smallest possible number of samples based on the estimated sparsity order on the fly, and can hence considerably reduce the sampling costs while achieving the desired sensing accuracy.
30 Smax=10 S Saving rate of sampling costs (%)
R EFERENCES
=20
max
25
Smax=30
20
15
10
5
0
0
5
10
15 Snz
20
25
30
Fig. 10. Saving rate η of sampling costs achieved by the proposed TS-CSS for different Snz and Smax when N = 1024.
As Fig. 10 indicates, the saving in sampling costs is most significant when Snz is far off from Smax . When Smax − Snz is large, the prior knowledge about sparsity order is coarse, which results in a large gap between the minimum number of samples for sparsity order estimation and that for sparse signal reconstruction; correspondingly, this is considerable room for saving the sampling costs, and the proposed TSCSS collects the benefit of this gap. The maximum saving in sampling costs stays unchanged when Snz Smax namely Snz ≤ Smax ×60% in Fig. 10. This is because the M1 samples collected at the first step in order to estimate Sˆnz are already sufficient to recover the sparse signal accurately. On the other hand, η approaches 0 as Snz gets close to Smax , which means diminished saving in sampling costs when the prior knowledge is close to the actual sparsity order.
[1] E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, Mar. 2008. [2] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, Apr. 2006. [3] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “A simple proof of the restricted isometry property for random matrices,” Const. Approx., vol. 28, no. 3, pp. 253–263, Dec. 2008. [4] D. M. Malioutov, S. Sanghavi, and A. S. Willsky, “Compressed sensing with sequential observations,” in Proc. 2008 IEEE ICASSP Conf., pp. 3357–3360. [5] D. Angelosante and G. B. Giannakis, “RLS-weighted lasso for adaptive estimation of sparse signals,” in Proc. 2009 IEEE ICASSP Conf., pp. 3245–3248. [6] M. J. Wainwright, “Information-theoretic limits on sparsity recovery in the high-dimensional and noisy setting,” IEEE Trans. Inf. Theory, vol. 55, no. 12, pp. 5728–5741, Dec. 2009. [7] W. Wang, M. J. Wainwright, and K. Ramchandran, “Informationtheoretic limits on sparse support recovery: dense versus sparse measurement,” IEEE Trans. Inf. Theory, vol. 56, no. 6, pp. 2967–2979, June 2010. [8] E. J. Candes, J. K. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006. [9] E. J. Candes and L. Demanet, “The curvelet representation of wave propagators is optimally sparse,” Comm. Pure Appl. Math., vol. 58, no. 11, pp. 1472–1528, Nov. 2005. [10] M. Mishali and Y. C. Eldar, “From theory to practice: sub-Nyquist sampling of sparse wideband analog signals,” IEEE J. Sel. Topics Signal Process., vol. 4, no. 2, pp. 375–391, Apr. 2010. [11] J. A. Tropp, J. N. Laska, M. F. Duarte, J. K. Romberg, and R. G. Baraniuk, “Beyond Nyquist: efficient sampling of sparse bandlimited signals,” IEEE Trans. Inf. Theory, vol. 56, no. 1, pp. 520–544, Jan. 2010. [12] J. Romberg, “Compressive sensing by random convolution,” SIAM J. Imag. Sci., vol. 2, no. 4, pp. 1098–1128, Dec. 2009. [13] X. Chen, Z. Yu, S. Hoyos, B. M. Sadler, and J. Silva-Martinez, “A subNyquist rate sampling receiver exploiting compressive sensing,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 58, no. 3, pp. 507–520, Mar. 2011. [14] C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd edition. Springer-Verlag, 2004. [15] R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. Roy. Statist. Soc. Ser. B, vol. 58, no. 1, pp. 267–288, 1996.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
[16] J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 14–20, Mar. 2008. [17] M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 83–91, Mar. 2008. [18] J. Trzasko, A. Manduca, and E. Borisch, “Highly undersampled magnetic resonance image reconstruction via homotopic 0 -minimization,” IEEE Trans. Med. Imag., vol. 28, no. 1, pp. 106–121, Jan. 2009. [19] Q. Ling and Z. Tian, “Decentralized sparse signal recovery for compressive sleeping wireless sensor networks,” IEEE Trans. Signal Process., vol. 58, no. 7, pp. 3816–3827, July 2010. [20] Z. Tian and G. B. Giannakis, “Compressed sensing for wideband cognitive radios,” in Proc. 2007 IEEE ICASSP Conf., pp. 1357–1360. [21] Y. Polo, Y. Wang, A. Pandharipande, and G. Leus, “Compressive wideband spectrum sensing,” in Proc. 2009 IEEE ICASSP Conf., pp. 2337– 2340. [22] J. A. Bazerque and G. B. Giannakis, “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,” IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1847–1862, Mar. 2010. [23] Y. Wang, Z. Tian, and C. Feng, “A two-step compressed spectrum sensing scheme for wideband cognitive radios,” in Proc. 2010 IEEE GLOBECOM, pp. 1–5. [24] F. Zeng, C. Li, and Z. Tian, “Distributed compressive spectrum sensing in cooperative multihop cognitive networks,” IEEE J. Sel. Topics Signal Process., vol. 5, no. 1, pp. 37–48, Feb. 2011. [25] Spectrum Policy Task Force, “Spectrum policy task force report,” Federal Communications Commission ET Docket 02-135, 2002. Yue Wang (M’11) received the B.S. degree in electrical engineering from Beijing Union University, Beijing, China, in 2004, and the Ph.D. degree in communication and information systems from Beijing University of Posts and Telecommunications, Beijing, China, in 2011. From September 2009 to March 2011, he was with the Department of Electrical and Computer Engineering, Michigan Technological University, Houghton, MI, USA, as a visiting Ph.D. student. In November 2011, he joined the Research Department of Hisilicon, Huawei Technologies Co., Ltd., Beijing, China, where he is currently a Research Engineer. His general interests are in the areas of wireless communications and signal processing. Current research focuses on cognitive radio and compressive sampling.
2125
Zhi Tian (M’98, SM’06) received the B.E. degree in Electrical Engineering from the University of Science and Technology of China, Hefei, China, in 1994, the M.S. and Ph.D. degrees from George Mason University, Fairfax, VA, in 1998 and 2000. Since August 2000, she has been on the faculty of Michigan Technological University, where she is currently a Professor. Dr. Tian’s general interests are in the areas of signal processing for wireless communications, estimation and detection theory. Current research focuses on cognitive radio networks and distributed wireless sensor networks. She served as Associate Editor for IEEE T RANSACTIONS ON W IRELESS C OMMUNICATIONS and IEEE T RANSACTIONS ON S IGNAL P ROCESSING. She received a CAREER award from the US National Science Foundation in 2003. Chunyan Feng received the B.S. degree in communications engineering, the M.S. and Ph.D. degrees in communication and information systems all from Beijing University of Posts and Telecommunications, Beijing, China. She is currently a Professor with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications. Her research interests are in the areas of broadband networks and wireless communication systems. Current research focuses on cognitive radio, key technology of B3G/4G systems, and green wireless communications.
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
Sparsity Order Estimation and its Application in Compressive Spectrum Sensing for Cognitive Radios Yue Wang, Member, IEEE, Zhi Tian, Senior Member, IEEE, and Chunyan Feng
Abstract—Compressive sampling techniques can effectively reduce the acquisition costs of high-dimensional signals by utilizing the fact that typical signals of interest are often sparse in a certain domain. For compressive samplers, the number of samples Mr needed to reconstruct a sparse signal is determined by the actual sparsity order Snz of the signal, which can be much smaller than the signal dimension N . However, Snz is often unknown or dynamically varying in practice, and the practical sampling rate has to be chosen conservatively according to an upper bound Smax of the actual sparsity order in lieu of Snz , which can be unnecessarily high. To circumvent such wastage of the sampling resources, this paper introduces the concept of sparsity order estimation, which aims to accurately acquire Snz prior to sparse signal recovery, by using a very small number of samples Me less than Mr . A statistical learning methodology is used to quantify the gap between Mr and Me in a closed form via data fitting, which offers useful design guideline for compressive samplers. It is shown that Me ≥ 1.2Snz log(N/Snz + 2) + 3 for a broad range of sampling matrices. Capitalizing on this gap, this paper also develops a two-step compressive spectrum sensing algorithm for wideband cognitive radios as an illustrative application. The first step quickly estimates the actual sparsity order of the wide spectrum of interest using a small number of samples, and the second step adjusts the total number of collected samples according to the estimated signal sparsity order. By doing so, the overall sampling cost can be minimized adaptively, without degrading the sensing performance. Index Terms—Sparsity order estimation, curve fitting, compressive sampling, cognitive radio, wideband spectrum sensing.
Manuscript received March 22, 2011; revised August 14, 2011; accepted September 23, 2011. The associate editor coordinating the review of this paper and approving it for publication was M. Morelli. This work was supported in part by the US National Science Foundation (NSF) under Grant ECS-0925881, the National Science Foundation of China (NSFC) under Grant 60902047, and the National High Technology Research and Development Program of China (“863” Program) under Grant 2009AA011802. Part of this work was presented at the IEEE GLOBECOM Conference, Miami, FL, USA, Dec. 2010. Y. Wang was with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing 100876, China. He is now with the Research Department of Hisilicon, Huawei Technologies Co., Ltd., Beijing 100095, China (e-mail: [email protected]). This work was performed when Y. Wang was a visiting Ph.D. student at Michigan Technological University. Z. Tian is with the Department of Electrical and Computer Engineering, Michigan Technological University, Houghton, MI 49931, USA (e-mail: [email protected]). C. Feng is with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing 100876, China (e-mail: [email protected]). Digital Object Identifier 10.1109/TWC.2012.050112.110505
I. I NTRODUCTION OMPRESSIVE sampling (CS) is an emerging framework for efficient acquisition and processing of sparse signals. It deals with the reconstruction of sparse signals from a small number of compressive samples collected from linear projections of the signals [1], [2]. Compared with the traditional Nyquist sampling theory, the number of samples required in CS techniques can be made much smaller than the signal dimension by exploiting the signal sparsity property exhibited in a certain domain. An Snz -sparse signal s ∈ RN can be expressed as a linear projection on some proper representation basis Ψ ∈ RN ×N in the form of s = Ψθ, where the vector of expansion coefficients θ ∈ RN only has Snz nonzero entries with Snz N . Essentially, the task of CS is to recover θ from a reduced set of M (< N ) samples x ∈ RM that are acquired by undersampling s with a linear compression matrix Φ ∈ RM×N , as follows:
C
x = Φs = ΦΨθ
(1)
where the number of rows in Φ decides the number of compressive samples collected. Theoretical work indicates that θ can be recovered exactly from x when the multiplication of Φ and Ψ satisfies the so-called restricted isometry property (RIP) [1]–[3]. The RIP offers a lower bound on the minimum number of samples needed for sparse signal reconstruction. For example, for Gaussian and Bernoulli random matrices, a log-like expression M ≥ αSnz log (N/Snz ) holds for some constant α, which means that the number of samples needed for recovering θ and then s can be much smaller than the number of unknowns N when Snz N [3]. As a result, the minimum sampling rate R = RN q M/N for signal reconstruction is decided by Snz , given N and the Nyquist rate RN q . The smaller the sparsity order Snz is, the smaller M and R need to be in order to recover the signal precisely. However, in practice, the actual sparsity order Snz is usually unknown a priori or even dynamically changing; instead, the maximum sparsity order Smax can be made available as a statistical upper bound of Snz . Thus, the practical sampling rate has to be set conservatively based on Smax in lieu of Snz , which may result in an unnecessarily high sampling rate as Smax deviates from Snz . A remedy to this wastage of the sampling resources is by adopting sequential sampling techniques [4], [5]. In [4], one sample is taken at a time, and the sparse signal is recovered at each time based on batch processing of all the past samples. The
c 2012 IEEE 1536-1276/12$31.00
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
sampling process stops once the signal estimate converges. This method entails high computational costs, because batchform sparse signal recovery is performed at every time step. In [5], a recursive sparse signal recovery algorithm is developed to considerably reduce the computational cost compared with batch processing; yet, it aims to reduce the computation per step and track slowly-varying signals, rather than minimizing the sampling rate and sensing costs. To alleviate the wasteful sampling costs at low reconstruction load, this paper aims to learn the signal sparsity order on the fly by introducing a novel concept of sparsity order estimation. The sparsity order estimation problem has not been studied in the CS literature. Indeed, most work on the minimum number of samples needed in CS can be divided into two categories: one for sparse signal estimation and the other for sparse support estimation. The former seeks to exactly reconstruct the sparse signal including both the positions and the amplitudes of all nonzero elements [1]– [3], while the latter aims to identify the positions of these nonzero elements [6], [7]. Evidently, both estimation tasks demand more sampling resources than the sparsity order estimation problem, because the latter only needs to identify the number of nonzero elements in the sparse signal, but not the position or the amplitude of each nonzero element. The tenet of CS suggests the possibility that signals be sampled at a rate proportional to the amount of information content to be extracted, which alludes to a lower sampling rate for sparsity order estimation. In this paper, a general methodology is provided to quantify the minimum number of samples required to estimate the sparsity order for different types of sampling matrices. To get around the difficulty in theoretical order statistics analysis for the sparsity order estimation problem, a simulation-based approach coupled with the curve fitting technique is employed to depict the fundamental limits of sampling rates in fitted closed-form expressions. The study confirms that the sparsity order estimation problem requires a smaller number of samples than that for sparse signal or support estimation problems. The nontrivial gap in terms of the sampling resources needed for sparsity order estimation versus sparse signal reconstruction can be useful in signal processing applications. This paper offers an illustrative example by developing a twostep compressive spectrum sensing (TS-CSS) algorithm for wideband cognitive radios (CRs). The TS-CSS first estimates the sparsity order Snz of the unknown wideband spectrum using a small number of samples collected at the first step, and then uses the estimated Sˆnz to decide the number of additional samples to be collected at the second step. Afterward, all samples from both steps are used to accurately reconstruct the wideband spectrum and eventually make the final decision on the spectrum occupancy. Overall, the TS-CSS adaptively adjusts its sampling rate according to Sˆnz estimated on the fly. As a result, the average sampling rate of the TS-CSS can be considerably lower than the traditional compressive sampling rate determined conservatively by Smax . Such a saving in acquisition costs is confirmed via simulations, given the desired sensing accuracy. Noticeably, reducing the sampling costs in CS applications can also lower the computational burden of the sparse signal recovery operations, which reduces the overall
2117
(hardware) implementation costs in practice. The rest of this paper is organized as follows. Section II describes the general principle to quantify the gap between the minimum number of samples needed for sparsity order estimation and that for sparse signal reconstruction. Section III illustrates the usefulness of such a gap by developing a twostep compressive spectrum sensing algorithm for wideband CRs. Performance evaluation results are presented in Section IV, followed by conclusions in Section V. II. F UNDAMENTAL L IMITS IN S PARSITY O RDER E STIMATION This section starts by introducing a novel concept of sparsity order estimation, which aims to identify the number of nonzero elements of a sparse signal without having to know their positions or amplitudes. To validate this concept, an empirical study is set up to evaluate the minimum number of samples needed for accurate sparsity order estimation. Then, a general methodology based on the curve fitting technique is developed to delineate the closed-form expressions for several widely-used random sampling matrices. This study provides quantitative results on the nontrivial gap between the minimum number of samples required for sparsity order estimation and that for sparse signal reconstruction. A. Concept of Sparsity Order Estimation The sparsity order of an N × 1 sparse signal vector θ is defined by its 0 -norm in the form Snz = θ0 , which refers to the number of nonzero elements of θ. The sparsity order estimation problem concerns estimating Snz from the compressive samples x in (1). The number of samples needed for sparsity order estimation depends on not only N and Snz , but also the matrices Φ and Ψ, which are discussed next. B. Options of Matrix Patterns The sample quality and estimation performance are determined by the matrices used in the sampling process. Here several widely-used matrix patterns are reviewed for later experimental use. 1) Sparse Representation Basis Ψ: In the context of CS, a signal may have a sparse representation when it is expressed on a proper basis characteristic of its specific application, e.g., Fourier basis for signals with underutilized frequency spectrum [8], wavelet basis and discrete curvelet transform basis for natural image signals [1], [9], etc. This paper will study frequency-sparse signals as an illustrative example in underutilized wideband CR networks; correspondingly the focus is on the Fourier basis in the form of an N × N discrete Fourier transform (DFT) matrix. For testing the universality of our methodology, the identity matrix will also be studied as a special case of the representation basis, which is employed when the signal sparsity appears in the same domain of the sampling process. 2) Random Sampling Matrix Φ: It has been shown that Gaussian and Bernoulli random sampling matrices have a high probability of satisfying the RIP property [3]. It means that samples acquired from such matrices have a better chance to
2118
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
perfectly recover sparse signals than other types of sampling matrices, given the same number of samples M . The entries {φi,j }i,j of these two random sampling matrices are independent and identically distributed according to Gaussian and Bernoulli probability distributions as follows (w.p. stands for with probability): Gaussian: Bernoulli:
φi,j ∼ N (0, 1/M ) ; √ +1/ M w.p. 1/2, √ φi,j = −1/ M w.p. 1/2.
(2) (3)
There are several options for hardware implementations of random sampling matrices [10]. In practical settings, structured random matrices are often employed for improved implementation affordability [11]–[13]. An exemplary sampler that will be tested in this paper is the so-called analogto-information converter (AIC) [11], which compresses and samples analog signals at uniform sub-Nyquist sampling rates. It has a serial concatenated structure with three main components: a wideband modulator made of a pseudo-random chipping sequence, a low-pass filter, and a uniform low-rate conventional analog-to-digital converter (ADC). The resulting sampling matrix Φ turns out to be a banded random matrix with a great deal of structure. Because of the special structure, the AIC alleviates the sampling burden on the ADC component, at the expense of slightly degraded recovery performance compared with those fully random and dense Gaussian or Bernoulli samplers. C. Experimental Design We now present a general methodology to determine the minimum number of samples Me that is required to reliably identify the sparsity order Snz of a sparse signal θ, while that needed to reconstruct the signal itself Mr is also provided as a benchmark. The developed approach is simulation-based, in which a large set of Monte Carlo simulation trials is performed following the steps below, and the trial results are used to fit into closed-form expressions quantifying Me and Mr as functions of the signal size N and the actual sparsity order Snz , for a given measurement matrix ΦΨ. Simulation-based quantification is a valid statistical learning tool that has been theoretically justified and fruitfully adopted in e.g., Markov Chain Monte Carlo (MCMC) methods for Bayesian networks [14]. In CS research, [11] adopts a simulation-based approach to quantify the performance bound of practical AIC samplers. Even when the RIP condition is derived analytically to come up the well-known expression Mr ≥ αSnz log(N/Snz ) for Gaussian and Bernoulli sampling matrices, the theoretical analysis is generally lacking concerning the precise value of α [1]–[3]. If a specific value for α is required, then one typically must determine it experimentally. Besides, theoretical results tend to work in an asymptotic regime where the probabilistic characteristics of the random matrices are valid, and hence are useful when the problem size is very large but may appear to be too conservative for practical problems of reasonably small to medium sizes. For these reasons, this paper develops experimental design along with the curve fitting technique to bypass the difficulty in theoretically analyzing the sparsity order estimation problem. The developed methodology results
in usable expressions that fit the simulation data into a closed form, and is general enough to be readily applied offline for a wide range of sampling matrices under practical operating regimes, without incurring high costs in computation and analysis. The following steps are performed in each Monte Carlo test. 1) Generate a Sparse Signal: An Snz -sparse signal vector θ of size N is generated randomly according to a generic stochastic model. That is, given a preset integer Snz ( N ), the sparse support of θ is a uniformly distributed random set of the Snz nonzero positions as follows: Supp (θ) = {{k1 , . . . , kSnz } |θki = 0, ki ∈ {1, . . . , N }} where
Pr(ki ∈ Supp(θ)) = Snz /N,
i = 1, . . . , N.
(4)
Meanwhile, the amplitude of each element of θ is set according to one of the following commonly-used distributions:1 = 1, ki ∈ Supp (θ) , Binary: θki (5) = 0, ki ∈ / Supp (θ) ; ∼ N (1, 0.01) , ki ∈ Supp (θ) , Gaussian: θki (6) = 0, ki ∈ / Supp (θ) ; ∼ U (0.8, 1.2) , ki ∈ Supp (θ) , or Uniform: θki (7) = 0, ki ∈ / Supp (θ) . 2) Collect Compressive Samples: With the randomly generated θ and the chosen Ψ and Φ, M measurements are generated from (1) to form a length-M sample vector x. ˆ of the 3) Reconstruct the Sparse Signal: An estimate θ original θ is computed using the basis pursuit (BP) algorithm for the noise-free case, which offers a tractable solution by minimizing the 1 -norm of θ subject to a linear projection constraint [1], as follows: ˆ = arg min θ s.t. x = ΦΨθ. θ 1
(8)
For the noisy case, the least absolute shrinkage and selection operator (LASSO) algorithm can be used to yield [15] ˆ = arg min θ s.t. x − ΦΨθ ≤ θ 1 2
(9)
where bounds the amount of noise energy in the sample vector. 4) Estimate the Sparsity Order: Having obtained the signal ˆ the sparsity order Snz can be estimated as estimate θ, Sˆnz =
N ˆ θ [i] ≥ λ
(10)
i=1
where λ > 0 is a preset decision threshold. Empirically, the threshold is set as λ = μ/2 for the noise-free case or λ = (μ + σn ) /2 for the noisy case, where μ is the average absolute value of all the Snz nonzero elements in θ, σn is the standard deviation of the noise component, and λ lies in the middle of μ and σn . As Fig. 1 indicates, the amplitudes of the ˆ are well separated from those of the nonzero elements in θ 1 This paper focuses on these simple and widely-used probabilistic distributions, while the amplitude setting can be extended to other distributions as well. The choice of the amplitude distribution has little impact on the simulation-based approach, because θ can be recovered accurately for any amplitudes, once the nonzero support is correctly identified.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
70
0.5
60
0
0
5
10
15 (a)
20
25
30
1 0.5 0
0
5
10
15 (b)
20
25
30
Minimum number of samples
1
50
40 min Me (Gaussian) min Me (Bernoulli) min Me (AIC) min Mr (Gaussian) min Mr (Bernoulli) min Mr (AIC) 1.2Slog(N/S+2)+3.2 1.32Slog(N/S+4)+2.6 1.71Slog(N/S+1)+0.9
30
20
1
=0.5
0.5 0
10
0 0
5
10
15 (c)
20
25
30
1
0
0
5
10
15 (d)
20
0
5
10 Sparsity order
15
20
Fig. 2. Minimum Me and Mr vs. Snz for sparsity order estimation and sparse signal reconstruction with different sampling matrices, using DFT representation basis and N = 128 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
=0.5
0.5
2119
25
30
D. Quantification Methodology Fig. 1. Illustration of a noise-free case when N = 30, Snz = 5, and λ = 0.5. (a) The sparse input signal generated from (4) and (5), (b) a successful reconstruction of the input signal, (c) a successful estimation of both the nonzero support and the sparsity order, although the amplitudes are not recovered accurately, and (d) a successful estimation of sparsity order, even though errors occur in both the amplitudes and the nonzero support.
zero elements, which makes the decision in (10) robust to the threshold value λ. This is because the 1 -norm penalty term on θ enforces a sparse solution by shrinking small estimates to zero. With (10), a successful identification of the sparsity order is declared when the estimated sparsity order is equal to the true value, i.e., Sˆnz = Snz . To identify the minimum number of samples Me needed for reliable sparsity order estimation, a range of values for M is tested, and Me is given by the smallest M that results in a success rate of 99% in correctly identifying Snz among 500 Monte Carlo trials. The simulation-based approach can also be used to empirically observe the minimum number of samples Mr needed to accurately reconstruct the sparse input signal. In that case, ˆ = θ to machine a successful recovery is declared when θ precision, and the smallest value Mr is reported at which the probability of empirical successes is greater than 99% in 500 Monte Carlo trials. Fig. 1 illustrates several simulation outcomes of the sparse signal recovery problem in Step 3) and the sparsity order estimation problem in Step 4). Apparently, given the same sampling resources, the probability of successful sparsity order estimation can be higher, because Sˆnz = Snz is possible even when the amplitudes or positions of those nonzero elements are recovered inaccurately.
According to the experimental design described above, we now analyze the empirical results on the minimum numbers of samples Me and Mr for the noise-free case. Specifically, we seek to delineate Me and Mr as functions of the sparsity order Snz and the signal size N , given the sampling matrix pattern. As suggested by both analysis [2], [3] and experiments [11], the desired function for Mr typically takes on a log-like form, which also conforms well with the empirical data we have observed for both Mr and Me . Hence, our quantification methodology is to adopt the general expression ASnz log(N/Snz + B) + C, and use experimental results and statistical curve fitting techniques to find the values of A, B and C that minimize the data fitting error. First, we fix N and find the expressions of Me and Mr as functions of Snz . Consider for example the case of N = 128 and using the DFT representation basis, while Snz varies from 1 to 20. For different sampling matrices including Gaussian, Bernoulli and AIC, the empirical data of Me and Mr as well as the fitted expressions using statistical regression are depicted in Fig. 2. The fitted expression of the minimum Me for sparsity order estimation is given by Me,min ≈ 1.2Snz log (N/Snz + 2) + 3.2. Meanwhile, the fitted expressions of the minimum Mr for sparse signal recovery are in the form of Mr,min ≈ 1.32Snz log (N/Snz + 4) + 2.6 for Gaussian and Bernoulli sampling matrices, and Mr,min ≈ 1.71Snz log (N/Snz + 1) + 0.9 for the AIC sampling matrix [11]. Apparently, in order to achieve the same reconstruction accuracy, a structured AIC random sampler requires a bit more samples than samplers using Gaussian and Bernoulli random matrices, in return for simplified hardware implementation. Alternatively, we fix Snz and quantify the two numbers Me and Mr as functions of N . Fig. 3 shows the empirical results and the fitted curves for Snz = 10, where N increases from 64 to 1024. The fitted expressions from both Fig. 2 and Fig. 3 are near identical.
2120
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
100
80 70 60
90 80 Minimum number of samples
Minimum number of samples
90
100 min Me (Gaussian) min Me (Bernoulli) min Me (AIC) min Mr (Gaussian) min Mr (Bernoulli) min Mr (AIC) 1.2Slog(N/S+2)+2.9 1.32Slog(N/S+4)+3.5 1.69Slog(N/S+1)+0.3
50 40
60 50 40 30
30 20
70
min Me (Gaussian) min Mr (Gaussian) 1.5Slog(N/S+2)+6.6 2.5Slog(N/S−1)+10.9
20
2
10
3
10
10 Signal dimension
Fig. 3. Minimum Me and Mr vs. N for sparsity order estimation and sparse signal reconstruction with different sampling matrices, using DFT representation basis and Snz = 10 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
E. Fundamental Limits in Sample Rates The following Results (1-3) summarize the curve fitting results, from the two sets of experiments depicted in Fig. 2 and Fig. 3. Result 1: Successful estimation of the sparsity order can be achieved with 99% probability when the number of samples Me satisfies the lower bound: Me ≥ 1.2Snz log (N/Snz + 2) + 3.
(11)
Result 2: To obtain successful sparse signal reconstruction with 99% probability, the lower bounds of Mr for different samplers are Mr ≥ 1.32Snz log (N/Snz + 4) + 3
(12)
for Gaussian/Bernoulli samplers, and Mr ≥ 1.7Snz log (N/Snz + 1)
(13)
for the AIC sampler [11]. Comparison of Result 1 and Result 2 leads to Result 3. Result 3: For any given Snz and N , it always holds that Mr ≥ Me . The expressions (11)-(13) are concluded when a size-N DFT matrix is used as the sparse representation basis Ψ. This basis matrix is relevant to the compressive spectrum sensing problem to be investigated in Section III. The general principle of statistical learning can be used to quantify the fundamental sampling costs for any matrix patterns. For example, when Ψ = I is an identity matrix, the empirical results and fitted curves are illustrated in Fig. 4. The required sampling costs appear to be higher than those in (11)-(13), even though the general trends are quite similar. In general, the proposed simulation-based approach can be conducted offline for any choices of matrix patterns, and offers simple expressions to guide practical sampling systems.
0
5
10 Sparsity order
15
20
Fig. 4. Minimum Me and Mr vs. Snz for sparsity order estimation and sparse signal reconstruction with Gaussian random matrix as sampler and identity matrix as the sparse representation basis, for N = 128 in noise-free case. The markers show the empirical results, while the curves delineate the curve fitting results.
III. A PPLICATION IN C OMPRESSIVE S PECTRUM S ENSING FOR W IDEBAND CR S As a new paradigm for efficient acquisition and processing of sparse signals, CS has found attractive applications in a broad range of signal processing problems such as image processing [16], [17], medical imaging [18], wireless sensor networks [19], cognitive spectrum sensing [20]–[24], to name a few. As an illustrative example, this section presents the application of sparsity order estimation in compressive spectrum sensing for wideband CRs. In CR networks where users dynamically share the spectrum, CRs need to rapidly and accurately identify spectrum opportunities over a very wide band. Subject to analog-todigital hardware limitations in the over-GHz regime, CS techniques are well motivated for wideband spectrum detection exploiting sparsity, because there are plenty of temporarily unused spectrum opportunities in wideband CR networks [25]. Several CS-based approaches have been developed to detect the frequency occupancy of primary users using sub-Nyquistrate samples [20], [21]. To determine the suitable sampling rate, they all implicitly assume that (the upper bound of) the sparsity order of the underutilized spectrum is known exactly. However, in practical CR applications, the actual sparsity order Snz corresponds to the instantaneous spectrum occupancy of wireless users which is dynamically varying. Hence, Snz is often unknown; instead, what is typically available is its upper bound Smax , which can be measured from the maximum spectrum utilization observed statistically over a time period. In practice, conservative determination of the sampling rate based on Smax can cause unnecessarily high sensing costs. To alleviate such wasteful acquisition costs and keep low reconstruction load, we develop a two-step compressive spectrum sensing (TS-CSS) algorithm for wideband CRs [23], which capitalizes on our key observation of the gap between the number of samples needed for sparsity order estimation and that for sparse signal reconstruction.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
A. System Model and Problem Statement
70
(14)
where xt of length M consists of the collected compressive samples, rt of length N is the discrete-time representation of r(t) at Nyquist rate with elements rt [n] = r (t) | t=n/RN q , n = 0, . . . , N − 1, F−1 is the N × N inverse DFT matrix, and Φ is the M × N (M < N ) AIC sampling matrix. Evidently, the actual sampling rate R = RN q M/N in collecting the sample vector xt is lower than the Nyquist rate RN q , while the recovery of the unknown rt or rf in (14) will enable spectrum sensing at full frequency resolution. B. Two-Step Compressive Spectrum Sensing (TS-CSS) In the presence of an unknown Snz , a TS-CSS algorithm is now developed for wideband spectrum sensing at low sampling costs. The basic idea is to utilize the nontrivial gap between the minimum numbers of samples for sparsity order estimation and for sparse signal reconstruction (c.f. Fig. 5). To do so, each sensing block is divided into two time slots of length T1 and T2 , used for sparsity order estimation and signal detection, respectively. A typical choice is to set T1 = T2 = T /2, while the effect of non-equal time division will be discussed later. The overall procedure of the TS-CSS algorithm is shown in Fig. 6, and details of each step are elaborated next.
Mc
60
Minimum number of samples
Suppose that a wideband primary network spans over a total of B Hz, which consists of N non-overlapping subchannels of equal bandwidth. At a given time and location, the network is underutilized with only Snz ( N ) channels occupied by active primary users, whereas the rest N − Snz channels are idle and offer spectrum opportunities to secondary CR users. Clearly, the sparse nature of the signal of interest exhibits in the frequency domain. However, neither the locations of the occupied Snz channels, nor the sparsity order Snz , are known a priori. Instead, only an upper bound Smax of the sparsity order can be available in practice. Given the received signal r (t) of bandwidth B Hz and a coarsely known Smax , the wideband spectrum sensing task for a CR detector is to detect the occupancy of the N spectrum bands. Those idle frequency bands present potential transmission opportunities for CRs, if and only if they can be detected correctly. A CR detector adopting CS techniques starts with an analog-to-digital sampler that collects M discrete-time samples from r (t). Here we utilize the AIC as the random sampler and the DFT matrix as the sparse representation basis. During spectrum sensing, each time window for sensing shall be small enough compared with the channel coherence time and user dynamics, such that transmission opportunities can be utilized before the spectrum occupancy profile has changed. It is hence justified to assume that the unknown frequency spectrum of interest is block stationary for each sensing window of time length T = N/RN q , given N and RN q . Spectrum sensing is performed block by block, aiming to recover the N × 1 unknown frequency response vector rf of invariant strengths and Snz in each block. Hence, the discrete-time representation of the CS sampling operation on the received r (t) at a CR detector can be expressed in the following general form: xt = Φrt = ΦF−1 rf
2121
M c M
Saving
50
^
M
`M
M M1
2
M1
For sparse signal reconstruction
40
30
20
For sparsity order estimation
10
0
S nz 0
2
4
6
8 10 12 Sparsity order
S max 14
16
18
20
Fig. 5. Basic idea of the TS-CSS utilizing the gap between the numbers of samples required for sparsity order estimation at the first step and for sparse signal reconstruction at the final step.
r (t )
Estimate sparsity order
M1
Reconstruct sparse signal
S nz
rˆ f
Make sensing decision
d
M2
J
S max
Fig. 6.
Predict sampling rate
Flow chart of the TS-CSS algorithm.
1) Estimate the Sparsity Order at the First Step: With the available Smax and using the solid curve for Me in Fig. 5, the first step is to decide M1 which is the number of samples needed for estimating the sparsity order within the time slot T1 . From Result 1, and substituting the unknown Snz in (11) by its upper bound Smax , M1 is decided as follows: M1 = 1.2Smax log (N/Smax + 1/2) .
(15)
Here, M1 is chosen conservatively according to Smax , such that it suffices to accurately estimate the sparsity order of any signal with Snz ≤ Smax . Since the M1 samples are acquired during the first half period of r(t), the sampling process at the first step can be described as xt1 = [Φ1 0] rt = [Φ1 0] F−1 rf
(16)
where Φ1 is the M1 × N/2 sampling matrix. Given xt1 , the TS-CSS aims to estimate Sˆnz of the sparse spectrum rf . This can be done by first recovering ˆrf via (8), where θ is replaced by rf and (16) is employed as the linear
2122
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
constraint. Then, following (10), Sˆnz is estimated as Sˆnz =
N
(|ˆrf [i]| ≥ λ) .
Algorithm: Two-Step Compressive Spectrum Sensing (TS-CSS) (17)
i=1
2) Sample Adaptively at the Second Step: At the second step of the TS-CSS, additional M2 samples are collected into a data vector xt2 over the second time period of r(t), using an M2 × N/2 sampling matrix Φ2 as follows: xt2 = [0 Φ2 ] rt = [0 Φ2 ] F−1 rf .
(18)
The primary issue at the second step is to decide the number M2 such that rf can be accurately estimated from a total of M = M1 + M2 compressive samples. This can be done by using the dashed curve in Fig. 5 along with the estimated Sˆnz . Specifically, replacing Snz in (13) by Sˆnz , M can be decided as M = 1.7Sˆnz log(N/Sˆnz + 1).
(19)
If the upper bound Smax is loose such that M1 ≥ M , then the M1 samples collected at the first step are already adequate for reliably recovering the sparse spectrum, and hence there is no need to sample further at the second step. Thus, M2 is decided as M2 = max(M − M1 , 0)
(20)
where M1 and M are calculated from (15) and (19), respectively. 3) Reconstruct the Signal Spectrum from All the Collected Samples: After the two-step sampling process, the unknown wideband spectrum rf can be reconstructed from all available samples using the 1 -norm minimization formula in (8). Concatenating the sample vectors xt1 and xt2 in (16) and (18) from both steps, (8) becomes ˆrf = arg min rf 1 s.t.
xt1 xt2
=
Φ1 0
0 Φ2
F−1 rf .
(21) 4) Make Spectrum Occupancy Decisions: Having obtained ˆrf from (21), the CR detector can finally make decisions on the spectrum occupancy by comparing the signal energy within each frequency band with a decision threshold γ, as follows: ˆ = |ˆrf |2 ≥ γ . d
(22)
ˆ of length N indicates the sensing The binary decision vector d outcomes on the monitored wideband spectrum. If an entry of ˆ is 0, then the corresponding frequency band is decided to d be idle, presenting a potential spectrum opportunity for CRs. If the entry is 1, then the corresponding frequency band is decided to be occupied by primary users. The above steps are summarized in Algorithm TS-CSS.
• • • • • •
set T1 = T2 = T /2; collect M1 samples from r(t) over the first time period of T1 seconds, where M1 is decided by (15); estimate Sˆnz from (17); collect M2 samples from r(t) over the second time period of T2 seconds, where M2 is decided by (20); estimate the sparse signal spectrum ˆrf from (21); ˆ from (22). make the spectrum occupancy decision d
C. Discussions Several remarks on this TS-CSS algorithm are in order. 1) Saving in Sampling Costs: It is of interest to compare the sampling cost of the TS-CSS with that of a traditional onestep CS approach, which would directly use Smax to decide the total number of samples M used for recovering the sparse signal within one step. Essentially, M is given by Result 2 after replacing the unknown Snz by its upper bound Smax , which is indicated in Fig. 5. Apparently, M exceeds the total number M used in the TS-CSS, resulting in higher samplingrate requirements for the one-step approach. The difference (M − M ) reflects the saving in sampling costs. Notice that if the prior knowledge is already very good such that Smax = Snz , then the saving decreases to 0. On the other hand, given a certain Snz , the largest amount of saving happens when M1 ≥ M , in which case there is a big gap between the prior knowledge and the actual sparsity order such that Smax Snz . At least, the saving is positive as long as Smax > Snz , which takes place in most cases. The quantitative saving in sampling costs will be presented in next section via computer simulations. 2) Time Ratio T1 /T : The ratio T1 /T may impact the sensing efficiency. Intuitively, an assignment that consumes fewer samples and shorter length in the first time slot results in more efficient sensing, since Snz can be quickly estimated to adaptively decide the number of additional samples required in the second time slot. Fig. 7 presents the minimum numbers of samples required in the first time slot for different values of the time ratio T1 /T . The observed data and corresponding fitted curves therein are obtained by the same methodology developed in Section II. Fig. 7 indicates that the lower the ratio is chosen, the fewer samples are needed to successfully estimate the sparsity order. This result can be explained by the higher peak sampling rate resulted from a shorter length of the first time slot, as shown in Fig. 8. However, the drawback of a small ratio is the limited range of sparsity order which can be estimated; for example, the range of predictable Snz is only from 1 to 12 for the ratio T1 /T = 1/4 in Fig. 8, which also explains why the curve corresponding to the ratio T1 /T = 1/4 has a shorter range than its counterparts in Fig. 7. Therefore, there exists a tradeoff between the sampling efficiency and the applicable range of the TS-CSS algorithm. For simplicity and without loss of generality, T1 /T = 1/2 is adopted. 3) Implementation of Sampling Rate Switching: In TSCSS, there are two different sampling rates for the two successive steps, each of which is a uniform sub-Nyquist rate. The rate switching can be conveniently incorporated in practical
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
80
Minimum number of samples
70 60 50 40 T1/T=1 T1/T=3/4 T1/T=1/2 T1/T=1/4 1.2Slog(N/S+2)+3 1.2Slog(N/S+1.5)+2.9 1.2Slog(N/S+0.5)+0.1 1.2Slog(N/S−1.5)+0.2
30 20 10 0
0
5
10
15 20 Sparsity order
25
30
35
Fig. 7. Minimum number of samples vs. sparsity order for sparsity order estimation, for different values of the ratio T1 /T when N = 128 in noise-free case.
1 0.9
Peak sampling rate / RNq
0.8 0.7 0.6
T1/T=1 T1/T=3/4 T1/T=1/2 T1/T=1/4 (1.2Slog(N/S+2)+3)/N 4/3(1.2Slog(N/S+1.5)+2.9)/N 2(1.2Slog(N/S+0.5)+0.1)/N 4(1.2Slog(N/S−1.5)+0.2)/N
0.4 0.3 0.2 0.1 0
Consider a wideband CR system with N = 1024, Smax = 20, and Snz = 15. The sparse nonzero support is chosen as in (4), and the amplitudes of the nonzero elements are Gaussian distributed as in (6). The AIC sampler is used and the LASSO algorithm is employed for sparse signal recovery in the presence of ambient noise. The key design parameters are: M1 is computed from (15), M is computed from (19), and M2 = M − M1 . Sensing performance metrics of interest are the probability of miss-detection Pm and that of false-alarm Pf . These two detection criteria indicate the level of interference to primary users and the wastefulness of available spectrum opportunities for CRs, respectively. Given the true state vector d ∈ {0, 1}N of the target wideband spectrum, Pm and Pf can be expressed as follows: ˆ ˆ dT (d = d) (1 − d)T (d = d) , Pf = E Pm = E 1T d N − 1T d (23) where E denotes expectation, and 1 denotes the length-N allone vector. Practical CR applications, e.g., IEEE 802.22, concern the quality of service (QoS) of coexisting primary systems. To this end, it is required that Pm should be no higher than a predefined QoS-specific threshold Pm,th , given a tolerable Pf . Hence, the sensing performance can be evaluated by the following probability: Pr (Pm ≤ Pm,th ) , given a fixed Pf .
0.5
0
5
10
15
20 25 Sparsity order
30
35
40
45
Fig. 8. Peak sampling rate vs. sparsity order for sparsity order estimation, for different values of the ratio T1 /T when N = 128 in noise-free case.
hardware implementations, without causing much complication. For example, a serial random sampler, such as the AIC sampler [11], can flexibly adjust its sampling rate at the last ADC block to yield the two different effective sampling rates for the two steps, specifically, R1 = RN q M1 /(N/2) at the first step, and R2 = RN q M2 /(N/2) at the second step. Alternatively, in a parallel random sampler implementation for general sampling matrices [10], a punctured switch can be used to selectively admit different numbers of parallel branches in the two steps, specifically M1 and M2 branches, respectively. IV. P ERFORMANCE E VALUATIONS This section presents simulation results to verify the sensing performance and sampling costs of the proposed TS-CSS, using the traditional one-step CS approach as a benchmark.
2123
(24)
In (24), the two sensing performance metrics Pm and Pf are jointly considered to reflect the best-effort QoS satisfaction. In the simulations, both Pf and Pm,th are set to be 0.1, which means that the achieved probability of detection shall exceed 90%. Fig. 9 depicts the sensing performance versus signal-tonoise ratio (SNR) in terms of the best-effort QoS satisfaction. The SNR is defined to be the total signal power of the wideband signal over the entire spectrum, scaled by the power of the white noise. For a given SNR, the total number of samples used is averaged over 500 Monte Carlo trials and marked along the curves to show the sampling costs. In Fig. 9, it is clear that the best-effort QoS satisfaction of the TS-CSS is very close to that of the traditional one-step CS approach, when both of them use the same prior knowledge of Smax . As SNR increases, both of the probabilities of QoS satisfaction approach 1. Meanwhile, the total number of samples used in the TS-CSS reduces to M = 108 corresponding to an average sampling rate of R = (108/1024)RN q , whereas the traditional one-step CS approach samples at a higher rate of R = (135/1024)RN q throughout. The reduction in samples means lower sensing costs and less energy consumption, both of which are important merits for CR detectors in practice. Fig. 10 evaluates the saving rate of the acquisition costs offered by the TS-CSS over the one-step approach, for various Snz and Smax . The saving in sampling costs is measured by the following rate: η=
M − max (M1 , M ) × 100%. M
(25)
2124
IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 6, JUNE 2012
V. C ONCLUSIONS
1
108 114.7 109.7 116.8 110.1 108 108 117.2
0.9 0.8
118.7 135
Pr(Pm Pm,th)
0.7 0.6
119.5
0.5 0.4 119.7 0.3 0.2 119.9 0.1
traditional one−step CS proposed TS−CSS
123.7121.9 126.1 120.7
0 −6
−4
−2
0
2 SNR
4
6
8
10
Fig. 9. Comparison of the proposed TS-CSS and the traditional one-step CS approach in terms of the probability of QoS satisfaction and the total number of samples used for different SNR cases when N = 1024, Smax = 20, and Snz = 15.
The cornerstone of this work is the fact that the number of samples required for estimating the signal sparsity order can be much smaller than that for reconstructing the unknown sparse signal itself. To reveal such a nontrivial gap, this paper first presents a general methodology to quantify the gap via simulation-based probabilistic evaluation and the curve fitting technique. The fundamental limits in sparsity order estimation are delineated by data-fitted closed-form expressions with loglike form, for some widely-used sampling matrices. Next, to illustrate the usefulness of this gap in practical applications, a two-step compressive spectrum sensing (TS-CSS) algorithm is developed for wideband CRs that only coarsely know the spectrum sparsity order. The TS-CSS algorithm adaptively switches the sampling rates for the two sequential steps of sparsity order estimation and sparse spectrum reconstruction, which effectively avoids collecting unnecessary samples caused by the discrepancy between actual and maximum sparsity orders. Compared with the traditional one-step CS-based spectrum sensing, the TS-CSS utilizes the smallest possible number of samples based on the estimated sparsity order on the fly, and can hence considerably reduce the sampling costs while achieving the desired sensing accuracy.
30 Smax=10 S Saving rate of sampling costs (%)
R EFERENCES
=20
max
25
Smax=30
20
15
10
5
0
0
5
10
15 Snz
20
25
30
Fig. 10. Saving rate η of sampling costs achieved by the proposed TS-CSS for different Snz and Smax when N = 1024.
As Fig. 10 indicates, the saving in sampling costs is most significant when Snz is far off from Smax . When Smax − Snz is large, the prior knowledge about sparsity order is coarse, which results in a large gap between the minimum number of samples for sparsity order estimation and that for sparse signal reconstruction; correspondingly, this is considerable room for saving the sampling costs, and the proposed TSCSS collects the benefit of this gap. The maximum saving in sampling costs stays unchanged when Snz Smax namely Snz ≤ Smax ×60% in Fig. 10. This is because the M1 samples collected at the first step in order to estimate Sˆnz are already sufficient to recover the sparse signal accurately. On the other hand, η approaches 0 as Snz gets close to Smax , which means diminished saving in sampling costs when the prior knowledge is close to the actual sparsity order.
[1] E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, Mar. 2008. [2] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, Apr. 2006. [3] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “A simple proof of the restricted isometry property for random matrices,” Const. Approx., vol. 28, no. 3, pp. 253–263, Dec. 2008. [4] D. M. Malioutov, S. Sanghavi, and A. S. Willsky, “Compressed sensing with sequential observations,” in Proc. 2008 IEEE ICASSP Conf., pp. 3357–3360. [5] D. Angelosante and G. B. Giannakis, “RLS-weighted lasso for adaptive estimation of sparse signals,” in Proc. 2009 IEEE ICASSP Conf., pp. 3245–3248. [6] M. J. Wainwright, “Information-theoretic limits on sparsity recovery in the high-dimensional and noisy setting,” IEEE Trans. Inf. Theory, vol. 55, no. 12, pp. 5728–5741, Dec. 2009. [7] W. Wang, M. J. Wainwright, and K. Ramchandran, “Informationtheoretic limits on sparse support recovery: dense versus sparse measurement,” IEEE Trans. Inf. Theory, vol. 56, no. 6, pp. 2967–2979, June 2010. [8] E. J. Candes, J. K. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006. [9] E. J. Candes and L. Demanet, “The curvelet representation of wave propagators is optimally sparse,” Comm. Pure Appl. Math., vol. 58, no. 11, pp. 1472–1528, Nov. 2005. [10] M. Mishali and Y. C. Eldar, “From theory to practice: sub-Nyquist sampling of sparse wideband analog signals,” IEEE J. Sel. Topics Signal Process., vol. 4, no. 2, pp. 375–391, Apr. 2010. [11] J. A. Tropp, J. N. Laska, M. F. Duarte, J. K. Romberg, and R. G. Baraniuk, “Beyond Nyquist: efficient sampling of sparse bandlimited signals,” IEEE Trans. Inf. Theory, vol. 56, no. 1, pp. 520–544, Jan. 2010. [12] J. Romberg, “Compressive sensing by random convolution,” SIAM J. Imag. Sci., vol. 2, no. 4, pp. 1098–1128, Dec. 2009. [13] X. Chen, Z. Yu, S. Hoyos, B. M. Sadler, and J. Silva-Martinez, “A subNyquist rate sampling receiver exploiting compressive sensing,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 58, no. 3, pp. 507–520, Mar. 2011. [14] C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd edition. Springer-Verlag, 2004. [15] R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. Roy. Statist. Soc. Ser. B, vol. 58, no. 1, pp. 267–288, 1996.
WANG et al.: SPARSITY ORDER ESTIMATION AND ITS APPLICATION IN COMPRESSIVE SPECTRUM SENSING FOR COGNITIVE RADIOS
[16] J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 14–20, Mar. 2008. [17] M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 83–91, Mar. 2008. [18] J. Trzasko, A. Manduca, and E. Borisch, “Highly undersampled magnetic resonance image reconstruction via homotopic 0 -minimization,” IEEE Trans. Med. Imag., vol. 28, no. 1, pp. 106–121, Jan. 2009. [19] Q. Ling and Z. Tian, “Decentralized sparse signal recovery for compressive sleeping wireless sensor networks,” IEEE Trans. Signal Process., vol. 58, no. 7, pp. 3816–3827, July 2010. [20] Z. Tian and G. B. Giannakis, “Compressed sensing for wideband cognitive radios,” in Proc. 2007 IEEE ICASSP Conf., pp. 1357–1360. [21] Y. Polo, Y. Wang, A. Pandharipande, and G. Leus, “Compressive wideband spectrum sensing,” in Proc. 2009 IEEE ICASSP Conf., pp. 2337– 2340. [22] J. A. Bazerque and G. B. Giannakis, “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,” IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1847–1862, Mar. 2010. [23] Y. Wang, Z. Tian, and C. Feng, “A two-step compressed spectrum sensing scheme for wideband cognitive radios,” in Proc. 2010 IEEE GLOBECOM, pp. 1–5. [24] F. Zeng, C. Li, and Z. Tian, “Distributed compressive spectrum sensing in cooperative multihop cognitive networks,” IEEE J. Sel. Topics Signal Process., vol. 5, no. 1, pp. 37–48, Feb. 2011. [25] Spectrum Policy Task Force, “Spectrum policy task force report,” Federal Communications Commission ET Docket 02-135, 2002. Yue Wang (M’11) received the B.S. degree in electrical engineering from Beijing Union University, Beijing, China, in 2004, and the Ph.D. degree in communication and information systems from Beijing University of Posts and Telecommunications, Beijing, China, in 2011. From September 2009 to March 2011, he was with the Department of Electrical and Computer Engineering, Michigan Technological University, Houghton, MI, USA, as a visiting Ph.D. student. In November 2011, he joined the Research Department of Hisilicon, Huawei Technologies Co., Ltd., Beijing, China, where he is currently a Research Engineer. His general interests are in the areas of wireless communications and signal processing. Current research focuses on cognitive radio and compressive sampling.
2125
Zhi Tian (M’98, SM’06) received the B.E. degree in Electrical Engineering from the University of Science and Technology of China, Hefei, China, in 1994, the M.S. and Ph.D. degrees from George Mason University, Fairfax, VA, in 1998 and 2000. Since August 2000, she has been on the faculty of Michigan Technological University, where she is currently a Professor. Dr. Tian’s general interests are in the areas of signal processing for wireless communications, estimation and detection theory. Current research focuses on cognitive radio networks and distributed wireless sensor networks. She served as Associate Editor for IEEE T RANSACTIONS ON W IRELESS C OMMUNICATIONS and IEEE T RANSACTIONS ON S IGNAL P ROCESSING. She received a CAREER award from the US National Science Foundation in 2003. Chunyan Feng received the B.S. degree in communications engineering, the M.S. and Ph.D. degrees in communication and information systems all from Beijing University of Posts and Telecommunications, Beijing, China. She is currently a Professor with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications. Her research interests are in the areas of broadband networks and wireless communication systems. Current research focuses on cognitive radio, key technology of B3G/4G systems, and green wireless communications.
