This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

Projected Barzilai-Borwein Methods Applied to Distributed Compressive Spectrum Sensing Le Thanh Tan∗ , Hyung Yun Kong∗ , Vo Nguyen Quoc Bao∗ ∗ School

of Electrical Engineering University of Ulsan, San 29 of MuGeo Dong, Nam-Gu, Ulsan, Korea 680-749 Email: {tanlh,hkong,baovnq}@mail.ulsan.ac.kr

Abstract—Cognitive radio allows unlicensed (cognitive) users to use licensed frequency bands by exploiting spectrum sensing techniques to detect whether or not the licensed (primary) users are present. In this paper, we present a compressed sensing applied to spectrum-occupancy detection in wide-band applications. The collected analog signals from each cognitive radio (CR) receiver at a fusion center are transformed to discrete-time signals by using analog-to-information converter (AIC) and then employed to calculate the autocorrelation. For signal reconstruction, we exploit a novel approach to solve the optimization problem consisting of minimizing both a quadratic (l2 ) error term and an l1 -regularization term. Specifically, we propose the Basic gradient projection (GP) and projected Barzilai-Borwein (PBB) algorithm to offer a better performance in terms of the mean squared error of the power spectrum density estimate and the detection probability of licensed signal occupancy.

I. INTRODUCTION The explosive growth of wireless communication has made the problem of spectrum utilization more critical. On the one hand, the increasing diversity and demand of high qualityof-service applications are the main reasons of the crowded nature of spectrum allocation. On the other hand, the spectrum scarcity is not the result of heavy usage of spectrum, but is due to the inefficiency of the static frequency allocation. For example, the typical spectrum utilization of around five percent or even less is reported in [1] . To solve this problem, the Federal Communications Commission (FCC) proposed the opening of licensed bands to unlicensed users and cognitive radio (CR) was born to improve the utilization of spectrum resource. In addition, the IEEE has established an IEEE 802.22 workgroup to build the standards of WRAN based on CR techniques [2]. Spectrum sensing for wide-band CR applications associates with considerable technical challenges. The radio terminals are required to employ a bank of tunable narrowband bandpass filters and search one narrow frequency band at a time. The large bandwidth operation requires an unfavorably large number of RF components. Furthermore, high speed processing units (DSPs or FPGAs) are needed to flexibly search over multiple frequency bands concurrently. In particular, input analog signals must be converted to discrete-time signals at Nyquist sampling rate or higher by using analog-to-digital converter (ADC). However, the extremely high sampling rates of wide-band ADCs are a main challenge with this model. Meanwhile, due to the timing requirements for rapid sensing,

only a limited number of measurements can be acquired from the received signal, which may not provide sufficient statistic when traditional linear signal reconstruction methods are employed. Compressed sensing (CS) [3], [4] builds on the surprising revelation that a signal having a sparse representation in one basis can be recovered from a small number of projections onto a second basis that is incoherent with the first. Especially, compressed sampling approach can get the sparse signal at the rates lower than the Nyquist sampling method; signal reconstruction which is a solution to a convex optimization problem called min-l1 with equality constraints makes use of the basic pursuit (BP) or some modified methods such as orthogonal matching pursuit (OMP), tree-based OMP (TOMP) [5], [6]. In [5], authors present a model of a single wide-band CR that uses CS based spectrum sensing schemes. This paper considers the problem of determining spectrum occupancy of a wide-band system. Our proposed CR system has a number of CRs and a centralized fusion center that collects data from individual CRs with different signal-to-noise ratios (SNRs) and decides these spectra to be available or not. Specifically, the same wide-band analog signal is transformed to discrete-time signal while preserving its salient information by using a so-called analog-to-information converter (AIC). For that kind of sparse input signal, an AIC can be acquired at sub-Nyquist rates (matching the information rate of the signal), and it also offers the capability of extracting the feature of data; therefore, the performance of the AIC system is comparable to that of ADC system in term of SNR or effective number of bits (ENOB) [7], [8], [9]. Then output signals of AIC are used to generate autocorrelation vectors that will be collected at the fusion. To obtain an estimate of the signal spectrum, we propose a novel version of distributed CS algorithm based on [10] by using a standard approach to minimize an objective function includes a quadratic (squared l2 ) error term combined with a sparseness-inducing (l1 ) regularization term. It can be easily observed that basic gradient projection (GP) can perform the high quality reconstruction [11]; however, the time-consuming property makes this approach not appropriate for spectrum sensing scheme which requires time-constraint. Thus, the novel approach, which exploits the projected Barzilai-Borwein (PBB) technique embedded in a continuation heuristic to recover the efficient performance [12], is introduced to improve both time-reduction and the

978-1-4244-5188-3/10/$26.00 ©2010 IEEE

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

performance of signal reconstruction. We also compare the performance of the distributed compressive spectrum sensing scheme with that of the scheme of [5] for a single CR to show the gains accrued from spatial diversity and exploiting the joint sparsity structure. We use (i) the mean squared error (MSE) between the reconstructed power spectrum density (PSD) estimate and the PSD based on Nyquist rate sampling, and (ii) the probability of detecting spectrum occupancy over the channels as performance measures. The rest of this paper is organized as follows. Section II contains the model of the wide-band analog signal compressed sensing. A compressive spectrum sensing scheme for single CR is presented in Section III. And an extention to collaborative compressed spectrum sensing for multiple CR is shown in Section IV. Section V demonstrates simulation results. Finally, concluding remarks are given in Section VI.

the received signal is wideband. It is natural to think about ways to avoid the high-speed ADC by applying CS to the analog signal directly. A related idea was first described in [8], where the analog signal was first demodulated with a pseudo-random chipping sequence p(t), then passed through an analog filter h(t), and the measurements were obtained in serial by sampling the filtered signal at sub-Nyquist rate. The serial sampling structure is appropriate for real-time processing. However, to achieve a satisfactory signal reconstruction quality, the order of the filter is usually higher than 10. In addition, because the measurements are obtained by sampling the output of the analog filter sequentially, they are no longer independent due to the convolution in the filter, which brings some redundancy in the measurements. Specifically, suppose that we have an analog signal x(t) which is K − sparse over some basis Ψ for t ∈ [0, T ] as in the following expression: x(t) =

II. SYSTEM MODEL OF WIDEBAND ANALOG SIGNAL COMPSESSED SENSING We consider the frequency range of interest to be comprised of maxI non-overlapping consecutive spectrum bands, a CR network consisting of J CRs and a centralized fusion center. Sensing is performed periodically at each CR and the results are sent to the fusion center, where a decision is made on whether or not there is a licensed signal present in each channel. B. Overview of compressed sensing According to Donoho [3], in CS theories, an N ×1 vector of discrete-time signal x = Ψs, where Ψ is the N × N sparsity basis matrix and s is the N × 1 vector with K  N non-zero (and large enough) entries si , can be used to reconstruct the signal from M measurements; especially, M depends on the reconstruction algorithm and is usually much less than N . This measurement can be done by projecting x on to an M × N basis matrix Φ that is incoherent with Ψ [13] (1)

The reconstruction is done by solving the following l1 -norm optimization problem as ˆ s= arg min s1 s.t. y = ΦΨs. s

si ψi (t),

(3)

i=1

A. Signal model

y = Φx = ΦΨs.

N 

(2)

Linear programming techniques, e.g., basis pursuit [14], or iterative greedy algorithms [15] can be used to solve (2). C. Compressed sensing of analog signals Because CS was proposed for discrete-time signal processing, we must use ADC sampling at Nyquist rate to discreterize the analog signal before applying the CS. After that, the compressed sensed data are sent to DSP blocks for further manipulation. While it is true that the data volume to be processed by DSP blocks is reduced due to the CS, a highspeed ADC sampling at Nyquist rate is still required when

where x is the N × 1 vector x = Ψs, Ψ is the N × N sparsity basis matrix Ψ = [ψ0 (t), ψ1 (t), . . . , ψN (t)] and s an N × 1 vector with K  N non-zero elements si . It has been shown that x can be recovered using M = KO(log N ) non-adaptive linear projection measurements on to an M × N basis matrix Φ that is incoherent with Ψ [13]. The received signal y can be viewed as the transmitted signal plus some additive noise y = Φx + n = ΦΨs + n.

(4)

There are several choices for the distribution of Φ such as Gaussian, Bernoulli. Reconstruction is achieved by solving the l1 -norm optimization problem as in (2). In this paper, the reconstruction problem, that has been highly interested in solving the convex unconstrained optimization problem, is a standard approach consisting in minimizing an objective function which includes a quadratic (squared l2 -norm) error term combined with a sparseness-inducing (l1 -norm) regularization term. So the problem can be given by 1 2 y − ΦΨs2 + τ s1 . (5) 2 Basic GP is able to solve a sequence of problems (5) efficiently for a sequence of values of τ . The gradient projection algorithms for solving a quadratic programming reformulation of a class of convex nonsmooth unconstrained optimization problems are significantly faster (in some cases by orders of magnitude), especially in large-scale settings. Instances of poor performance have been observed when the regularization parameter is small, but in such cases the gradient projection methods can be embedded in a simple continuation heuristic to recover their efficient practical performance. The new algorithms are easy to implement, work well across a large range of applications, and do not appear to require applicationspecific tuning. min s

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

u

x(t)

Analog filter h(t)

pc(t)

Fig. 1.

yk

k/M quantizer

Autocorrelation

Denote that φ∗i,j is the (i, j)-th element of ΦA , the M × N ¯ A has its (i, j)-th element given by matrix Φ 0 i = 1, j = 1, ..., N,   ¯A (14) Φ = i,j φM +2−i,j i = 1, j = 1, ..., N,

ry

AIC block

CS acquisition at individual CR sensing receiver.

III. COMPRESSIVE SPECTRUM SENSING AT SINGLE CR We begin by describing the CS acquisition and recovery scheme for a single CR (J = 1) case. Fig. 1 depicts the acquisition at a single CR sensing receiver. The analog baseband signal x(t) is sampled using an AIC. Following the circuit implementation of the AIC system in series of previous works [7], [8], [9], an AIC may be conceptually viewed as an ADC operating at the Nyquist rate, followed by compressive sampling. Denote the N × 1 stacked vector at the input of the ADC by T  k = 0, 1, 2 ..., xk = xkN xkN +1 ... xkN +N +1 (6) and the M × N compressive sampling matrix by ΦA . The output of the AIC denoted by the M × 1 vector T  k = 0, 1, 2 ..., yk = ykM ykM +1 ... ykM +M +1 (7) is given by (8) yk = ΦA xk . The respective N × N and M × M autocorrelation matrices of the compressed signal and the input signal vectors in (6) and (7) are related as follows:   H (9) R y = E yk y H k = ΦA Rx ΦA , where subscript H denotes the Hermitian. The elements of the matrices in (9) are given by: [Ry ]ij = ry (i − j) = ry∗ (j − i), [Rx ]ij = rx (i − j) = rx∗ (j − i). The respective 2N × 1 and 2M × 1 autocorrelation vectors corresponding to (6) and (7) can be expressed as follows:  T , rx = 0 rx (−N + 1) ... rx (0) ... rx (N − 1) (10)  T , ry = 0 ry (−M + 1) ... ry (0) ... ry (M − 1) (11) here the first zero values are artificially inserted. And these above vectors represent the first column and row of the respective autocorrelation matrices. To obtain the CS recovery like the formula (5), we must to make the relation between the autocorrelation vectors in (10) and (11). Using operations in matrix algebra, we can derive as

 Φ=

¯ A Φ1 Φ

¯ A Φ2 Φ

ΦA Φ3

ΦA Φ4

where hankel(c, r) is a hankel matrix (i.e., symmetric and constant across the anti-diagonals), note that c is the first column and r is the last row of this matrix. toeplitz(c, r) is a toeplitz matrix (i.e., symetric and constant across the diagonals), note that c is the first column, and r is the first row of this matrix. And 0a×b is the a × b zero matrix. We also know that using the wavelet-based edge detection in [16-17], the band boundaries (locations) can be recovered from 2N -1 local maxima of the wavelet modulus zs and the band number is determined by the number of local peaks; as an experiment when N  M in [5], zs can be recovered under the sparseness constraint, and therefore there is a linear transformation equality linking zs to the compressed measurement vector ry . And rx has a sparse representation in the edge spectrum domain [5], that is rx = Gzs ,

(15) −1

where zs is the discrete 2N × 1 vector, and G = (ΓFW) . The 2N × 2N matrices W and F represent respectively a wavelet-based smoothing and a Fourier transform. The 2N × 2N matrix Γ is a derivative operation given by ⎡ ⎤ 1 0 ··· 0 ⎢ ⎥ . ⎢ ⎥ ⎢ −1 . . · · · 0 ⎥ ⎥. Γ=⎢ ⎢ ⎥ ⎢ 0 ... ... ... ⎥ ⎣ ⎦ 0

···

−1

1

Combining (12) and (15), we can formulate the CS reconstruction of the edge spectrum as a convex unconstrained optimization problem: 1 2 ry − ΦGzs 2 + τ zs 1 (16) 2 To solve the above problem, we use the GP approach which is described in the Section IV for an individual CR case. The spectrum estimate can be evaluated as a cumulative T sum  zˆs (1) zˆs (2) · · · zˆs (2N ) . The of terms ˆ zs = min zs

discrete components of the PSD estimate are given by

ry = Φrx , note that

and the N × N matrices Φ1 , Φ2 , Φ3 , Φ4 are defined 

∗ ∗ , as Φ1 = hankel [0N ×1 ] , 0 φ1,1 ... φ1,N −1   

 φ∗1,1 ... φ∗1,N , φ∗1,N 01×(N −1) , Φ2 = hankel 

 Φ3 , = toeplitz [0N ×1 ] , 0 φ1, N ... φ1, 2   

 φ1,1 ... φ1,N , φ1,1 01×(N −1) Φ4 = toeplitz ,

(12) Sˆx (n) =

 .

(13)

n  k=1

zˆs (k)

(17)

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

k,1

1

•

y,1

•

J

The following subsections represent two algorithms to solve the above problem coresponding to two different ways of (k) (k) choosing αj and λj .

y,J

k,J

(k)

Step 3. Choose the second scalar λj ∈ [0, 1].   (k+1) (k) (k) (k) (k) Step 4. Set: pj = pj + λj wj − pj .

.

Fig. 2. Distributed compressive spectrum sensing scheme for multiple CRs.

B. Basic Gradient Projection:The GP - Basic Algorithm (k)

In this algorithm, we search from each iterate pj

IV. COLLABORATIVE COMPRESSED SPECTRUM SENSING Let xj (t) be the wide-band analog baseband signal received at the j-th CR sensing receiver. Each CR sensing receiver processes the received signal to obtain an 2M × 1 autocorrelation vector ry, j of the compressed signal, as in the CS acquisition step described in Section III, these vectors are then sent to the fusion center. The fusion center applies a GP algorithm to jointly reconstruct J received PSDs Sˆx,j , j = 1, . . . , J and then obtains an average PSD. The average PSD is then used to determine the spectrum occupancy. A. Overview of GP approach We now describe the GP algorithm used for reconstruction of J PSDs. We write A = ΦG in terms of its columns   (18) A = a1 a2 · · · a2N . At the j-th CR, we introduce vectors uj and vj and make the substitution zs,j = uj − vj , uj ≥ 0, vj ≥ 0, j = 1, · · · , J. + Here uj (i) = (zs,j (i)) = max {0, zs,j (i)} and vj (i) = + (−zs,j (i)) = max {0, −zs,j (i)} for all i = 1, . . . , 2N . Note + that (a) = max {0, a}. Therefore, we have zs,j  = 1T2N uj + T T 12N vj , where 12N = [1, 1, . . . , 1] is the vector consisting of 2N ones. The problem (16) can be modified as min u,v

s.t.

1 2

2

ry,j −A (uj − vj )2 +τ 1T2N uj +τ 1T2N vj uj ≥ 0 vj ≥ 0

(19)

Problem (19) can be written in more standard boundconstrained quadratic programming (BCQP) form as   min cT pj + 12 pTj Bpj ≡ F pj , p (20) s.t. pj ≥ 0     uj −b T where pj = , c = τ 14N + , b = AT ry,j vj b   AT A −AT A and B = . The next step is solving the −AT A AT A (k)

problem (20) by using a GP technique. From iterate pj (k+1) , we must follow the below steps, iterate pj • •

(k)

Step 1. Choose the scalar parameters αj > 0. + 

(k) (k) (k) (k) Step 2. Set: wj = pj − αj ∇F pj .

to

along the

(k) pj

, projecting onto the nonneganegative gradient −∇F tive orthant, and performing a backtracking line search until sufficient decrease is attained in F . We define the vector g(k) as       (k) ∇F pj (k) i , if pj,i > 0 or ∇F pj (k) i < 0 (k) gj,i = 0, otherwise. (21) where i = 1, . . . , 2N . The procedure of this algorithm is described as follows: 1) Input:   (0)

(0)

(0)

p2 · · · pJ . a) An initial p(0) = p1 b) A 2M × J data matrix R = ry,1 ry,2 · · · ry,J received from J CR sensing receivers. c) Choose parameters β ∈ (0, 1) and μ ∈ (0, 1/2). d) Set k = 0. 2) Output: A 2N × J reconstruction matrix Zs =  zs,1 zs,2 · · · zs,J , the average of J PSD estiˆ (J) mate vectors S x . 3) Procedure: a) Step 1. Compute α0,j as the following expression [12]: T

(k) (k) gj gj . (22) α0,j =

T (k) (k) gj Bgj Note that α0,j is solved from the expression:   (k) (k) (23) α0,j = arg min F pj − αj gj . αj

Then to guarantee that α0,j is not too small or too large, we replace α0,j by mid (αmin , α0,j , αmax ). Here mid (α1 , α2 , α3 ) is defined to be the middle value of three scalar values. b) Step 2. Backtracking line search: choose (k) αj to be the first number in the sequence α0,j , βα0,j , β 2 α0,j , . . . and satisfy the following inequality 

+ 



(k) (k) (k) (k) − ≤ F pj F pj −αj ∇F pj 

T +



(k) (k) (k) (k) (k) μ∇F pj pj − pj −αj ∇F pj (24)

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

and update the new set of values +

(k+1) (k) (k) (k) = pj − αj ∇F pj . pj c) Step 3. Termination test: (k) • Condition: We now convert pj =  

T

T T (k) (k) , vj to an approximate uj (k) zs,j,GP

(k) uj

T

(k) (k) (k) Bγj = 0, we choose λj = 1. Note that if γj (k)

3) Step 3. Update αj : Denote T

(k) (k) (k) Bγj . ξj = γj (k)

If ξj

(k) vj .

= − And scanning solution through the entire J CR to check the termination condition, i.e. the convergence gradient (CG) iteration is terminated when satisfying 2

2

ry,j − Azs,j 2 ≤ εD ry,j − Azs,j,GP 2 , (25) where εD is a small positive parameter. However, the termination iteration is also performed when the number of CG steps reaches to maxiterD. • We perform the convergence test and terminate (k+1) with approximation solution pj if it is satisfied the above conditions; otherwise we increase k to k + 1 and go back to step 1. = d) Step 4. Store the results: Store pj  T T T (uj ) , (vj ) , and calculate the reconstruction vector ˆ zs,j = uj − vj . The j-th PSD estimate n  ˆ x,j (n) = zˆs,j (k). And the average vector is S k=1

of J PSD estimate vectors is J 1 ˆ ˆ (J) Sx,j . S x = J j=1

(26)

C. Projected Barzilai-Borwein Reconstruction Algorithm The improvement of this algorithm is updating the step by the following formula [18] 

(27) γ (k) = −Hk−1 ∇F p(k) ,  (k)  . where Hk is an approximation to the Hessian of F p The procedure of this algorithm is similar to the basic GP algorithm except the following steps: (k) 1) Step 1. Compute step γj for the j-th CR as the following expression: +



(k) (k) (k) (k) (k) γj = pj − αj ∇F pj − pj . (28) (k)

(k)

2) Step 2. Line search: The scalar λj , (λj ∈  (k) (k) (k) [0, 1]) will be found to minimize F pj + λj γj (k+1)

and update the new set of values pj =

+

(k) (k) (k) pj − αj ∇F pj . Because F is quadratic, the (k) λj

can be evaluated by the line search parameter following closed-form expression: ⎧

T  ⎫

(k) (k) ⎪ ⎪ ⎨ ⎬ ∇F pj γj (k) , 1 λj = mid 0,

.  T ⎪ ⎪ (k) (k) ⎩ ⎭ γj Bγj

(29)

(k+1)

= 0, let αj (k+1)

αj

= αmax , otherwise ⎧ ⎫ 2 (k) ⎪ ⎪ ⎨ ⎬ γj 2 = mid αmin , , α . max (k) ⎪ ⎪ ξj ⎩ ⎭

D. Performances 1) MSE Performance: The normalized MSE of estimated PSD is computed by ⎧ ⎫ 2 (J) ⎪ ⎪ ˆ (J) ⎨ S ⎬ x − Sx (J) 2 , (30) MSE = E 2 ⎪ ⎪ ⎩ ⎭ Sx (J) 2

ˆ (J) where S and Sx (J) denote the average of the J PSD x estimate vectors based on our compressed sensing approach and the periodogram using the signals sampled at the Nyquist rate, respectively. 2) Detection performances: We evaluate the probability ˆ (J) of detection Pd based on the averaged PSD estimate S x . The detection analysis to follow, strictly speaking, holds only for samples collected at Nyquist rate. We however use this as a simple way to analyze the detection performance in the compressive sampling case as well. The decision of the presence of licensed transmission signals in the certain channel is made by an energy detector using the estimated frequency response over that channel, i.e., the test statistic is (J)

EI

IK 

=

ˆ (J) S x (i),

I = 1, 2, . . . , maxI,

(31)

i=(I−1)K+1

where I is the channel index, maxI is the number of channels, and K is the number of PSD samples of each channel. The PSD estimate of the j-th CR node can be evaluated as H  2 ˆx,j (i) = 1 |Xh,j (i)| , S H

(32)

h=1

whereXh,j (i) is the Fourier transform of the h-th block of the received signal xh,j (n), j representing the CR node index, n representing the time sample index, each block containing 2N time samples, and H denoting the number of blocks. Substituting (26) and (32) to (31), the test static can be obtained by (J)

EI

=

1 JH

IK 

J  H 

2

|Xh,j (i)| .

(33)

i=(I−1)K+1 j=1 h=1

The decision rule is chosen as H1 (J) > < H0

EI

μ,

I = 1, 2, . . . , maxI,

(34)

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

Parameters Elementary period T Number of carriers K Value of carrier number Kmin Value of carrier number Kmax Duration of symbol part TU Carrier spacing 1/TU Spacing between carriers Kmin and Kmax

1 Basic GP, 5CRs PBB, 5 CRs Basic GP, 1 CR PBB, 1 CR Ref. [5], 1 CR

0.9

Reconstruction error

0.8 0.7 0.6

2k mode 7/64µs 1,705 0 1,704 2,048× T 224µs 4,464 Hz 7.61 MHz

TABLE I T HE OFDM PARAMETERS FOR THE 2 K MODE .

0.5 0.4 0.3 0.2

V. SIMULATION RESULTS 0.1

0.2 0.3 Compression rate [M/N]

0.4

0.5

Fig. 3. Reconstruction error (MSE) for Basic GP and PBB approaches versus compression rate M/N for various number of collaborating CRs (SNRs of active channels varying from 8dB to 10dB).

1 0.9

Performance

0.8 0.7

Pd, 5 CRs

0.6

Pd, 1 CR P , 1CR, Ref. [5] d

0.5

Pfa, 5 CRs

0.4

Pfa, 1 CR Pfa, 1 CR, Ref. [5]

0.1 0.05 0

0.1

0.2 0.3 Compression rate [M/N]

0.4

0.5

Fig. 4. Probability of detection Pd and probability of a false alarm Pf a for PBB versus compression rate M/N for various number of collaborating CRs (SNRs of active channels varying from 8dB to 10dB).

where H0 , H1 represent the hypotheses of the absence and presence of primary signals, respectively, and μ is the decision (J)  threshold. Under H0 , EI / σn2 / (JH) ∼ χ22JKH has a central χ2 distribution with 2JKH degrees of freedom. The (J) probability of a false alarm Pf a can be obtained by (J) Pf a

  μ Γ JH, JH , =1− Γ (JH)

(35)

where Γ(., .) is the upper incomplete gamma function [19, Sec. (8.350)], Γ(.) is the gamma function [19, Sec. (13.10)]. Under (J) H1 , the probability of detection Pd is evaluated by (J)

Pd

=

Ia ! " 1  (J) Pr EI > μ, a

(36)

I=I1

where Ii , i = 1, . . . , a denote the indices of a active channels.

The model for simulation can be briefly described in this section. We consider at baseband, a wide frequency band of interest ranging from -38.05 to 38.05 MHz, containing maxI = 10 non-overlapping channels of equal bandwidth of 7.61 MHz. Our simulations will focus in the 2k mode of the DVB-T standard. This particular mode is intended for mobile reception of standard definition DTV. The structure of signal is followed an OFDM frame. Each frame has a duration of TF , and consists of 68 OFDM symbols. Four frames constitute one super-frame. Each symbol is constituted by a set of C = 1,705 carriers in the 2k mode and transmitted with a duration TS . A useful part with duration TU and a guard interval with a duration Δ (choosen to 0) compose TS . The over-sampling factor is 2. The occupancy ratio of the total 76.1 MHz band is 50%. The received signal is damaged by additive white Gaussian noise (AWGN) with a variance of σn2 = 1. The received SNRs on the a = 5 active channels are randomly varying from 8dB to 10dB. A Gaussian wavelet function is used for smoothing. For compressive sampling, 2N is 4096, the compressed rate M/N is varying from 5% to 50% and H = 160 is the number of blocks. The compressive sampling matrix ΦA has a Gaussian distributed function with zero mean and variance 1/M . The number of PSD samples of each channel is K = 25. We set αmin = 10−30 , αmax = 1030 for PBB algorithm, and use β = 0.5, μ = 0.1, and τ = 0.1 AT ry,i ∞ for both Basic and PBB algorithms. Fig. 3 illustrates MSE performance for Basic GP and PBB algorithms and compares with the performance result in [5]. In comparison with [5], our proposed approach in case of 1 CR slightly decreases the MSE performance because of the reduced mutual incoherent of Φ in (12), however, our approach can reduce the hardware cost due to AIC acquisition at the lower sampling rate. The results show that in comparison with Basic GP version, the PBB algorithm achieves the same performances while the Basic GP version takes a lot of time to get convergence [12]. So the following results are implemented by using the novel PBB algorithm. This figure also shows the performances of signal recovery quality in which MSE decreases when the value of compression rate M/N increases. However, as considering the effects of multiple CRs in spectrum sensing scheme, it is easily to observe that MSE also decreases as the number of CRs J increases; therefore, we can obtain the lower compression rate but not degrade the

This full text paper was peer reviewed at the direction of IEEE Communications Society subject matter experts for publication in the IEEE DySPAN 2010 proceedings

recovery performance by using more CRs in networks. For example, to receive MSE at 0.37, we must compress wideband signals at the rate 0.4 in the 1 CR case, while networks with 5 CRs can be implemented at the lower compression rate 0.3. For detection performance, Fig. 4 depicts the probability (J) of detection Pd with respect to both the compression ratio (J) M/N and the number of CRs J =1, 5, under a fixed Pf a of 0.01. This figure demonstrates that in order to obtain reliable performances, the joint collaboration and compression is necessary. Especially, collaboration among CRs can avoid the hardware cost of each CR by reducing the compression rate M/N while remaining the high detection performance. For instance, the probability of detection in case of 1 CR is ≈ 1 at the compression rate M/N over 0.2, while the collaboration among 5 CRs requires the compression rate M/N from the lower value 0.15. Especially, analyzing the results reveals the interesting conclusion, i.e., in Fig. 4, the detection performances under both our method and the approach in [5] over the examined range of compression rates are similar while in Fig.3, the MSE performances of these approaches have a bit differences. VI. CONCLUSION In this paper, we presented a distributed compressive spectrum sensing scheme for CR networks. To avoid the high speed ADC systems, the alternative converters called AICs are exploited to acquire the salient information of received signals at sub-Nyquist rates. Moreover, the GP approach is used for joint cooperation and compressive sensing. The major barrier of GP method, which takes a lot of time to reach the convergence, can be solved by modifying the backtracking line search for updating parameters. Among new fast CS techniques, PBB algorithm, which is used to update the step of the iterations in the recovery stage, demonstrates its outperformance, i.e., it not only achieves high quality of signal recovery but also increases the speed to quickly reach convergence. ACKNOWLEDGMENT This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (No. 2009-0073895) R EFERENCES [1] Proc. 1st IEEE Int’l. Symp. New Frontiers in Dynamic Spectrum Access Networks, Nov. 2005. [2] http://www.ieee802.org/22/. [3] Donoho, D.L., ”Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289-1306, 2006. [4] Candes, E.J., and Tao, T., ”Near-Optimal Signal Recovery From Random Projections: Universal Encoding Strategies?,” IEEE Transactions on Information Theory, vol. 52, no. 12, pp. 5406-5425, 2006. [5] T. Zhi, and G. B. Giannakis, ”Compressed Sensing for Wideband Cognitive Radios”, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. IV-1357-IV-1360, 2007. [6] M. F. Duarte, M. A. Davenport, M. B. Wakin et al., ”Sparse Signal Detection from Incoherent Projections”, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. III-III, 2006.

[7] S. Kirolos, T. Ragheb, J. Laska et al., ”Practical Issues in Implementing Analog-to-Information Converters”,in The 6th International Workshop on System-on-Chip for Real-Time Applications, pp. 141-146, 2006. [8] J. N. Laska, S. Kirolos, M. F. Duarte et al., ”Theory and Implementation of an Analog-to-Information Converter using Random Demodulation”, in IEEE International Symposium on Circuits and Systems (ISCAS), pp. 1959-1962, 2007. [9] J. L. S. Kirolos, M. Wakin, M. Duarte, D. Baron, T. Ragheb, Y. Massoud, and R. Baraniuk, ”Analog-to-information conversion via random demodulation,” in Proc. of the IEEE Dallas Circuits and Systems Workshop (DCAS), 2006. [10] M. F. Duarte, S. Sarvotham, D. Baron et al., ”Distributed Compressed Sensing of Jointly Sparse Signals,” in Conference Record of the ThirtyNinth Asilomar Conference on Signals, Systems and Computers, pp. 1537-1541, 2005. [11] J. J. More, ”Gradient projection techniques for large-scale optimization problems,” in Proceedings of the 28th IEEE Conference on Decision and Control, pp. 378-381 vol.371, 1989. [12] Y.-H. D. a. R. Fletcher, ”Projected Barzilai-Borwein methods for largescale box-constrained quadratic programming,” Numerische Mathematik, vol. 100, pp. 21-47, March, 2005. [13] E. C. a. J. Romberg, ”Sparsity and Incoherence in Compressive Sampling,” Inverse Problems, 23(3), pp. 969-985, June 2007. [14] D. L. D. S. S. Chen, and M. A. Saunders, ”Atomic decomposition by basis pursuit,” SIAMJ Sci. Comput., vol. 20, no. 1, pp. 33-61, 1999. [15] C. L. a. M. Do, ”Signal reconstruction using sparse tree representations,” SPIE Wavelets XI, vol. 5914, pp. 59140W.1-59140W.11, Sept 2005. [16] Z. Tian, and G. B. Giannakis, ”A Wavelet Approach to Wideband Spectrum Sensing for Cognitive Radios,” in 1st International Conference on Cognitive Radio Oriented Wireless Networks and Communications, pp. 1-5, 2006. [17] S. Mallat, and W. L. Hwang, ”Singularity detection and processing with wavelets,” IEEE Transactions on Information Theory, vol. 38, no. 2, pp. 617-643, 1992. [18] J. Barzilai and J. M. Borwein, ”Two-point step size gradient methods,” IMA J. Numer. Anal., vol.8, no.1, pp. 141-148, January 1, 1988. [19] I. S. Gradshtein, I. M. Ryzhik, A. Jeffrey and D. Zwillinger, ”Table of integrals, series and products,” 7th ed., Oxford: Academic, 2007.