Phased Array-Based Sub-Nyquist Sampling for Joint Wideband Spectrum Sensing and Direction-of-Arrival Estimation
aa r X i v : . [ ee ss . SP ] O c t Phased Array-Based Sub-Nyquist Sampling forJoint Wideband Spectrum Sensing andDirection-of-Arrival Estimation
Feiyu Wang, Jun Fang, Huiping Duan, and Hongbin Li,
Senior Member, IEEE
Abstract — In this paper, we study the problem of joint wide-band spectrum sensing and direction-of-arrival (DoA) estimationin a sub-Nyquist sampling framework. Specifically, consideringa scenario where a few uncorrelated narrowband signals spreadover a wide (say, several GHz) frequency band, our objective is toestimate the carrier frequencies and the DoAs associated with thenarrowband sources, as well as reconstruct the power spectra ofthese narrowband signals. To overcome the sampling rate bottle-neck for wideband spectrum sensing, we propose a new phased-array based sub-Nyquist sampling architecture with variable timedelays, where a uniform linear array (ULA) is employed and thereceived signal at each antenna is delayed by a variable amountof time and then sampled by a synchronized low-rate analog-digital converter (ADC). Based on the collected sub-Nyquistsamples, we calculate a set of cross-correlation matrices withdifferent time lags, and develop a CANDECOMP/PARAFAC (CP)decomposition-based method for joint DoA, carrier frequencyand power spectrum recovery. Perfect recovery conditions forthe associated parameters and the power spectrum are analyzed.Our analysis reveals that our proposed method does not requireto place any sparse constraint on the wideband spectrum, onlyneeds the sampling rate to be greater than the bandwidth of thenarrowband source signal with the largest bandwidth among allsources. Simulation results show that our proposed method canachieve an estimation accuracy close to the associated Cram´er-Rao bounds (CRBs) using only a small number of data samples.
Index Terms — Joint wideband spectrum sensing and direction-of-arrival (DoA) estimation; compressed sensing; CANDE-COMP/PARAFAC (CP) decomposition.
I. I
NTRODUCTION
Wideband spectrum sensing, which aims to identify thefrequency locations of a few narrowband transmissions thatspread over a wide frequency band, has been of a growinginterest in signal processing and cognitive radio communi-cations [1], [2]. To perform wideband spectrum sensing, aconventional receiver requires to sample the received signalat the Nyquist rate, which may be infeasible if the spectrumunder monitoring is very wide, say, reaches several GHz. Also,
Feiyu Wang, and Jun Fang are with the National Key Laboratory of Scienceand Technology on Communications, University of Electronic Science andTechnology of China, Chengdu 611731, China, Email: [email protected] Duan is with the School of Electronic Engineering, University ofElectronic Science and Technology of China, Chengdu 611731, China, Email:[email protected] Li is with the Department of Electrical and Computer Engineering,Stevens Institute of Technology, Hoboken, NJ 07030, USA, E-mail: [email protected] work was supported in part by the National Science Foundation ofChina under Grant 61522104, and the National Science Foundation underGrant ECCS-1408182 and Grant ECCS-1609393. a high sampling rate results in a large amount of data whichplace a heavy burden on subsequent storage and processing.To alleviate the sampling rate requirement, a variety of sub-Nyquist sampling schemes, e.g. [3]–[6], were developed. Therationale behind such schemes is to exploit the inherent spar-sity in the frequency domain and formulate wideband spectrumsensing as a sparse signal recovery problem which, accordingto the compressed sensing theory [7], [8], can perfectly recoverthe signal of the entire frequency band based on compressedmeasurements or sub-Nyquist samples. Furthermore, in [9]–[11], it was shown that it is even possible to perfectlyreconstruct the power spectrum without placing any sparseconstraint on the wideband spectrum under monitoring.In some applications such as electronic warfare, one neednot only conduct wideband spectrum sensing, but also iden-tify the carrier frequencies and directions-of-arrival (DoAs)associated with the narrowband signals that live within thewide frequency band [12]. Besides, in massive MIMO ormillimeter wave systems where signals are transmitted viabeamforming techniques, the DoA information would allowa cognitive radio to more efficiently exploit the vacant bands[13]. In [14], [15], ESPRIT-based methods were proposed forjoint carrier frequency and DoA estimation. These methods,however, require the signal to be sampled at the Nyquist rate.Recently, with the advent of compressed sensing theories,the sparsity inherent in the spectral and spatial domains wasutilized to devise sub-Nyquist sampling-based algorithms forjoint wideband spectrum sensing and DoA estimation. Specifi-cally, in [16], a compressed sensing method was developed in aphased array framework, where a multicoset sampling schemeis executed at each antenna to collect non-uniform samples.In practice, the multicoset sampling may be implementedusing multiple channels, with each channel delayed by adifferent time offset and then sampled by a low-rate analog-digital converter (ADC). Since the multicoset sampling hasto be performed at each antenna, the scheme [16] involvesa high hardware complexity. In [17], [18], a simplified sub-Nyquist receiver architecture was proposed, in which eachantenna output is connected with only two channels, i.e. adirect path and a delayed path. An ESPRIT-based algorithmwas then developed for joint DoA, carrier frequency, andsignal reconstruction. In addition to the above time delay-based sub-Nyquist receiver architectures, an alternative sub-Nyquist sampling approach, referred to as phased array-basedmodulated wideband converter (MWC), was proposed in [13],[19] for carrier and DoA estimation. The receiver utilizes an
L-shaped array, and all sensors have the same sampling patternimplementing a single channel of the MWC. Perfect recoveryconditions were analyzed, and reconstruction algorithms basedon compressed sensing techniques were developed in [19].In this paper, we propose a new sub-Nyquist receiverarchitecture, referred to as the phased-array based sub-Nyquistsampling architecture with variable time delays, for joint wide-band spectrum sensing and DoA estimation. Similar to [16]–[18], the proposed receiver architecture employs a uniformlinear array. The received signal at each antenna is delayedby a pre-specified time shift and then sampled at a sub-Nyquist sampling rate. Compared with existing sub-Nyquistreceiver architectures, our proposed sub-Nyquist scheme issimpler and easier to implement: it requires only one ADCfor each antenna output, thus leading to a lower hardwarecomplexity. Meanwhile, in our proposed architecture, the timedelays for different antennas can be arbitrary as long as theysatisfy a mild condition, which relaxes the requirement on theaccuracy of time delay lines. From the collected sub-Nyquistsamples, we calculate a set of cross-correlation matrices withdifferent time lags, based on which a third-order tensor thatadmits a CANDECOMP/PARAFAC (CP) decomposition canbe constructed. We show that the DoAs and the carrierfrequencies, along with the power spectra associated withthe sources, can be recovered from the factor matrices. Theperfect recovery condition is analyzed. Our analysis showsthat, to perfectly recover the power spectrum of the widefrequency band and the associated parameters, we only needthe sampling rate to be greater than the bandwidth of thenarrowband source signal with the largest bandwidth amongall sources. In addition, our proposed method does not needto impose any sparse constraint on the wideband spectrum.We also derive the Cram´er-Rao bound (CRB) results for ourestimation problem. Simulation results show that our proposedmethod, with only a small number of data samples, can achievean estimation accuracy close to the associated CRBs.We notice that a CP decomposition-based approach wasproposed in [13] for joint DoA and carrier frequency es-timation. Different from our work, the construction of thetensor in [13] has to rely on an L-shaped array and exploitsthe cross-correlations between the two mutually perpendicularsub-arrays. In addition, the PARAFAC analysis in [13] canonly help extract the DoA and carrier frequency information,while in our proposed method, the DoA, carrier frequency,and power spectrum associated with each source can besimultaneously recovered from the CP decomposition.The rest of the paper is organized as follows. In SectionII, we provide notations and basics on the CP decomposition.The signal model and related assumptions are discussed inSection III. In Section IV, we propose a new phase-arraybased sub-Nyquist receiver architecture. A CP decomposition-based method for joint wideband spectrum sensing and DoAestimation is developed in Section V. The uniqueness of theCP decomposition is discussed in Section VI, and the CRBanalysis is conducted in Section VII. Simulation results areprovided in Section VIII, followed by concluding remarks inSection IX. = + + … (1)1 a (2)1 a (3)1 a (1) a R (2) a R (3) a R R l · l · χ I I I Fig. 1. Schematic of CP decomposition.
II. P
RELIMINARIES
To make the paper self-contained, we provide a brief reviewon tensors and the CP decomposition. More details regardingthe notations and basics on tensors can be found in [20].Simply speaking, a tensor is a generalization of a matrixto higher-order dimensions, also known as ways or modes.Vectors and matrices can be viewed as special cases of tensorswith one and two modes, respectively. Throughout this paper,we use symbols ⊗ , ◦ , and ⊙ to denote the Kronecker, outer,and Khatri-Rao product, respectively.Let X ∈ C I × I ×···× I N denote an N th-order tensor with its ( i , . . . , i N ) th entry denoted by X i ··· i N . Here the order N ofa tensor is the number of dimensions. Fibers are higher-orderanalogues of matrix rows and columns. The mode- n fibers of X are I n -dimensional vectors obtained by fixing every indexbut i n . Slices are two-dimensional sections of a tensor, definedby fixing all but two indices. Unfolding or matricization isan operation that turns a tensor into a matrix. The mode- n unfolding of a tensor X , denoted as X ( n ) , arranges the mode- n fibers to be the columns of the resulting matrix. The CPdecomposition decomposes a tensor into a sum of rank-onecomponent tensors (see Fig. 1), i.e. X = R X r =1 λ r a (1) r ◦ a (2) r ◦ · · · ◦ a ( N ) r (1)where a ( n ) r ∈ C I n , the minimum achievable R is referred to asthe rank of the tensor, and A ( n ) , [ a ( n )1 . . . a ( n ) R ] ∈ C I n × R denotes the factor matrix along the n -th mode. Elementwise,we have X i i ··· i N = R X r =1 λ r a (1) i r a (2) i r · · · a ( N ) i N r (2)The mode- n unfolding of X can be expressed as X ( n ) = A ( n ) Λ (cid:16) A ( N ) ⊙ · · · A ( n +1) ⊙ A ( n − ⊙ · · · A (1) (cid:17) T (3)where Λ , diag ( λ , . . . , λ R ) .III. S IGNAL M ODEL
Consider a scenario in which K uncorrelated, wide-sensestationary, and far-field narrowband signals spreading over awide frequency band impinge on a wideband uniform lineararray (ULA) with N receiver antennas, where we assume N >K . Let s ( t ) denote the combination of the K narrowbandsignals in the time domain. s ( t ) can be expressed as s ( t ) = K X k =1 s k ( t ) e jω k t (4) Fig. 2. Proposed Phased-Array based Sub-nyquist Sampling Architecturewith variable Time delays (PASSAT). where s k ( t ) and ω k ∈ R + denote the complex basebandsignal and the carrier frequency (in radians per second) ofthe k th source signal, respectively. Each source signal s k ( t ) is associated with an unknown azimuth DoA θ k ∈ [0 , π ) . Wehave the following assumptions regarding the source signals:A1 The K source signals { s k ( t ) } are assumed to be mutuallyuncorrelated, wide-sense stationary, and bandlimited to [ − B/ , B/ , i.e. B k ≤ B, ∀ k , where B k denotes thebandwidth of the k th source signal.A2 Sources either have distinct carrier frequencies { ω k } ordistinct DoAs { θ k } , i.e. for any two source signals, wehave ( θ i , ω i ) = ( θ j , ω j ) , ∀ i = j .A3 The multi-band signal s ( t ) is bandlimited to F =[0 , f nyq ] , and we assume f nyq ≫ B .Assumption A2 is assumed to make signals distinguishedfrom one another. Note that this assumption is less restrictivethan the one made in other works, e.g. [17], [19], which, inorder to remove the source ambiguity, require the quantity ω k cos( θ k ) to be mutually different for different signals, i.e. ω i cos( θ i ) = ω j cos( θ j ) ∀ i = j (5)After collecting the received signal at the array, our ob-jective is to jointly estimate the DoAs { θ k } , the carrierfrequencies { ω k } , as well as the power spectra associatedwith the K source signals. To accomplish this task, we, inthe following, propose a new phased-array based sub-Nyquistreceiver architecture.IV. P ROPOSED S UB -N YQUIST R ECEIVER A RCHITECTURE
A. Proposed Receiver Architecture
In our receiver architecture, the received signal at eachantenna is delayed by a pre-specified factor ∆ n and thensampled by a synchronous ADC with a sampling rate of f s =1 /T s , where f s ≪ f nyq . We have the following assumptionsregarding the delay factors and the sampling rate:A4 The time delay factors { ∆ n } can take arbitrary values aslong as the following condition holds valid (∆ n +2 − n +1 + ∆ n ) f nyq < (6)for some n ∈ { , . . . , N − } .A5 The sampling rate f s is no less than the bandwidthof the narrowband source signal which has the largestbandwidth among all sources, i.e. f s ≥ B .As will be shown later in our paper, Assumption A4 isessential to identify the unknown carrier frequencies. Also, in practice, the time delay factors { ∆ n } can be chosen to be ofthe same order of magnitude as the Nyquist sampling intervalsuch that the narrowband approximation in (7) holds valid. Theproposed receiver architecture, termed as the Phased-Arraybased Sub-Nyquist Sampling Architecture with variable Timedelays (PASSAT), is illustrated in Fig. 2. The analog signalobserved by the n th antenna can be expressed as x n ( t ) = K X k =1 s k ( t − ( n − τ k − ∆ n ) × e jω k ( t − ( n − τ k − ∆ n ) + w n ( t ) ≈ K X k =1 s k ( t ) e jω k ( t − ( n − τ k − ∆ n ) + w n ( t ) (7)where the approximation is due to the narrowband assumption, w n ( t ) represents the additive white Gaussian noise with zeromean and variance σ , and τ k denotes the delay between twoadjacent sensors for a plane wave arriving in the direction θ k and is given by τ k = d cos θ k C (8)Here d denotes the distance between two adjacent antennasand we assumeA6 The distance between two adjacent antennas d satisfies d < C/f nyq , where C is the speed of light.We will show later that this assumption is essential for therecovery of the DoAs.In practice, only the real part of x n ( t ) is observed and sam-pled. Nevertheless, the corresponding imaginary part ℑ [ x n ( t )] can be retrieved from the real part ℜ [ x n ( t )] by passingthe signal through a finite impulse response (FIR) Hilberttransformer. The complex analytic signal can also be roughlyapproximated by computing the discrete Fourier transform(DFT) of the output of each antenna and throwing away thenegative frequency portion of the spectrum [12]. B. Relation to and Distinction from Existing Architectures
We notice that a time delay-based sub-Nyquist architecturewas also introduced in [16]–[18]. Nevertheless, there are twokey distinctions between our architecture and theirs. Firstly,our architecture has a simpler structure with only N delaychannels, whereas the architecture proposed in [17] (see Fig.3(a)) requires N channels in total, in which each antennaoutput passes through two channels, namely, a direct path anda delayed path. As a consequence, the number of requiredADCs for the architecture [17] is twice the number of ADCsfor our architecture. In [18], a modified architecture was pro-posed based on [17]. It, however, still requires N channels,with an N -channel delay network added to the first antenna.Secondly, for our proposed architecture, the time delays cantake arbitrary values as long as the mild condition (6) issatisfied. In contrast, for other architectures, e.g. [16]–[18],a precise time control is required such that the time delaysacross different channels are strictly identical [17], [18], or thetime delays must be integer multiples of the Nyquist samplinginterval [16]. Due to the inaccuracy caused by the time shift (a)(b)Fig. 3. Existing sub-Nyquist receiver architectures. (a) Sub-Nyquist samplingarchitecture proposed in [17]. (b) Sub-Nyquist sampling architecture proposedin [19]. elements, maintaining accurate time delays on the order of theNyquist sampling interval is difficult. The inaccuracy in thesedelays will impair the recovery performance. Our architectureis free from this issue because it allows more flexible timedelays and we can use the actual time delays measured inpractice for our proposed recovery algorithm.In [13], [19], a phased-array MWC-based sub-Nyquistsampling architecture (see Fig. 3(b)) was proposed for jointwideband spectrum sensing and DoA estimation, in which anL-shaped array composed of N + 1 sensors is adopted, andthe output of each sensor is multiplied by a same periodicpseudo-random sequence, low-pass filtered and then sampledat a low rate. Compared to the phased-array MWC-basedsub-Nyquist sampling architecture, our proposed delay-basedscheme is much simpler to implement. In [21], it is argued thatthe delay-based architectures suffer two major disadvantageswhich include the need for high-precision delay lines as well asspecialized ADCs with high analog bandwidth. Nevertheless,as discussed above, our proposed architecture, different fromother delay-based schemes [16]–[18], has a relaxed require-ment on the precision of delay lines. Regarding the latter issue,it is known that there is an inherent bandwidth limitation forpractical ADCs, termed analog (full-power) bandwidth, whichdetermines the highest frequency that can be handled by thedevice. In spite of that, we notice that the inherent bandwidthof some inexpensive, low-end commercial ADCs such asADC12DC105 can reach up to 1GHz, while some high-endADCs with affordable prices, such as ADC12D500, have aninherent bandwidth up to 2.7GHz, which may accommodatemost wideband spectrum sensing applications. V. P ROPOSED
CP D
ECOMPOSITION -B ASED M ETHOD
Let δ ( · ) denote the indicator function defined as δ ( x ) = (cid:26) , x = 00 , x = 0 . (9)We first calculate the cross-correlation between two sensoroutputs x m ( t ) and x n ( t ) . Recalling Assumption A1, wehave R xmn ( t , t ) = E [ x m ( t ) x ∗ n ( t )]= K X k =1 R sk ( t , t ) a mk a ∗ nk + R wmn ( t , t ) (10)where R sk ( t , t ) , E (cid:2) s k ( t ) e jω k t s ∗ k ( t ) e − jω k t (cid:3) (11)denotes the autocorrelation of the k -th modulated sourcesignal, R wmn ( t , t ) , E [ w m ( t ) w ∗ n ( t )] = σ δ ( m − n ) δ ( t − t ) (12)represents the autocorrelation of the additive noise and a nk , e − j (( n − τ k ω k +∆ n ω k ) (13)Since the source signals are wide-sense stationary, theautocorrelation R sk ( t , t ) depends only on the time difference t − t . As a result, the cross-correlation of the sensor outputs R xmn ( t , t ) depends on the time difference t − t as well.Let T s denote the sampling interval of the ADCs. The timedifference has to be an integer multiple of the samplinginterval, i.e. t − t = lT s for l = − L, . . . , L . For notationalconvenience, we define r xm,n ( l ) , R xmn ( t + lT s , t ) r sk ( l ) , R sk ( t + lT s , t ) r wm,n ( l ) , R wmn ( t + lT s , t ) We can therefore express (10) as a discrete-time form: r xm,n ( l ) = K X k =1 r sk ( l ) a mk a ∗ nk + r wm,n ( l ) (14)for l = − L, . . . , L and m, n = 1 , . . . , N .Our objective is to recover the DoAs { θ k } , the carrierfrequencies { ω k } , as well as the power spectra associatedwith the K source signals based on the second-order statistics { r xm,n ( l ) } . For each time lag l , we can construct a correlationmatrix R x ( l ) with its ( m, n ) th entry given by r xm,n ( l ) . Also,it can be easily verified that R x ( l ) = K X k =1 r sk ( l ) a k a Hk + R w ( l ) (15)where R w ( l ) denotes the cross-correlation matrix of theadditive noise with its ( m, n ) th entry given by r wm,n ( l ) , and a k , [ a k a k . . . a Nk ] T (16)Since a set of cross-correlation matrices { R x ( l ) } Ll = − L areavailable, we can naturally express this set of correlation matrices by a third-order tensor R x ∈ C (2 L − × N × N whosethree modes respectively stand for the time lag l and theantenna indices, and its ( l, m, n ) -th entry given by r xm,n ( l ) .Notice from (15) that each slice of the tensor R x , R x ( l ) , isa weighted sum of a common set of rank-one outer products.The tensor R x thus admits a CP decomposition which de-composes a tensor into a sum of rank-one component tensors,i.e. R x = K X k =1 r k ◦ a k ◦ a ∗ k + R w (17)where ◦ denotes the outer product, R w ∈ C (2 L − × N × N withits ( l, m, n ) -th entry given by r wm,n ( l ) , and r k , [ r sk ( − L ) . . . r sk ( L )] T (18)Define R , [ r . . . r K ] (19) A , [ a . . . a K ] (20)The three matrices { R , A , A ∗ } are referred to as factormatrices associated with the noiseless version of R x . We seethat the information about the parameters { θ k , ω k } as well asthe power spectra can be extracted from the factor matrices.Motivated by this observation, we propose a two-stage methodwhich consists of a CP decomposition stage whose objectiveis to estimate the factor matrices and a parameter estimationstage whose objective is to jointly recover the DoAs, carrierfrequencies, and the power spectra of sources based on theestimated factor matrices. A. CP Decomposition
We first consider the scenario where the number of sources, K , is known or estimated a priori . Clearly, the CP de-composition can be accomplished by solving the followingoptimization problem min ˆ R , ˆ A (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) R x − K X k =1 ˆ r k ◦ ˆ a k ◦ ˆ a ∗ k (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F (21)where ˆ R = [ ˆ r . . . ˆ r K ] , ˆ A = [ ˆ a . . . ˆ a K ] , and k · k F denotes the Frobenius norm. On the other hand, note that theCP decomposition is unique under a mild condition. Thereforewe can use a new variable ˆ b k to replace ˆ a ∗ k , which leads to min ˆ R , ˆ A , ˆ B (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) R x − K X k =1 ˆ r k ◦ ˆ a k ◦ ˆ b k (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F (22)where ˆ B , [ ˆ b . . . ˆ b K ] . The above optimization can beefficiently solved through an alternating least squares (ALS)procedure which alternatively updates one of the factor matri-ces to minimize the data fitting error while keeping the other two factor matrices fixed: b R ( t ) = arg min R (cid:13)(cid:13)(cid:13)(cid:13) ( R x (1) ) T − ( b B ( t − ⊙ b A ( t − ) R T (cid:13)(cid:13)(cid:13)(cid:13) F (23) b A ( t ) = arg min A (cid:13)(cid:13)(cid:13)(cid:13) ( R x (2) ) T − ( b B ( t − ⊙ b R ( t ) ) A T (cid:13)(cid:13)(cid:13)(cid:13) F (24) b B ( t ) = arg min B (cid:13)(cid:13)(cid:13)(cid:13) ( R x (3) ) T − ( b A ( t ) ⊙ b R ( t ) ) B T (cid:13)(cid:13)(cid:13)(cid:13) F (25)where R x ( n ) denotes the mode- n unfolding of R x .If the knowledge of the number of sources, K , is unavail-able, more sophisticated CP decomposition techniques (e.g.[22]–[24]) can be employed to jointly estimate the modelorder and the factor matrices. The basic idea is to use lowrank-promoting priors or functions to automatically determinethe CP rank of the tensor. In [22], when the CP rank, K , isunknown, the following optimization was employed for CPdecomposition min ˆ R , ˆ A , ˆ B k R x − X k F + µ (cid:16) tr ( ˆ R ˆ R H ) + tr ( ˆ A ˆ A H ) + tr ( ˆ B ˆ B H ) (cid:17) s.t. X = ˆ K X k =1 ˆ r k ◦ ˆ a k ◦ ˆ b k (26)where ˆ K ≫ K denotes an overestimated CP rank, µ is aregularization parameter to control the tradeoff between low-rankness and the data fitting error, ˆ R = [ ˆ r . . . ˆ r ˆ K ] , ˆ A =[ ˆ a . . . ˆ a ˆ K ] , and ˆ B = [ ˆ b . . . ˆ b ˆ K ] . The above optimization(26) can still be solved by an ALS procedure as follows ˆ R ( t ) = arg min ˆ R (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:20) ( R x (1) ) T (cid:21) − " ˆ B ( t − ⊙ ˆ A ( t − √ µ I ˆ R T (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F ˆ A ( t ) = arg min ˆ A (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:20) ( R x (2) ) T (cid:21) − " ˆ B ( t − ⊙ ˆ R ( t ) √ µ I ˆ A T (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F ˆ B ( t ) = arg min ˆ B (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:20) ( R x (3) ) T (cid:21) − " ˆ A ( t ) ⊙ ˆ R ( t ) √ µ I ˆ B T (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F The true CP rank of the tensor, K , can be estimated byremoving those negligible rank-one tensor components afterconvergence. B. Joint DoA, Carrier Frequency and Power Spectrum Esti-mation
We discuss how to jointly recover the DoAs, carrier fre-quencies, and power spectra of sources based on the estimatedfactor matrices. As shown in the next subsection, the CPdecomposition is unique up to scaling and permutation ambi-guities under a mild condition. More precisely, the estimatedfactor matrices and the true factor matrices are related as ˆ R = R Λ Π + E (27) ˆ A = A Λ Π + E (28) ˆ B = A ∗ Λ Π + E (29)where { Λ , Λ , Λ } are unknown nonsingular diagonal ma-trices which satisfy Λ Λ Λ = I ; Π is an unknown per-mutation matrix; and E , E , and E denote the estimation errors associated with the three estimated factor matrices,respectively. The permutation matrix Π can be ignored as itis common to all three factor matrices. Also, since we haveprior knowledge that columns of A / √ N have unit norm, theamplitude ambiguity can be estimated and removed, in whichcase we can write ˆ R = R ˜Λ + ˜ E (30) ˆ A = A ˜Λ + ˜ E (31) ˆ B = A ∗ ˜Λ + ˜ E (32)where ˜Λ , ˜Λ , ˜Λ are unknown nonsingular diagonal matriceswith their diagonal elements lying on the unit circle.Notice that the k th column of A is characterized by theDoA and carrier frequency associated with the k th source.We now discuss how to estimate { ω k } and { τ k } from theestimated factor matrix ˆ A . Note that ˆ B is also an estimate of A . Therefore either ˆ A or ˆ B can be used to estimate { ω k } and { τ k } . Let ˆ a k denote the k -th column of ˆ A , and write ˜Λ = diag { e − jϕ , . . . , e − jϕ K } (33)where { ϕ k } ∈ [0 , π ) are unknown parameters. To simplifyour exposition, we ignore the estimation errors ˜ E , ˜ E , and ˜ E .Write z = re jϕ , and define arg( z ) , mod ( ϕ, π ) arg( z ) ∈ [0 , π ) (34)where mod ( a, b ) is a modulo operator which returns theremainder of the Euclidean division of a by b . Recalling (13),we have η nk , mod ( − arg(ˆ a nk ) , π )= mod (( n − τ k ω k + ∆ n ω k + ϕ k , π ) (35)where ˆ a nk denotes the n th entry of ˆ a k . Let η k , [ η k . . . η Nk ] T and let D p denote a difference matrix defined as D p , − . . . − . . . ... ... . . . . . . ... . . . − ∈ R ( p − × p To recover ω k , we conduct a two-stage difference operationas follows β (1) k = mod ( D N η k , π ) (36) β (2) k = mod ( D N − β (1) k , π ) (37)It can be easily verified that entries of β (1) k and β (2) k arerespectively given as β (1) nk = mod ( τ k ω k + (∆ n +1 − ∆ n ) ω k , π ) , n = 1 , . . . , N − (38) β (2) nk = mod ((∆ n +2 − n +1 + ∆ n ) ω k , π ) , n = 1 , . . . , N − (39) From (39), we can see that the information about thecarrier frequency ω k is extracted after performing the two-stage difference operation. By properly devising the time delayfactors { ∆ n } , we can ensure that for some n ∈ { , . . . , N } ,the condition (6) holds valid, i.e. (∆ n +2 − n +1 + ∆ n ) f nyq < (40)The above condition implies (∆ n +2 − n +1 + ∆ n ) ω max < π (41)where ω max , max { ω , . . . , ω K } . Therefore ω k can simplybe estimated as ˆ ω k = β (2) n ,k ∆ n +2 − n +1 + ∆ n (42)In fact, for a careful selection of time delay factors { ∆ n } , thecondition (6) (i.e. (40)) may be satisfied for different choicesof n . As a result, we can obtain multiple estimates of ˆ ω k . Toimprove the estimation performance, a final estimate of ˆ ω k can be chosen as the average of these multiple estimates.Under Assumption A6, that is, d < C/f nyq , we have τ k ω max < π . Thus, substituting the estimated ˆ ω k back into(38), τ k can be obtained as ˆ τ k = mod (cid:16) β (1) nk − (∆ n +1 − ∆ n )ˆ ω k , π (cid:17) ˆ ω k (43)Note that for each β (1) nk , n = 1 , . . . , N − , we can obtainan estimate of τ k . Therefore multiple estimates of τ k can becollected. Again, an average operation can be conducted toyield a final estimate of τ k . Based on ˆ τ k , an estimate of theassociated DoA θ k can be readily obtained from (8).We now discuss how to recover the power spectra of thesources { s k ( t ) } . Let ˜ r sk ( τ ) , R sk ( t + τ, t ) , where τ ∈ R canbe any real value. The power spectrum of the k th source canthus be expressed as the Fourier transform of ˜ r sk ( τ ) , i.e. ˜ S k ( ω ) = Z + ∞−∞ ˜ r sk ( τ ) e − jωτ d τ (44)Let S k ( ω ) denote the discrete-time Fourier transform (DTFT)of the autocorrelation sequence { r sk ( l ) } + ∞ l = −∞ , i.e. S k ( ω ) = ∞ X l = −∞ r sk ( l ) e − jωlT s (45)According to the sampling theorem, ˜ S k ( ω ) and S k ( ω ) arerelated as follows S k ( ω ) = 1 T s + ∞ X n = −∞ ˜ S k (cid:18) ω + n πT s (cid:19) (46)Under Assumption A5, i.e. f s ≥ B ≥ B k , the power spectrum ˜ S k ( ω ) can be perfectly recovered by filtering S k ( ω ) with abandpass filter, i.e. ˜ S k ( ω ) = (cid:26) T s S k ( ω ) , ω ∈ [ ω k − πf s , ω k + πf s ]0 , ω / ∈ [ ω k − πf s , ω k + πf s ] . (47) Given the estimated factor matrix ˆ R , the DTFT of theautocorrelation sequence { r sk ( l ) } can be approximated as ˆ S k ( ω ) = L X l = − L ˆ r sk ( l ) e − jωlT s (48)When L is chosen to be sufficiently large, the estimation errordue to the time lag truncation is negligible. Also, althoughthere exists a phase ambiguity between the estimated auto-correlation sequence ˆ r k and the true autocorrelation sequence r k , this phase ambiguity can be removed by noting that thepower spectrum S k ( ω ) is real and non-negative. In addition,the power spectrum of each source is automatically paired withits associated DoA and carrier frequency due to the reason thatboth ˆ R and ˆ A experience a common permutation operation.VI. U NIQUENESS OF
CP D
ECOMPOSITION
We see that the uniqueness of the CP decomposition iscrucial to our proposed method. It is well known that theessential uniqueness of CP decomposition can be guaranteedby Kruskal’s condition [25]. Let k X denote the k-rank of amatrix X , which is defined as the largest value of k X suchthat every subset of k X columns of the matrix X is linearlyindependent. We have the following theorem concerning theuniqueness of CP decomposition. Theorem 1:
Let ( X , Y , Z ) be a CP solution which decom-poses a third-order tensor X ∈ C d × d × d into p rank-onearrays, where X ∈ C d × p , Y ∈ C d × p , and Z ∈ C d × p .Suppose the following Kruskal’s condition k X + k Y + k Z ≥ p + 2 (49)holds and there is an alternative CP solution ( c X , b Y , b Z ) whichalso decomposes X into p rank-one arrays. Then we have c X = X ΠΛ x , b Y = Y ΠΛ y , and b Z = Z ΠΛ z , where Π isa unique permutation matrix and Λ x , Λ y , and Λ z are uniquediagonal matrices such that Λ x Λ y Λ z = I . Proof:
A rigorous proof can be found in [26].Note that Kruskal’s condition cannot hold when R = 1 .However, in that case the uniqueness has been proven byHarshman [27]. Kruskal’s sufficient condition is also necessaryfor R = 2 and R = 3 , but not for R > [26].From the above theorem, we know that if k R + k A + k A ∗ ≥ K + 2 (50)then the CP decomposition of R x is essentially unique. Since A ∗ is the complex conjugate of A , we only need to examinethe k-ranks of A and R .Note that the ( n, k ) th entry of A is given by a nk = e − j (( n − τ k ω k +∆ n ω k ) (51)which is a function of the time delay factor ∆ n . It is notdifficult to design a set of time delay factors { ∆ n } such that k A = K . For example, we divide N antennas into two groups S = { , . . . , K } and S = { K + 1 , . . . , N } . We set the delayfactors in the first group to be linearly proportional to n − ,i.e. ∆ n = ( n − ν for n ∈ S , where ν ≥ is a constant. Inthis case, the first K rows of A form a Vandermonde matrix: A [1: K, :] = Vand ( τ ω + νω , . . . , τ K ω K + νω K ) (52) where Vand ( φ , . . . , φ K ) is defined asVand ( φ , . . . , φ K ) , e − j (0 φ ) . . . e − j (0 φ K ) e − j (1 φ ) . . . e − j (1 φ K ) ... . . . ... e − j (( K − φ ) . . . e − j (( K − φ K ) Thus A is full column rank with k A = K as long as { ω k τ k + νω k } are distinct from each other. If we set ν = 0 , we onlyneed { ω k τ k } , i.e. { ω k cos θ k } , are distinct from each other. Forthe case where the quantities { ω k cos θ k } for different sourcesignals may be identical, we can set ν = 0 , in which case westill have k A = K provided that the carrier frequencies { ω k } are mutually different. In other words, as long as AssumptionA2 is satisfied, we can always set an appropriate value of ν to ensure k A = K . For other more general choices of timedelay factors { ∆ n } , it can be numerically checked that the k-rank of A still equals to K with a high probability, althougha rigorous proof is difficult.Since we have k A = K , we only need k R ≥ in orderto satisfy Kruskal’s condition. This condition k R ≥ is metif every two columns of R are linearly independent. Notethat the k th column of R , r k , is a truncated autocorrelationsequence of the k th modulated signal s k ( t ) e jω k t . Clearly,if the baseband signals { s k ( t ) } have distinct power spectra,then any two columns of R are linearly independent, whichimplies k R ≥ . In practice, since source signals usually havedifferent bandwidths, the diverse power spectra condition canbe easily satisfied. Even if the baseband signals { s k ( t ) } haveidentical power spectra, the autocorrelation sequences of anytwo modulated signals { s k ( t ) e jω k t , s k ( t ) e jω k t } could stillbe linearly independent as long as their carrier frequenciessatisfy mod {| ω k − ω k | , πf s } 6 = 0 (53)The above condition ensures that autocorrelation sequences { r k } of different modulated signals have distinct exponentialterms { e jω k lT s } (see (11)). Due to the randomness of locationsof the carrier frequencies, the condition (53) is very likely tobe satisfied in practice. As a result, we have k R ≥ .VII. CRB A NALYSIS
In this section, we develop Cram´er-Rao bound (CRB)results for the joint DoA, carrier frequency, and power spectraestimation problem considered in this paper. As is well known,the CRB is a lower bound on the variance of any unbiasedestimator [28]. It provides a benchmark for evaluating theperformance of our proposed method. In addition, the CRBresults illustrate the behavior of the resulting bounds, whichhelps understand the effect of different system parameters,including the noise power σ , the number of antennas N andthe number of samples N s , on the estimation performance. A. Signal Model
Recall that the analog signal at each antenna is sampledwith a sampling rate f s = 1 /T s . The sampled signal at the n th antenna can be written as (cf. (7)) x n ( lT s ) = K X k =1 s k ( lT s ) e jω k ( lT s − ( n − τ k − ∆ n ) + w n ( lT s )= K X k =1 a nk s k ( lT s ) e jω k ( lT s ) + w n ( lT s ) (54)The above signal model can be rewritten in a vector-matrixform as x l = As l + w l , l = 0 , . . . , N s − (55)where A is defined in (20), x l , [ x ( lT s ) . . . x N ( lT s )] T , w l , [ w ( lT s ) . . . w N ( lT s )] T , and s l , h s ( lT s ) e jω ( lT s ) . . . s K ( lT s ) e jω K ( lT s ) i T Suppose we collect a total number of N s ( l = 0 , . . . , N s − )samples. The received signal can thus be expressed as X = AS + W (56)where X , [ x . . . x N s − ] S , [ s . . . s N s − ] W , [ w . . . w N s − ] . Let x , vec ( X T ) , where vec ( Z ) denotes a vectorizationoperation which stacks the columns of Z into a single columnvector. We have x = ¯ As + w (57)where x , vec( X T ) , w , vec( W T ) s , vec( S T ) , ¯ A , A ⊗ I N s (58)in which I n denotes an n × n identity matrix. We assumethat w ∼ CN ( , σ I N · N s ) and s ∼ CN ( , R s ) followa circularly-symmetric complex Gaussian distribution, where R s denotes the source covariance matrix which needs tobe estimated along with other parameters. Note that in ourproposed algorithm, the additive noise w and the sourcesignal s are not restricted to be circularly-symmetric complexGaussian. Here we make such an assumption in order tofacilitate the CRB analysis.Under the assumption that w and s are circularly-symmetriccomplex Gaussian random variables, we can readily verifythat x also follows a circularly-symmetric complex Gaussiandistribution, i.e. x ∼ CN ( , R x ) , where R x , E (cid:2) xx H (cid:3) = E h ¯ Ass H ¯ A H i + E (cid:2) ww H (cid:3) = ¯ AR s ¯ A H + σ I NN s (59)From Assumption A1, we know that R s is a block diagonalmatrix, i.e. R s = diag( P , . . . , P K ) . (60) where P k , E (cid:2) ˜ s k ˜ s Hk (cid:3) denotes the autocorrelation matrix ofthe k th signal, and ˜ s k is the transpose of the k th row of S ,i.e. ˜ s k , h s k (0 T s ) e jω k (0 T s ) . . . s k (( N s − T s ) e jω k (( N s − T s ) i T Also, in Assumption A1, each source is assumed to be wide-sense stationary. Therefore the autocorrelation matrix P k isa Hermitian-Toeplitz matrix. Here Toeplitz means that it hasdiagonal-constant entries, i.e. each descending diagonal fromleft to right is constant. Let p k denote the constant for elementslocated on the main diagonal, and p kl , l ≥ , denote theconstant for elements located on the l th diagonal below themain diagonal of P k . Let T l , (cid:18) I N s − l (cid:19) ∈ R N s × N s and T − l , (cid:18) I N s − l (cid:19) ∈ R N s × N s The autocorrelation matrix P k can thus be expressed as P k = p k I N s + L X l =1 (cid:2) p kl T − l + ( p kl ) ∗ T l (cid:3) (61)where L is chosen to be sufficiently large to ensure p kl = 0 for l > L . From (61), we can see that P k is characterized byparameters p k , (cid:2) p k ℜ ( p k ) . . . ℜ ( p kL ) ℑ ( p k ) . . . ℑ ( p kL ) (cid:3) (62)As a result, R s is characterized by parameters p , [ p . . . p K ] (63)On the other hand, notice that ¯ A is a parameterized matrix,with each column of A determined by the DoA and the carrierfrequency of each source, i.e. { θ k , ω k } . Unfortunately, thevalue ranges for the DoA and the carrier frequency differ byorders of magnitude, which may cause numerical instabilityin computing the CRB matrix. To address this difficulty, we,instead, analyze the CRB for the following two parameters { ξ k , ψ k } defined as ξ k , ω k τ k ψ k , ω k /c (64)where c is a parameter appropriate chosen (e.g. c = 10 ) suchthat values of ξ k and ψ k roughly have the same scale. Also,we define ξ , [ ξ . . . ξ K ] ψ , [ ψ . . . ψ K ] We see that the complete set of parameters to be estimatedinclude α , (cid:2) ξ ψ p σ (cid:3) (65)Recall that x follows a complex Gaussian distribution withzero mean and covariance matrix R x . Therefore the log-likelihood function of α can be expressed as L ( α ) ∝ − ln | R x | − x H R − x x (66) B. Calculation of The CRB Matrix
Since the random vector x follows a circularly-symmetriccomplex Gaussian distribution, we can resort to the Slepian-Bangs formula [29], [30] to compute the Fisher informationmatrix (FIM). According to the Slepian-Bangs formula, the ( i, j ) th element of the FIM Ω is calculated as Ω ij = tr (cid:18) R − x ∂ R x ∂α i R − x ∂ R x ∂α j (cid:19) (67)where α i and α j denote the i th and the j th entries of α ,respectively.By utilizing the structures of ¯ A and R s (cf. (58) and (60)), R x can be expressed as R x =[ a ⊗ I N s . . . a K ⊗ I N s ] · diag( P , . . . , P K ) · [ a ⊗ I N s . . . a K ⊗ I N s ] H + σ I N · N s =[ a ⊗ P . . . a K ⊗ P K ] · [ a ⊗ I N s . . . a K ⊗ I N s ] H + σ I N · N s = K X k =1 ( a k a Hk ) ⊗ P k + σ I N · N s (68)where a k , defined in (16), is the k th column of A .We first compute the partial derivative of R x with respectto ξ k and ψ k . From (13) and the definition of { ξ k , ψ k } , wecan write a nk = e − j (( n − ξ k + c ∆ n ψ k ) (69)Thus we have ∂ a k ∂ξ k = − j · diag(0 , . . . , N − · a k (70) ∂ a k ∂ψ k = − j · c · diag(∆ , . . . , ∆ N ) · a k (71)Combining (68) and (70)–(71), we have ∂ R x ∂ξ k = (cid:18) ∂ a k ∂ξ k a Hk + a k ∂ a Hk ∂ξ k (cid:19) ⊗ P k (72) ∂ R x ∂ψ k = (cid:18) ∂ a k ∂ψ k a Hk + a k ∂ a Hk ∂ψ k (cid:19) ⊗ P k . (73)Similarly, we can obtain the partial derivatives with respect toother parameters as follows ∂ R x ∂ ( p k ) = (cid:0) a k a Hk (cid:1) ⊗ (cid:18) ∂ P k ∂ ( p k ) (cid:19) = (cid:0) a k a Hk (cid:1) ⊗ I N s (74) ∂ R x ∂ ( ℜ ( p kl )) = (cid:0) a k a Hk (cid:1) ⊗ (cid:18) ∂ P k ∂ ( ℜ ( p kl )) (cid:19) = (cid:0) a k a Hk (cid:1) ⊗ ( T − l + T l ) (75) ∂ R x ∂ ( ℑ ( p kl )) = (cid:0) a k a Hk (cid:1) ⊗ (cid:18) ∂ P k ∂ ( ℑ ( p kl )) (cid:19) = (cid:0) a k a Hk (cid:1) ⊗ ( j T − l − j T l ) (76)and ∂ R x ∂ ( σ ) = I N · N s . (77)After obtaining the FIM Ω , the CRB can be calculated as [28] CRB( α ) = Ω − . (78) Frequency (Hz) × D O A Real: Frequency-DOAEstimated: Frequency-DOA (a) True and estimated carrier frequencies and DoAs.
Frequency (Hz) × P o w e r s pe c t r u m (b) Original power spectra of sources. Frequency (Hz) × P o w e r s pe c t r u m (c) Estimated power spectra of sources.Fig. 4. True and estimated carrier frequencies, DoAs, and power spectra ofsources, SNR = 5 dB. VIII. S
IMULATION R ESULTS
In this section, we carry out experiments to illustrate theperformance of our proposed method. In our simulations, weset f nyq = 1 GHz. The distance between two adjacent antennas, d , is set equal to d = 0 . × C/f nyq in order to meet thecondition in Assumption A6. The number of antennas is setto N = 8 , and for simplicity, the time delay factors are set as ∆ n = (cid:26) s , n = 1 , . . . , N/ − s , n = N/ , . . . , N (79)With this setup, the condition (6) can be satisfied for n = N/ − . The signal-to-noise ratio (SNR) is defined asSNR , E [ | s ( t ) | ] σ (80)We first consider the case in which K = 3 uncorrelated,wide-sense stationary sources spreading over the wide fre-quency band (0 , MHz impinge on a ULA of N antennas.The DoAs of these three sources are given respectively by θ = 2 . , θ = 1 . , and θ = 0 . . The carrier frequen-cies and bandwidths associated with these sources are set to f = 152 MHz, f = 323 MHz, f = 432 MHz, B = 20 MHz, B = 20 MHz, and B = 15 MHz. The complex baseband Frequency (Hz) × D O A Real: Frequency-DOAEstimated: Frequency-DOA (a) True and estimated carrier frequencies and DoAs.
Frequency (Hz) × P o w e r s pe c t r u m Source 1Source 2 (b) Original power spectra of sources.
Frequency (Hz) × P o w e r s pe c t r u m Source 1Source 2 (c) Estimated power spectra of sources.Fig. 5. Estimated carrier frequencies, DoAs and power spectra for sourcesthat have partial spectral overlap, SNR = 20 dB. signals are generated by passing the complex white Gaussiannoise through low-pass filters with different cutoff frequencies.Also, the number of data samples used for calculating thecorrelation matrices is set to N s = 10 . The sampling rate f s is chosen to be f s = 28 MHz, which is slightly higherthan the minimum sampling rate f s ≥ B = max { B , B , B } required for perfect recover of the power spectrum of the widefrequency band. The SNR is set to 5dB. Fig. 4(a) shows thetrue (marked with ‘ (cid:3) ’) and the estimated (marked with ‘ + ’)carrier frequencies and DoAs for the three sources. We cansee that the estimated carrier frequencies and DoAs coincidewith the groundtruth well. Fig. 4(b) and Fig. 4(c) respectivelydepict the original power spectrum and the estimated powerspectrum of the wide frequency band. It can be observed thatour proposed method, even with a low SNR and a samplingrate far below the Nyquist rate, is able to accurately identifythe locations of the occupied bands.Next, we examine the scenario where frequency bands of thenarrowband sources overlap each other. Set K = 2 . The DoAsof these two sources are given respectively by θ = 2 . and θ = 0 . . The carrier frequencies and bandwidths associatedwith these two sources are set to f = 151 . MHz, f = 161 . MHz, B = 20 MHz, and B = 10 MHz. The powerspectra associated with the two sources are shown in Fig. 5(b),from which we can see that the two sources partially overlap inthe frequency domain. The number of data samples N s and thesampling rate f s remain the same as in the previous example.The SNR is set to 20dB. The estimated carrier frequencies,DoAs, and the power spectra of the two sources are plotted inFig. 5(a) and Fig. 5(c). We see that our proposed method workswell for sources with partially overlapping frequency bands.This example shows that our proposed method not only canperform wideband spectrum sensing, but also has the abilityto blindly separate power spectra of sources that have partialspectral overlap.To better evaluate the performance of our proposed method,we calculate the mean square errors (MSEs) for the followingsets of parametersMSE ( ψ ) = K X k =1 | ψ k − ˆ ψ k | MSE ( ξ ) = K X k =1 | ξ k − ˆ ξ k | MSE ( θ ) = K X k =1 | θ k − ˆ θ k | Recalling that in our analysis, instead of concerning { θ k , ω k } ,we define two new parameters ξ k , ω k τ k and ψ k , ω k /c andderive the CRB for { ξ k , ψ k } in order to avoid the numericalinstability issue. The MSEs of the sets of parameters { ξ k , ψ k } are also included to compare with their associated CRB results.The estimation accuracy of the carrier frequencies is quantifiedby the normalized mean square error (NMSE) defined asNMSE ( ω ) = K X k =1 | ω k − ˆ ω k | | ω k | In this example, we set the number of sources K = 2 . Theparameters associated with these two sources are given as: f = 152 MHz, f = 437 MHz, B = 126 KHz, B = 63 KHz, θ = π/ , and θ = π/ . The sampling rate is set to f s = 1 . MHz. Fig. 6 and Fig. 7 depict the MSEs/NMSEs ofrespective sets of parameters vs. the number of samples N s ,where we set SNR = 5 dB and SNR = 15 dB, respectively.MSE/NMSE results are averaged over 1000 independent runs,where the baseband complex source signals are randomlygenerated for each run. We see that our proposed method canachieve an estimation accuracy close to the CRBs by usingonly a small number of data samples, e.g. N s = 200 . In Fig.8, we plot the MSEs/NMSEs of different sets of parametersas a function of the SNR, where N s = 300 data samples areused. We see that under a moderately high SNR, say, SNR =15 dB, our proposed method attains an accurate estimate ofthe DoAs/carrier frequencies with the MSE (NMSE) as lowas − . IX. C ONCLUSIONS
We considered the problem of joint wideband spectrumsensing and DoA estimation in this paper. To overcome
50 100 150 200 250 300
Number of samples per antenna (5dB) -4 -3 -2 -1 M SE MSE: ξ CRB: ξ (a)
50 100 150 200 250 300
Number of samples per antenna (5dB) -3 -2 -1 M SE MSE: ψ CRB: ψ (b)
50 100 150 200 250 300
Number of samples per antenna (5dB) -3 -2 -1 M SE MSE: θ (c)
50 100 150 200 250 300
Number of samples per antenna (5dB) -3 -2 -1 N M SE NMSE: ω (d)Fig. 6. MSEs and NMSE vs. the number of samples per antenna, where N = 8 and SNR = 5 dB. the sampling rate bottleneck, we proposed a phased-arraybased sub-Nyquist sampling architecture (termed as PASSAT)that is simpler in structure and easier for implementationas compared with existing sub-Nyquist receiver architectures.Based on the proposed receiver architecture, we developeda CP decomposition-based method for joint DoA, carrierfrequency, and power spectrum estimation. The conditionsfor exact recovery of the parameters and the power spectrumwere analyzed. Our analysis suggests that the perfect recoverycondition for our proposed method is mild: to recover thepower spectrum of the wide frequency band, we only need
50 100 150 200 250 300
Number of samples per antenna (15dB) -5 -4 -3 -2 M SE MSE: ξ CRB: ξ (a)
50 100 150 200 250 300
Number of samples per antenna (15dB) -4 -3 -2 -1 M SE MSE: ψ CRB: ψ (b)
50 100 150 200 250 300
Number of samples per antenna (15dB) -4 -3 -2 M SE MSE: θ (c)
50 100 150 200 250 300
Number of samples per antenna (15dB) -4 -3 -2 N M SE NMSE: ω (d)Fig. 7. MSEs and NMSE vs. the number of samples per antenna, where N = 8 and SNR = 15 dB. the sampling rate to be greater than the bandwidth of thenarrowband source signal which has the largest bandwidthamong all sources. In addition, even for the case where sourceshave partial spectral overlap, our proposed method is stillable to extract the DoA, carrier frequency, and the powerspectrum associated with each source signal. CRB analysis forour estimation problem was also carried out. Simulation resultsshow that our proposed method, with only a small number ofdata samples, can achieve an estimation accuracy close to theassociated CRBs. -10 -5 0 5 10 15 SNR (dB) -5 -3 -1 M SE MSE: ξ CRB: ξ (a) -10 -5 0 5 10 15 SNR (dB) -3 -1 M SE MSE: ψ CRB: ψ (b) -10 -5 0 5 10 15 SNR (dB) -4 -2 M SE MSE: θ (c) -10 -5 0 5 10 15 SNR (dB) -4 -2 N M SE NMSE: ω (d)Fig. 8. MSEs and NMSE vs. SNR (dB), where N = 8 and N s = 300 . R EFERENCES[1] E. Axell, G. Leus, E. G. Larsson, and H. V. Poor, “Spectrum sensingfor cognitive radio: State-of-the-art and recent advances,”
IEEE SignalProcessing Magazine , vol. 29, no. 3, pp. 101–116, May 2012.[2] H. Sun, A. Nallanathan, C.-X. Wang, and Y. Chen, “Wideband spectrumsensing for cognitive radio networks: a survey,”
IEEE Wireless Commu-nications , vol. 20, no. 2, pp. 74–81, Apr. 2013.[3] M. Mishali and Y. C. Eldar, “Blind multiband signal reconstruction:Compressed sensing for analog signals,”
IEEE Trans. Signal Processing ,vol. 57, no. 3, pp. 993–1009, Mar. 2009.[4] ——, “From theory to practice: Sub-Nyquist sampling of sparse wide-band analog signals,”
IEEE J. Sel. Topics Signal Process. , vol. 4, no. 2,pp. 375–391, Apr. 2010.[5] M. Mishali, Y. C. Eldar, O. Dounaevsky, and E. Shoshan, “Xampling: Analog to digital at sub-Nyquist rates,”
IET Circuits, Devices andSystems , vol. 5, no. 1, pp. 8–20, Jan. 2011.[6] M. Wakin, S. Becker, E. Nakamura, M. Grant, E. Sovero, D. Ching,J. Yoo, J. Romberg, A. Emami-Neyestanak, and E. Cand`es, “A nonuni-form sampler for wideband spectrally-sparse environments,”
IEEE Jour-nal on Emerging and Selected Topics in Circuits and Systems , vol. 2,no. 3, pp. 516–529, Sept. 2012.[7] E. Cand´es, J. Romberg, and T. Tao, “Robust uncertainty principles: exactsignal reconstruction from highly incomplete frequency information,”
IEEE Trans. Information Theory , vol. 52, no. 2, pp. 489–509, Feb. 2006.[8] D. L. Donoho, “Compressive sensing,”
IEEE Trans. Information Theory ,vol. 52, no. 4, pp. 1289–1306, Apr. 2006.[9] D. D. Ariananda and G. Leus, “Compressive wideband power spectrumestimation,”
IEEE Trans. Signal Processing , vol. 60, no. 9, pp. 4775–4789, Sept. 2012.[10] C.-P. Yen, Y. Tsai, and X. Wang, “Wideband spectrum sensing based onsub-Nyquist sampling,”
IEEE Trans. Signal Processing , vol. 61, no. 12,pp. 3028–3040, June 2013.[11] D. Cohen and Y. C. Eldar, “Sub-Nyquist sampling for power spectrumsensing in cognitive radios: A unified approach,”
IEEE Trans. SignalProcessing , vol. 62, no. 15, pp. 3897–3910, Aug. 2014.[12] M. D. Zoltowski and C. P. Mathews, “Real-time frequency and 2-D angleestimation with sub-Nyquist spatio-temporal sampling,”
IEEE Trans.Signal Processing , vol. 42, no. 10, pp. 2781–2794, Oct. 1994.[13] S. Stein, O. Yair, D. Cohen, and Y. C. Eldar, “Joint spectrum sensingand direction of arrival recovery from sub-Nyquist samples,” in
IEEEInternational Workshop on Signal Processing Advances in WirelessCommunications , Stockholm, Sweden, 2015, pp. 331–335.[14] A. N. Lemma, A.-J. van der Veen, and E. F. Deprettere, “Joint angle-frequency estimation using multi-resolution ESPRIT,” in
IEEE Interna-tional Conference on Acoustics, Speech, and Signal Processing , Seattle,Washington, USA, 1998, pp. 1957–1960.[15] ——, “Analysis of joint angle-frequency estimation using ESPRIT,”
IEEE Trans. Signal Processing , vol. 51, no. 5, pp. 1264–1283, May2003.[16] D. D. Ariananda and G. Leus, “Compressive joint angular-frequencypower spectrum estimation,” in
IEEE European Signal ProcessingConference , Marrakech, Morocco, 2013, pp. 1–5.[17] A. A. Kumar, S. G. Razul, and C.-M. S. See, “An efficient sub-Nyquistreceiver architecture for spectrum blind reconstruction and directionof arrival estimation,” in
IEEE International Conference on Acoustics,Speech, and Signal Processing , Florence, Italy, 2014, pp. 6781–6785.[18] ——, “Spectrum blind reconstruction and direction of arrival estimationat sub-Nyquist sampling rates with uniform linear array,” in
IEEEInternational Conference on Digital Signal Processing , Singapore, 2015,pp. 670–674.[19] S. S. Ioushua, O. Yair, D. Cohen, and Y. C. Eldar, “CaSCADE: Com-pressed carrier and DOA estimation,”
IEEE Trans. Signal Processing ,vol. 65, no. 10, pp. 2645–2658, May 2017.[20] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,”
SIAM review , vol. 51, no. 3, pp. 455–500, 2009.[21] A. Lavrenko, F. R¨omer, S. Stein, D. Cohen, G. D. Galdo, R. S. Thom¨a,and Y. C. Eldar, “Spatially resolved sub-Nyquist sensing of multibandsignals with arbitrary antenna arrays,” in
IEEE International Workshopon Signal Processing Advances in Wireless Communications , Edinburgh,UK, 2016, pp. 1–5.[22] J. A. Bazerque, G. Mateos, and G. B. Giannakis, “Rank regularizationand bayesian inference for tensor completion and extrapolation,”
IEEETrans. Signal Processing , vol. 61, no. 22, pp. 5689–5703, November2013.[23] P. Rai, Y. Wang, S. Guo, G. Chen, D. Dunson, and L. Carin, “ScalableBayesian low-rank decomposition of incomplete multiway tensors,” in
Proc. of the 31st Inter. Conf. on Mach. Learning (ICML-14) , vol. 32,Beijing, China, 2014, pp. 1800–1808.[24] Q. Zhao, L. Zhang, and A. Cichocki, “Bayesian CP factorization ofincomplete tensors with automatic rank determination,”
IEEE Trans.Pattern Anal. Mach. Intell. , vol. 37, no. 9, pp. 1751–1763, Sept. 2015.[25] J. B. Kruskal, “Three-way arrays: rank and uniqueness of trilineardecompositions, with application to arithmetic complexity and statistics,”
Linear Algebra and its Applications , vol. 18, no. 2, pp. 95–138, 1977.[26] A. Stegeman and N. D. Sidiropoulos, “On Kruskal’s uniqueness condi-tion for the Candecomp/Parafac decomposition,”
Linear Algebra and itsApplications , vol. 420, no. 2-3, pp. 540–552, Jan. 2007.[27] R. A. Harshman, “Determination and proof of minimum uniquenessconditions for PARAFAC1,”
UCLA Working Papers in Phonetics , no. 22,pp. 111–117, 1972. [28] S. M. Kay, Fundamentals of Statistical Signal Processing: EstimationTheory . Upper Saddle River, NJ: Prentice Hall, 1993.[29] P. Stoica and A. Nehorai, “Performance study of conditional and uncon-ditional direction-of-arrival estimation,”
IEEE Trans. Acoust., Speech,Signal Processing , vol. 38, no. 10, pp. 1783–1795, Oct. 1990.[30] P. Stoica, E. G. Larsson, and A. B. Gershman, “The stochastic CRB forarray processing: a textbook derivation,”