ROBUST PARAMETRIC POWER SPECTRAL DENSITY (PSD) CONSTRUCTION

Abstract

This disclosure describes techniques for constructing power spectral density (PSD) maps representative of the distribution of radio frequency (RF) power as a function of both frequency and space (geographic location). For example, the disclosure describes techniques for construction PSD maps using robust basis pursuit forms of signal expansion.

Description

TECHNICAL FIELD

The invention relates to wireless communication and, more specifically, to spectrum cartography.

BACKGROUND

All wireless transmissions use a portion of the radio frequency spectrum. Cellular phones, broadcast television, satellite, and short-distance wireless networks such as Bluetooth and wireless local area networks (WLAN) utilize different portions of the Wi-Fi, for example, typically use wireless frequency spectrum. Often it is important to coordinate the use of the various technologies and frequency ranges to ensure that the technologies do not interfere with each other or with planned future services.

SUMMARY

This disclosure describes techniques for constructing power spectral density (PSD) maps representative of the distribution of radio frequency (RF) power as a function of both frequency and space (geographic location). For example, the disclosure describes techniques for construction PSD maps using robust basis pursuit form of signal expansion.

To further address the challenges encountered with this multi-dimensional sensing vision, the present paper presents a collaborative sensing scheme whereby CRs cooperate to localize the actively transmitting primary users (Pus) and estimate their PSD across space in the presence of model uncertainties.

One example of the techniques described herein is a novel sensing scheme based on a parsimonious system model accounting for the scarce presence of active PUs in the same frequency band(s), in the monitored area, due to mutual interference. Using a virtual grid-based approach for the potential PU transmitter locations, a form of spatial-domain sparsity emerges because actual PU transmitters are present in only few of the potential (grid) locations. A basis expansion model is then adopted to approximate the PU transmit-PSD distribution in frequency, which renders the sensing objective tantamount to estimating the PSD basis coefficients corresponding to each grid point. Since individual PU transmissions are narrow-band relative to the large swath of frequencies a CR can sense, only few of the PSD basis coefficients are nonzero—a fact giving rise to frequency-domain sparsity.

This parsimonious system model thus entails a form of hierarchical dual-domain sparsity in the PSD basis coefficients that are to be estimated, in the sense that groups of coefficients corresponding to locations with no PUs will be collectively zero. In addition, some of the basis coefficients within groups corresponding to active PU locations will be zero. Capitalizing on this form of hierarchical sparsity, a group sparse regression problem is formulated, which is to be solved centrally by either a fusion center (FC), or, individual CRs using their measured PSDs. A novel low-complexity algorithm for solving such a problem is developed using the alternating direction method of multipliers (ADMoM).

One challenge addressed by the techniques is acquiring the grid-to-CR channel gains present in the underlying regression matrix. One way to acquire such information is through the channel gain cartography approach. However, possible inaccurate channel gains or adoption of a shadowing-agnostic path loss-only model could deteriorate the performance of the sensing algorithm. Also, a grid-based approach introduces itself possible model offsets, as the actual PU locations may not coincide with points of the grid. To account for these uncertainties, a robust version of the group sparse (GS) least-absolute-shrinkage-and-selection operator (Lasso) is described herein. One technique described herein is an extension of the sparse total least-squares (TLS) framework to incorporate the hierarchical sparsity inherent to this sensing application. Combining the merits of Lasso, group Lasso, and TLS, the proposed group sparse (GS-)TLS approach described herein yields hierarchically-sparse PSD estimates that are also robust to model uncertainties induced by the random channel, grid offsets, and basis approximation errors. In spite of the non-convexity of the proposed GS-TLS criterion, an iterative solver with guaranteed convergence to at least a locally-optimal solution is described.

Additional factors that may compromise accuracy of PSD estimates at the CRs, are abrupt changes in shadow fading that may be due to, e.g., moving obstacles or moving CRs, and, possible failures of the sensing modules themselves. A robust GS-TLS formulation is proposed here, that is capable of discerning and removing such so-called model outliers, which in turn leads to reliable PSD estimates. However, sorting out unreliable measurements not only promotes estimation accuracy, but also leads to self-healing and re-organization mechanisms for the CRs network.

This disclosure introduces the basis expansion model, and describes the PSD observations used for the model fitting approach. A centralized algorithm for solving GS-Lasso problems is described, whereas perturbations in the channel (regression) matrices are considered. An outlier-resilient sensing algorithm is devised, numerical results are provided.

In one example, a method comprises sensing local radio-frequency (RF) interference spectrum at each of a plurality of sensors positioned at a plurality of locations within a geographic region, and computing a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location. The basis expansion model is computed as a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band. Moreover, the coefficients are computed of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso).

In another example, a system comprises a plurality of sensors to sense a local radio-frequency (RF) interference spectrum at each of a plurality of locations within a geographic region. A processor computes a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location. The processor computes the basis expansion model is computed as a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band. The processor computing the coefficients of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso).

In another example, a mobile device comprises a sensor to sense a local radio-frequency (RF) interference spectrum at locations within a geographic region. The mobile device further comprises a processor that computes a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location. The processor computes the basis expansion model is computed as a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band. The processor computes the coefficients of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso).

The details of one or more examples are set forth in the accompanying drawings and the description below.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a block diagram illustrating a system in which a plurality of sensors and a centralized computing center cooperate to construct PSD maps.

FIG. 2 is a block diagram illustrating a system in which a plurality of sensors executed distributed techniques to construct PSD maps.

FIG. 3 illustrates a basis expansion model with overlapping raised cosine pulses.

FIG. 4 shows an example CR network topology.

FIGS. 5A, 5B represent PSD maps.

FIG. 6 illustrates evolution of the entries of Θ^(j).

FIG. 7 illustrates a true PSD map with shadowing propagation effects.

FIGS. 8A, 8B illustrate estimated PSD maps of FIG. 7

FIG. 9 illustrates a true PSD map with deep shadowing caused by an obstacle.

FIGS. 10A, 10B illustrate estimated PSD maps of FIG. 9

FIG. 11 illustrates an example CR topology and PU transmitter locations for numerical testing of the described techniques.

FIG. 12 is a graph illustrating interference power using the batch technique.

FIG. 13 is a graph illustrating interference power using the online technique.

FIG. 14 is a graph illustrating average normalized RMSE.

FIG. 15 is a graph illustrating one-time slot prediction.

FIG. 16 shows a detailed example of various devices that may be configured to execute program code to practice some embodiments in accordance with the current disclosure.

DETAILED DESCRIPTION

FIG. 1 is a block diagram illustrating a system 10 in which a plurality of sensors 12 (e.g., cognitive radios (CR) in this example) deployed within a spatial region. Each of CRs 12 sense the ambient interference spectrum from other RF-enabled devices 15 within its surrounding region and communicate the sensed observations to one another via messages 14. CRs 12 collaborate to construct one or more power spectral density (PSD) maps based on the exchanged PSD observations. In general, the PSD maps are representative of the distribution of radio frequency (RF) power as a function of both frequency and spatial location within the geographic region. In the example of FIG. 2, system 10 includes a centralized fusion center (FC) 16 that performs the techniques described herein to compute the PSD maps 17 based on the sensed observations relayed to the fusion center via the CRs. In one example of system 10, FC 16 includes and maintains PSD maps 17 within a database along with location data for each of CRs 12. Each location may, for example, be represented as a position vector within the geographic region. A dedicated control channel may be established by which CRs 12 exchange PSD observations via messages 14, which are ultimately related to FC 16 by those CRs in communication with the FC.

In general, FC 16 is a computing system of one or more computing devices that applies the techniques described herein to compute one or more PSD maps for the geographic region. In one example, FC 16 applies robust parametric basis pursuit as a form of signal expansion for construction PSD maps based on the observations. FC may, for example, compute the basis expansion model from the sensed RF interference spectrum at each of the CRs 12 to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location. In general, a basis expansion model consists of a superposition of shifted and scaled versions of reference basis functions. The reference basis functions represent the frequency distribution of the RF power, i.e., RF power present at different frequency slots. The basis expansion model may be computed as a plurality of reference basis functions having corresponding coefficients, each of the reference basis functions representing a power emitted by an RF-enabled device 15 on a corresponding frequency band. As used herein, a paramentric basis expansion model refers to a basis expansion model where scaling coefficients of the reference basis functions are computed as scalars.

Moreover, as described below, FC 16 may apply a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso) in computing the coefficients of the basis expansion model. In general, least-absolute-shrinkage-and-selection operator (LASSO) is one example of a compressive sensing operator. Group sparse (GS) LASSO techniques are described in which a group sparsity structure for the basis expansion model is encoded and enforced. In this way, as further described below, a parametric PSD map representative of the power distribution of the RF power can be recovered, i.e., computed, with fewer measurements. The GS Lasso operator, in effect, is applied to select which coefficients in the model are non-zero and are to be estimated by FC 16. For example, in accordance with the techniques described herein, the group sparse Lasso (GS Lasso) identifies the locations in space where and RF-enabled device 15 is actively transmitting, and the frequency band(s) where transmission occurs. Specifically, by application of the GS Lasso described herein, non-zero coefficients in the basis expansion model for a transmitter's location correspond to occupied bands. In this way, FC 16 may apply the GS Lasso operator described herein to enforce group sparsity within the basis expansion model, thereby reflecting sparsity within system 10 that may arise due to geographic sparsity of sensors RF-enabled devices 15 within the measured environment and sparsity in the frequency domain since individual transmissions may occupy small portions of the spectrum. In some cases, FC 16 may use a most parsimonious sparse signal expansion using an overcomplete basis set may be used to constructing the PSD maps.

Consider an incumbent PU system comprising N_s transmitters (sources) located in a geographical area A⊂R². Their activity over a frequency band B is to be monitored via cooperation of N_r CRs, also located in A. Let S:={x_sεA}_s=1^NS denote the PU locations. The sensing objective is to localize the PU sources, and reveal available portions of B for the CRs to transmit opportunistically.

Let u_s(t) be the (unknown) signal transmitted by PU s at time t. Then, the signal received at CR position x_r at time t can be expressed as

$y_r\left(t\right):=\sum _{s=1}^{N_s}\sum _{l=0}^{L_{\mathrm{sr}}-1}h_{x_s\to x_r}\left(t;l\right)u_s\left(t-l\right)+v\left(t\right)$

where h_xs→xr(t;l) is the l-th tap of the time-varying channel impulse response of the link x_s→x_r, and v(t) denotes the additive white noise. Regarding signals {u_s(t)}, the following is assumed.

(As1) Sources {u_s(t)} are stationary, mutually uncorrelated, independent of {h_xs→xr(t;l)}, with vanishing correlation per channel coherence interval.

Consider approximating the transmit-PSD of PU s using the following basis expansion model:

$\begin{array}{cc}{\Phi }_s\left(f\right)=\sum _{v=1}^{N_b}{\Theta }_{\mathrm{sv}}b_v\left(f\right),s=1,2,\dots ,N_s& \left(1\right)\end{array}$

where N_b is assumed sufficiently large, and {Θ_sv}_v=1^Nb are nonnegative coefficients. Possible choices of {b_v(f)}_v=1^Nb include the set of non-overlapping rectangles of unit height spanning the bandwidth B of interest. In this case, each Θ_sV represents the power emitted by source s on the frequency band corresponding to the basis function b_v(f). Alternatively, overlapping raised cosine bases can be employed with support B_v=[f_v−(1+/2T_s),f_v+(1+/2T_s)], where is the roll-off factor and T_s the symbol period, can be employed; see also FIG. 3.

Channel {h_xs→xr(t;l)} can be decomposed as

$h_{x_s\to x_r}\left(t;l\right)={\left({\gamma }_{x_s\to x_r}s_{x_s\to x_r}\right)}^{\frac{1}{2f}}{\hspace{0.17em}}_{x_s\to x_r}\left(t;l\right),$

where γ_xs→xr stands for the path loss, s_xs→xr the temporally- and spatially-colored shadowing, and {f_xs→xr(t;l)} for the multi-path fast time-varying fading. The latter satisfies the following. (As2) Variables {f_xs→xr(t;l)} are complex Gaussian with zero mean and variance σ_f,sr², stationary with respect to t, and uncorrelated across the lag variable l and the spatial variables x_s and x_r. Without loss of generality, assume that

$\sum _{1=1}^{L_{\mathrm{sr}}}{\sigma }_{f,\mathrm{sr}}^2=1$

for every s and r.

Received samples {y_r(t)} are parsed into N-dimensional blocks, where N is chosen equal to (or smaller than) the coherence interval of the small-scale fading, over which h_xs→xr(t;l) remains approximately invariant with respect to (wrt) t. These data blocks are hereafter indexed by n, so that t=nN+m, with m=0, 1, . . . , N−1.

Shadowing and small-scale fading are characterized by different dynamics. The following is assumed regarding channel propagation and modeling.

(As3) The coherence interval of shadow fading exceeds that of f_xs→xr(t;l). Also, shadowing variations are sufficiently slower than the coherence interval of the PU signals.

Based on (As1)-(As3), it is possible to express the PSD measured at location x due to N_s simultaneous PU transmissions as

$\begin{array}{cc}{\Phi }_x\left(f\right)=\sum _{s=1}^{N_s}g_{x_s\to x}{\Phi }_s\left(f\right)+{\sigma }_v^2=\sum _{s=1}^{N_s}g_{x_s\to x}\sum _{v=1}^{N_b}{\Theta }_{\mathrm{sv}}b_v\left(f\right)+{\sigma }_v^2& \left(2\right)\end{array}$

where σ_v² denotes noise variance at the receiver, and g_xs→xr:=γ_xs→xrs_xs→xr the averaged channel gain.

As neither the number of PU sources nor their locations are known to the CRs, a set of N_g candidate transmit-PUs is postulated on a grid of locations G:={x_gεA}_g=1^Ng. Without prior knowledge of the area(s) where PU activity is more likely, the set of candidate locations G can be simply formed by discretizing A to the set of grid points G.

Define the N_b×1 vector Θ_g:=[Θ_g1, . . . , Θ_gNb]^T collecting the basis coefficients that correspond to location x_g, and let Θ:=[Θ₁^T, . . . , Θ_Ng^T]^T; also, let B_xr be the N×N_bN_g matrix B_xr:=[b_xr(f₁), . . . , b_xr(f_N)]^T, with b_x(f_k) having entries {g_xg→xr·b_v(f_k)}. Then, upon defining φ_xr:=[Φ_xr(f₁), . . . , Φ_xr(f_N)]^T, the received PSD at CR location x_r, sampled at frequencies {f_k}_k=1^N, can be compactly written as:

φ_xr=B_xrΘ+σ_v²1_N.  (3)

The sensing objective of revealing PU locations and the available portions (sub-bands) of B is tantamount to estimating Θ. To this end, CRs rely on the periodogram estimate of Φ_xr(f) at the sampling locations {x_r}_r=1^Nr, and N frequency bins {f_k}_k=1^N. The fast Fourier transform of samples {y_r(t)}, namely Y_r,N(n,f), and the periodogram ̂φ_xr,N(n;f):=(1/N)|Y_r,N(n,f)|² are computed per data block n. To average out small-scale fading effects, and allow for tracking of shadow fading as well as possible variations of the PUs' power spectra, the periodogram estimate at CR r is formed using an exponentially weighted moving average operation as

$\begin{array}{cc}{\hat{\Phi }}_{x_r}\left(\tau ,f\right):=\sum _{n=1}^{\tau }{\alpha }^{\tau -n}{\varphi }_{x_r},N^{\left(n,f\right)}& \left(4\right)\end{array}$

with αε(0,1] denoting the so-called forgetting factor. As shown gives an estimate of the PSD measured at point x_r, and frequency f as

̂Φ_xr(τ;f)=Φ_xr(τ;f)+e_xr,N(τ;f)  (5)

with (asymptotic) variance bounded as

$\begin{array}{c}\lim \\ N\to \infty \end{array}\mathrm{var}\left[e_{r,N}\left(\tau ;f\right)\right]\le \frac{3}{2}\left(1-\alpha \right){\Phi }_{x_r}^2\left(f\right).$

After dropping τ for notational brevity, let ̂φ_xr:=[̂Φ_x(f₁), . . . , ̂Φ_x(f_N)]^T.

Based on the linear model (3), the sensing objective is to estimate Θ from the received-PSD estimate

$\hat{\phi }:={\left[{\hat{\phi }}_{x_1},\dots ,{\hat{\phi }}_{x_{N_r}}\right]}^T$

gathered at CR locations {x_r}_r=1^Nr.

The following describes spectrum sensing via Group Sparse Lasso. The number of active PUs transmitting over the same spectral band in a given area is naturally limited by mutual interference. As a consequence, the number of PU sources (N_s) is far smaller than N_g, for a sufficiently dense grid. Absence of PU sources in most grid locations gives rise to a group sparsity of the vector Θ, since Θ_g=0_Nb for each of the locations x_g that are not occupied by a PU transmitter. In addition to space, sparsity in the vector Θ is also manifested in the frequency domain because of the parsimonious linear model (2). Compared to the possibly large swath of frequencies that the CRs can sense, individual PU transmissions typically occupy small portions of the spectrum (say, in the order of MHz). Sparsity in the frequency domain implies that individual entries within each group Θ_g are zero.

One criterion for estimating Θ is the (non-negative) least-squares (LS). However, LS fails to provide a parsimonious model estimate involving only the prominent variables. The Lasso and the so-called group Lasso on the other hand, were proposed to overcome such a limitation of LS. In the Lasso criterion, the LS cost is augmented with the l₁-norm |Θ|₁ to encourage sparsity at the single-coefficient level; while in the group Lasso, the regularization term

$R_G\left(\Theta \right):=\sum _{g=1}^{N_g}{{\Theta }_g}_2$

is referred to herein as a group sparse (GS) penalty that enforces group sparsity within the model.

Combining Lasso with group Lasso, the so-called group sparse (GS-)Lasso provides a parsimonious model estimate, where sparsity is accounted for both at the group- and at the single-coefficient levels. This hierarchical sparsity is possible by regularizing the conventional LS cost with the term

$R_G\left(\Theta \right):=\sum _{g=1}^{N_g}{{\Theta }_g}_2$

combined with |Θ|₁.

Taking also into account the non-negativity of PU power spectra, Θ can be estimated by solving the following sparse regression problem, where Θ indicates the location in space and the frequency band for any active transmitter (RF-enabled device 15):

$\begin{array}{cc}\hat{\theta }=\arg \underset{\theta \ge 0_{N_bN_g}}{\min }\left[\frac{1}{2}\sum _{r=1}^{N_r}{{\varphi }_{x_r}-B_{x_r}\theta -{\sigma }_v^21_N}_2^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)\right]& \left(6\right)\end{array}$

where the coefficient λ₁≧0 enforces sparsity at individual entries, whereas λ₂≧0 promotes group sparsity. For λ₁=0 (λ₂=0), (6) reduces to the Lasso (group Lasso) based estimate. In this way, the group sparse Lasso (GS Lasso) identifies the locations in space where and RF-enabled device 15 is actively transmitting, and the frequency band(s) where transmission occurs. Specifically, non-zero coefficients in the basis expansion model for a transmitter's location correspond to occupied bands. PU localization and PSD estimation was viewed as a sparse linear regression model; here, the formulation is considerably broadened by taking into account both individual and group sparsity.

To obtain regression matrices {B_xr}, the channel gains {g_xg→xr} need to be estimated. To this end, CRs can simply neglect shadowing, and in resort to the distance-dependent path loss model ĝ_xg→xr=min{1,(∥x_g−x_r∥₂/d₀)^−η}, where d₀ and η are preselected constants depending on the propagation environment. Alternatively, more sophisticated techniques can be employed. Perturbations in the regression matrices {B_xr} arising due to inaccurate channel estimation and grid-mismatch effects are dealt with below.

With respect to the PSD atlas, it is worth re-iterating that identifying the support of the vector Θ reveals not only the primary sub-bands occupied, but also the locations where the active PU transmitters reside. Complementing this information with either the PUs' channel gain maps or a simple path loss-based propagation model, CRs can readily reconstruct the PSD atlas; that is, estimate PSD maps at any location of the monitored area as

$\begin{array}{cc}{\hat{\Phi }}_x\left(f\right)=\sum _{g=1}^{N_g}{\hat{g}}_{x_g\to x}\sum _{v=1}^{N_b}{\hat{\Theta }}_{\mathrm{gv}}b_v\left(f\right),\forall x\in A& \left(7\right)\end{array}$

with ĝ_xg→x the estimate of g_xg→x. Having available estimates of the PSD map across space per frequency band (hence the term atlas), CRs can adjust their transmit power to prevent harmful interference inflicted to the PUs. In fact, the positions of potential PU receivers can be deduced from the PSD atlas; and thus, CR transmission powers can be properly adapted.

ADMoM-based Solver: In this section, a reduced-complexity algorithm attaining the optimal solution of GS-Lasso problems will be developed using the alternating direction method of multipliers (ADMoM). The disclosure show that (6) admits an equivalent reformulation that can be solved via ADMoM. Before doing so, the following lemmas are needed.

Lemma 1 Consider the following convex minimization problem in the variable yεR^N

$\begin{array}{cc}{y}^{*}={\mathrm{argmin}}_y\left[\frac{c}{2}y^Ty-y^Ta+\lambda {y}_2\right].& \left(8\right)\end{array}$

Albeit non differentiable, (8) admits a closed-form solution. Specifically, the global minimizer y* is given by the following soft-thresholding vector operation expressed in terms of [a]₊:=max{0,a} as

$\begin{array}{cc}{y}^{*}={\frac{a}{c{a}_2}\left[{a}_2-\lambda \right]}_{+}.& \left(9\right)\end{array}$

Proof. It will be argued that the solver of (8) takes the form y=za for some scalar z≧0. This is because among all y with the time l₂-norm, the Cauchy-Schwarz inequality implies that the maximizer of a^Ty is collinear with (and in the same direction of) a. Substituting y=za into (8) renders the problem scalar z≧0, with solution z⁸=(∥a∥−λ)₊/(c∥a∥₂), which completes the proof.

Lemma 2 Consider the following non-smooth convex problem in the vector variable yεR^N

$\begin{array}{cc}{y}^{*}={\mathrm{argmin}}_y\left[\frac{c}{2}y^Ty-y^Ta+\lambda {y}_1\right].& \left(10\right)\end{array}$

Using operator T_λ(·) defined as

T _λ(a):=[sgn(a₁)[|a₁|−λ]₊, . . . , sgn(a_N)[|a_N|−λ]₊]^T

the global minimized of (10) can be written as

$\begin{array}{cc}{y}^{*}=\frac{1}{c}T_{\lambda }\left(a\right).& \left(11\right)\end{array}$

If non-negativity of the entries in y is imposed, the solution of (10) subject to y≧0_N is obtained by using the vector soft-thresholding operator

$\begin{array}{cc}{T}_{\lambda }^{+}\left(a\right):={\left[\max \left\{0,a_1-\lambda \right\},\dots ,\max \left\{0,a_N-\lambda \right\}\right]}^T\mathrm{as}y^{*}=\frac{1}{c}T_{\lambda }^{+}\left(a\right).& \left(12\right)\end{array}$

Proof. Note first that (12) can be solved element-wise; specifically, each entry y_i of y is found by solving the scalar problem y*_i=arg min_yC(y):=(c/2y²−a_iy+λ|y|), which has a non-differentiable cost. The necessary and sufficient condition for y*_i to minimize C(y) is [27, p. 92]

$\begin{array}{cc}\{\begin{array}{cc}a_i\le \lambda ,& \mathrm{if}y_i^{*}=0\\ {\mathrm{cy}}_i^{*}-a+\lambda \frac{y_i^{*}}{y_i^{*}}=0,& \mathrm{if}y_i^{*}\ne 0\end{array}& \left(13\right)\end{array}$

which is satisfied by y*_i=sgn(a_i[|a_i|−λ]₊; see also [21]. When y is enforced to be non-negative, solution (12) can be easily derived from (13) element-wise. □

Consider now the N_bN_g×1 auxiliary vector variables γ and ξ, and neglect irrelevant terms to re-write the GS-Lasso problem as:

$\begin{array}{cc}\left\{\hat{\theta },\hat{\gamma },\hat{\xi }\right\}=\arg \underset{\theta ,\gamma ,\xi }{\min }\left[\frac{1}{2}R\theta -{\theta }^Tt+{\lambda }_1{\xi }_1+{\lambda }_2{}_G\left(\gamma \right)\right]s.t.\xi \succcurlyeq 0_{N_bN_g}\theta =\gamma ,\theta =\xi & \left(14\right)\end{array}$

where R:=Σ_r=1^NrB_xr^TB_xr, r:=Σ_r=1^NrB_xr^Tφ_xr, and φ_xr:={circumflex over (φ)}_xr−σ_v^L1_N. For simplicity, σ_v² is assumed to be known; however, it could be incorporated in (14) and estimated as the intercept.

Letting η and μ denote the Lagrange multipliers associated with the equality constraints Θ=γ and Θ=ξ, respectively, the quadratically augmented Lagrangian function of problem (14) is

$\begin{array}{cc}ℒ\left(\theta ,\gamma ,\xi ,\eta ,\mu \right)=\frac{1}{2}{\theta }^TR\theta -{\theta }^Tr+{\lambda }_1{\xi }_1+{\lambda }_2{}_G\left(\gamma \right)+{\eta }^T\left(\theta -\gamma \right)+{\mu }^T\left(\theta -\xi \right)+\frac{c_1}{2}{\theta -\gamma }_2^2+\frac{c_2}{2}{\theta -\xi }_2^2& \left(15\right)\end{array}$

where c₁, c₂>0 are arbitrary constants. Then, for any initial vectors γ⁽⁰⁾, ξ⁽⁰⁾, η⁽⁰⁾, μ⁽⁰⁾, the ADMoM algorithm entails the following primal-dual iterative updates

$\begin{array}{cc}{\theta }^{\left(j\right)}=\arg \underset{\theta }{\min }ℒ\left(\theta ,{\gamma }^{\left(j-1\right)},{\xi }^{\left(j-1\right)},{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(16a\right)\\ \left\{{\gamma }^{\left(j\right)},{\xi }^{\left(j\right)}\right\}=\arg \underset{\gamma ,\xi \succcurlyeq 0}{\min }ℒ\left({\theta }^{\left(j\right)},\gamma ,\xi ,{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(16b\right)\\ {\eta }^{\left(j\right)}={\eta }^{\left(j-1\right)}+c_1\left({\theta }^{\left(j\right)}-{\gamma }^{\left(j\right)}\right)& \left(16c\right)\\ {\mu }^{\left(j\right)}={\mu }^{\left(j-1\right)}+c_2\left({\theta }^{\left(j\right)}-{\xi }^{\left(j\right)}\right)& \left(16d\right)\end{array}$

where j=1, 2, . . . is the iteration index. The first step updates the primal vector Θ^(j) by using the values of the auxiliary variables and the Lagrange multipliers obtained at the previous iteration j−1; since L(·) is quadratic in Θ, the convex optimization problem (16a) can be solved in closed form as

θ^(j)=(R+(c₁+c₂)I_NDNg)⁻¹×(r+c₁γ^(j-1)+c₂ξ^(j-1)−η^(j-1)−μ^(j-1)).  (17)

Next, variables γ and ξ can be updated using the newly computed vector θ^(j), with the Lagrange multipliers fixed from the previous iteration. Inspection of the function L(θ^(j), γ, ξ, η^(j-1), μ^(j-1)) reveals that (16b) can be split into two sub-problems, where minimization over γ and ξ can be performed separately. After neglecting irrelevant terms, minimization of (16b) wrt γ reduces to the following non-differentiable convex problem

$\begin{array}{cc}{\gamma }^{\left(j\right)}=\arg \underset{\gamma }{\min }\left[\frac{c_1}{2}{\gamma }^T\gamma -{\gamma }^T\left(c_1{\theta }^{\left(j\right)}+{\eta }^{\left(j-1\right)}\right)+{\lambda }_2{}_G\left(\gamma \right)\right]& \left(18\right)\end{array}$

which, in turn, can be separated in following N_g sub-problems

$\begin{array}{cc}{\gamma }_g^{\left(j\right)}=\arg \underset{{\gamma }_g}{\min }\left[\frac{c_1}{2}{\gamma }_g^T{\gamma }_g-{\gamma }_g^T\left(c_1{\theta }_g^{\left(j\right)}+{\eta }_g^{\left(j-1\right)}\right)+{\lambda }_2{{\gamma }_g}_2\right],g=1,\dots ,N_g& \left(19\right)\end{array}$

where γ_g and η_g are N_b×1 sub-vectors of γ and η, respectively, collecting elements {γ_n}_n=Ng(g-1)+1^Ngg and {η_n}_n=Ng(g-1)+1^Ngg.

From Lemma 1, the global minimizer of each sub-problem (19) is given by

$\begin{array}{cc}{\gamma }_g^{\left(j\right)}=\left({\Theta }_g^{\left(j\right)}+c_1^{-1}{\eta }_g^{\left(j-1\right)}\right)\frac{{\left[{\mathrm{—c}}_1{\Theta }_g^{\left(j\right)}+{\eta }_g^{\left(j-1\right)}{—}_2-{\lambda }_2\right]}_{+}}{{\mathrm{—c}}_1{\Theta }_g^{\left(j\right)}+{\eta }_g^{\left(j-1\right)}g{—}_2}.& \left(20\right)\end{array}$

Upon neglecting constant terms, minimization of (16b) wrt the non-negative variable ξ can be obtained after solving the following non-smooth convex problem

$\begin{array}{cc}{\xi }^{\left(j\right)}={\mathrm{argmin}}_{\xi \succcurlyeq 0}\left[\frac{c_2}{2}T_{\xi -\xi }T\left({\mu }^{\left(j-1\right)}+c_2{\Theta }^{\left(j\right)}\right)+{\lambda }_1{\mathrm{—\xi —}}_1\right]& \left(21\right)\end{array}$

which, from Lemma 2, admits the following closed-form solution

$\begin{array}{cc}{\xi }^{\left(j\right)}=\frac{1}{c_2}T_{{\lambda }_1}^{+}\left(c_2{\Theta }^{\left(j\right)}+{\mu }^{\left(j-1\right)}\right).& \left(22\right)\end{array}$

The overall ADMoM-based solver for GS-Lasso problems is tabulated as Algorithm 1.

Algorithm 1 ADMoM-based GS-Lasso solver
Initialize γ(0) = 0NbNg, ξ(0) = 0NbNg, η(0) = 0NbNg, and μ(0) =
0NbNg
Form R and r
for j = 0, 1, . . . do
Update θ(j)  via (17)
Update γg(j)  via (20) for all g = 1, . . . , Ng
Update ξ(j)  via (22)
Update η(j) = η(j−1) + c1 (θ(j) − γ(j))
Update μ(j) = μ(j−1) + c2 (θ(j) − ξ(j))
end for

The distinct feature of the proposed ADMoM-based algorithm for solving the GS-Lasso problem (14) is its computationally affordable implementation, offered by the closed-form expressions for the primal variable updates; as well as the simple updates of the dual variables η^(j) and μ^(j). Furthermore, since ADMoM has provable convergence to the global minimizer when the considered problem is convex, convergence of the proposed algorithm to ̂Θ in (6) is ensured as stated next.

Proposition 1

For any c₁,c₂>0 and any initializing vectors γ⁽⁰⁾, ξ⁽⁰⁾, η⁽⁰⁾ and μ⁽⁰⁾, the iterates (17) for Θ^(j), ( ) for {γ_g^(j)}, ( ) for ξ^(j), and (16c)(16d) for η^(j) and μ^(j), respectively, are convergent. Also, Θ^(j) converges to the solution of the GS-Lasso (6); i.e.,

$\begin{array}{c}\lim \\ j\to +\infty \end{array}{\Theta }^{\left(j\right)}\widehat{=}\Theta .$

Shadow fading as well as possible slow temporal variations of the PU transmit-PSDs lead to time-varying {Φ_xr(f)}. Following the lines, time-varying PSDs can be tracked by employing the following time-weighted version of the GS-Lasso.

$\begin{array}{cc}\hat{\theta }\left(t\right)=\arg \underset{\theta \succcurlyeq {\theta }_{N_bN_g}}{\min }\left[\frac{1}{2}\sum _{r=1}^t{\beta }_{\tau ,t}\sum _{r=1}^{N_{\tau }}{{\ddot{\varphi }}_{x_r}\left(\tau \right)-B_{x_r}\left(\tau \right)\theta }_2^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)\right]& \left(23\right)\end{array}$

where β_τ,tε(0,1] is the so-called forgetting factor, and index τ=1, . . . , t emphasizes the temporal variability of channels and received PSDs. Also, to address the need for real-time processing, the estimation of 0 in (23) can be performed on-line, where each iteration of the ADMoM algorithm is performed after acquiring new estimates {φ_xr(τ)}. In this case, the ADMoM iteration index j coincides with the temporal index τ.

Algorithm 1 is centralized, meaning that the whole set of PSD estimates {φ_xr}_r=1^Nr are available at either an FC or a CR cluster head. To reduce the considerable message-passing overhead associated with globally sharing PSD measurements across CRs, and to address scalability and robustness concerns (FC constitutes an isolated point of failure), a distributed counterpart of Algorithm 1 can be derived along the lines.

Techniques for Spectrum Sensing under Channel Uncertainties are now described, in including a Group Sparse Total Least-Squares (TLS) technique.

Uncertainty in the matrices {B_xr} is manifested because of (i) errors in the estimates of {g_xg→xr} (with or without accounting for shadowing; (ii) grid offsets when PUs are located between grid points; and, (iii) basis expansion approximation errors. To cope with these perturbations, a robust version of the GS-Lasso is developed in this section. That is, TLS may be applied as described herein to address situations having perturbed compressive sampling, which may arise whenever any of the uncertainties result in a mismatch between actual and estimated bases reference functions computed for the model.

TLS is the workhorse used for estimating non-sparse vectors obeying an over-determined linear system of equations with uncertainty present in both the regression matrix and the observations (fully-perturbed model). Sparsity in the estimate was taken into account, where the TLS framework was extended to solve sparse under-determined fully-perturbed linear systems. The sparse TLS approach is broadened here to account for sparsity present both at individual entries, and also at groups of entries.

Define

$B:={\left[B_{x_1}^T,\dots ,B_{x_{N_r}}^T\right]}^T,\mathrm{and}\mathrm{let}\hspace{0.17em}{\hspace{0.17em}}^{-}\phi :={\left[{\phi }_{x_1}^T,\dots ,\Phi _{x_{N_r}}^T\right]}^T,$

and E a NN_r×N_gN_b matrix capturing perturbations corrupting the matrix B. Consider now estimating Θ as follows:

$\begin{array}{cc}\left\{\hat{\theta },\hat{E}\right\}=\arg \underset{\underset{\theta \succcurlyeq 0_{N_bN_g}}{E}}{\min }\left[\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+E\right)\theta }_2^2+\frac{1}{2}{E}_F^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)\right].& \left(24\right)\end{array}$

Relative to the classical TLS, the cost in the group sparse (GS-)TLS problem (24) is augmented with the regularization terms accounting for the two forms of sparsity inherent to Θ. Compared to, problem (24) includes also the term λ₂R_G(Θ).

Problem (24) is generally non-convex due to the presence of the product EΘ; thus, it is in general difficult to obtain a globally-optimal solution. However, a novel reduced-complexity algorithm with provable convergence to a stationary point of (24) will be developed in the ensuing section.

An Alternating Descent Solver is no described. The cost in (24) will be optimized here iteratively using a block coordinate descent algorithm, which cyclically minimizes it wrt E (keeping Θ fixed), and wrt Θ after fixing E. Specifically, the following two steps are performed at the i-th iteration:

$\begin{array}{cc}\left(i1\right)\mathrm{Fix}E={\hat{E}}^{\left(i-1\right)},\mathrm{and}\mathrm{update}{\hat{\theta }}^{\left(i\right)}\mathrm{as}{\hat{\theta }}^{\left(i\right)}=\arg \underset{\theta \succcurlyeq 0_{N_bN_g}}{\min }\left[\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+{\hat{E}}^{\left(i-1\right)}\right)\theta }_2^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)\right].& \left(25\right)\\ \left(i2\right)\mathrm{Fix}\theta ={\hat{\theta }}^{\left(i\right)},\mathrm{and}\mathrm{obtain}{\hat{E}}^{\left(i\right)}\mathrm{as}{\hat{E}}^{\left(i\right)}=\arg \underset{E}{\min }\frac{1}{2}{\stackrel{\_}{\varphi }B{\hat{\theta }}^{\left(i\right)}-E{\hat{\theta }}^{\left(i\right)}}_2^2+\frac{1}{2}{E}_F^2.& \left(26\right)\end{array}$

By fixing ̂E⁽ⁱ⁾, (24) boils down to a GS-Lasso problem; thus, ̂Θ⁽ⁱ⁾ can be computed by using the ADMoM-based solver after replacing (17) with the following update, where j still represents the index for the (inner) ADMoM iterations:

$\begin{array}{cc}{\hat{\theta }}^{\left(i,j\right)}={\left[{\left(B+{\hat{E}}^{\left(i-1\right)}\right)}^T\left(B+{\hat{E}}^{\left(i-1\right)}\right)+\left(c_1+c_2\right)l_{N_bN_g}\right]}^{-1}\times \hspace{1em}\left[{\left(B+{\hat{E}}^{\left(i-1\right)}\right)}^T\stackrel{\_}{\varphi }+c_1{\gamma }^{\left(j-1\right)}+c_2{\xi }^{\left(j-1\right)}-{\eta }^{\left(j-1\right)}-{\mu }^{\left(j-1\right)}\right].& \left(27\right)\end{array}$

The quadratic convex problem (26) admits the following closed form solution

Ê ⁽ⁱ⁾=(1+∥{circumflex over (θ)}⁽ⁱ⁾∥₂²)⁻¹[ φ−B{circumflex over (θ)}⁽ⁱ⁾]{circumflex over (θ)}^(i)T  (28)

which can be obtained after equating the derivative of the cost in (26) with zero. The overall solver for GS-TLS is tabulated as Algorithm 2.

Algorithm 2 GS-TLS
Initialize Ê(0) = 0NbNg×NbNg
while Stopping criterion is not satisfied (i iteration index) do
Initialize γ(0) = 0NbNg, ξ(0) = 0NbNg, η(0) = 0NbNg, μ(0) =
0NbNg
While Stopping criterion is not satisfied (j iteration index) do
Update θ(i,j) via (27)
Update γgj via (20) with θ(i,j) in place of θ(j), for all g =
1, . . . , Ng
Update ξ(j) via 22) with θ(i,j) in place of θ(j)
Update η(j) = η(j−1) + c1 (θ(i,j) − γ(j))
Update μ(j) = μ(j−1) + c2 (θ(i,j) − ξ(j))
end while
Update {circumflex over (θ)}(i) = θ(i,j)
UpdateÊ(i) = (1+ ∥{circumflex over (θ)}(i)∥22)−1[ φ − B{circumflex over (θ)}(i)] {circumflex over (θ)}(i)T
end while

Under certain conditions, the block coordinate descent algorithm is known to converge (at least) to a local optimum point, as asserted next.

Algorithm 3 Robust GS-TLS
Initialize ô(0)  = 0NbNg and Ê(0)  = 0NbNg×NbNg
while stopping criterion is not satisfied (i iteration index) do
Initialize γ(0) = 0NbNg , ξ(0) = 0NbNg , η(0) = 0NbNg , and μ(0) =
0NbNg
while Stopping criterion is not satisfied (j iteration index ) do
Update θ(i,j) via (37)
Update γg(j) via (20) with θ(i,j) in place of θ(j) for all g =
1, . . . , Ng
Update o(i,j) via (38)
Update ξ(j)  via (22) with θ(i,j)  in place of θ(i)
Update η(j) = η(j−1) + c1 (θ(i,j) − γ(j))
Update μ(j) = μ(j−1) + c2 (θ(i,j) − ξ(j))
end while
Update {circumflex over (θ)}(i) = θ(i,j)  and ô(i) = o(i,j)
Update Ê(i) = (1+ ∥{circumflex over (θ)}(i)∥22)−1[ φ − B{circumflex over (θ)}(i)+ô(i)] {circumflex over (θ)}(i)T
end while

Proposition 2. For any initialization {{circumflex over (θ)}⁽⁰⁾, Ê⁽⁰⁾}, the Iterates {{circumflex over (θ)}⁽ⁱ⁾, Ê⁽ⁱ⁾} (25)-(26) converge monotonically to a stationary point of problem (4).

The following describes outlier-resilient spectrum sensing techniques, including Robust GS-TLS.

The problem dealt with in the previous section accounts for uncertainty in the entries of the regression matrix B. However, due to particularly abrupt local shadow fading, failures of the sensing modules, or unexpected narrow-band impulsive noise and/or interference, CRs observations may be affected by abundant errors. This section develops schemes for discerning and removing the observations that largely deviate from the underlying model (a.k.a. outliers).

A simple heuristic to detect unreliable data could be to estimate Θ via (24), compute the residuals and, then, reject the PSD observations whose residuals exceed a certain threshold. A systematic method that accounts for possible outliers can be found, where the underlying linear regression model is augmented by an auxiliary outlier vector. Using this model, the receiver PSD at the CR locations φ can be expressed as

φ=(B+E)Θ+o+e  (29)

where the nonzero entries of the NN_r×1 real vector o capture outliers; and e is a proper vectorization of the periodogram estimation errors. Since few outliers are expected compared to the total number of data collected in the vector o is sparse.

Capitalizing on the three forms of sparsity emerging from (i) the grid-based model (group sparsity), (ii) the PSD basis expansion (single-coefficient sparsity) and (iii) the outliers (single-coefficient sparsity), and accounting for perturbations in the regression matrix, the following robust GS-TLS is considered. In this way, a robust GS-TLS technique is described that addresses errors introduced by uncertainties, such as estimate errors, grid offsets, and BEM approximation errors, and also handles measurement corruption by inclusion of outliers.

$\begin{array}{cc}\left[\hat{\theta },\hat{E},\hat{o}\right]=\arg \underset{\underset{\underset{\theta \ge 0_{N_bN_g}}{E}}{o}}{\min }\hspace{1em}\hspace{1em}[\hspace{1em}\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+E\right)\theta +o}_2^2+\frac{1}{2}{E}_F^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)+{\lambda }_3{o}_1]& \left(30\right)\end{array}$

where λ₃≧0 promotes the (single-coefficient) sparsity of vector ô. The nature of the perturbations captured in E and o is in general different: o collects unmodeled errors, whereas E describes (small) perturbations. It is also worth noticing that the support of ô reveals the unreliable CR data.

An alternating descent algorithm is described. Although (30) is generally a non-convex problem, a block coordinate descent algorithm can still be employed. In this case, the cost in (30) will be iteratively minimized wrt E and {Θ,o}; that is, the following two updates are performed at the i-th iteration:

$\begin{array}{cc}\left(i1\right)\mathrm{Fix}E={\hat{E}}^{\left(i-1\right)}\mathrm{and}\mathrm{solve}& \\ \left\{{\hat{\theta }}^{\left(i\right)},{\hat{o}}^{\left(i\right)}\right\}=\arg \underset{\underset{\theta \succcurlyeq 0_{N_bN_g}}{o}}{\min }\left[\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+{\hat{E}}^{\left(i-1\right)}\right)\theta +o}_2^2+{\lambda }_1{\theta }_1+{\lambda }_2{}_G\left(\theta \right)+{\lambda }_3{o}_1\right].& \left(31\right)\\ \left(i2\right)\mathrm{Fix}\theta ={\hat{\theta }}^{\left(i\right)}\mathrm{and}o={\hat{o}}^{\left(i\right)}\mathrm{and}\mathrm{update}{\hat{E}}^{\left(i\right)}\mathrm{as}& \\ {\hat{E}}^{\left(i\right)}=\arg \underset{E}{\min }\frac{1}{2}{\stackrel{\_}{\varphi }-B{\hat{\theta }}^{\left(i\right)}+E{\hat{\theta }}^{\left(i\right)}+{\hat{o}}^{\left(i\right)}}_2^2+\frac{1}{2}{E}_F^2.& \left(32\right)\end{array}$

The quadratic problem (32) can be solved in closed form, to obtain

Ê ⁽ⁱ⁾=(1+∥{circumflex over (θ)}⁽ⁱ⁾∥₂²)⁻¹[ φ−B{circumflex over (θ)}⁽ⁱ⁾+ô⁽ⁱ⁾]{circumflex over (θ)}^iT.  (33)

As for (31), the ADMoM can be employed to find its optimal solution. To this end, (31) can be re-formulated as:

$\begin{array}{cc}\left\{{\hat{\theta }}^{\left(i\right)},\gamma ,\xi ,{\hat{o}}^{\left(i\right)}\right\}=\arg \underset{\gamma ,\xi ,o}{\underset{0\succcurlyeq 0_{N_bN_g}}{\min }}\left[\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+{\hat{E}}^{\left(i-1\right)}\right)\theta +o}_2^2+{\lambda }_1{\xi }_1+{\lambda }_2{}_G\left(\gamma \right)+{\lambda }_3{o}_1\right]s.t.\xi \succcurlyeq 0_{N_bN_g}\theta =\gamma ,\theta =\xi & \left(34\right)\end{array}$

with γ and ξ denoting auxiliary vector variables. Letting again η and μ denote the Lagrange multipliers associated with the constraints Θ=γ and Θ=ξ, respectively, the quadratically augmented Lagrangian function (34) is given by

$\begin{array}{cc}ℒ\left(\theta ,\gamma ,\xi ,o,\eta ,\mu \right)=\frac{1}{2}{\stackrel{\_}{\varphi }-\left(B+{\hat{E}}^{\left(i-1\right)}\right)\theta +o}_2^2+{\lambda }_1{\xi }_1+{\lambda }_2{}_G\left(\gamma \right)+{\lambda }_3{o}_1+{\eta }^T\left(\theta -\gamma \right)+{\mu }^T\left(\theta -\xi \right)+\frac{c_1}{2}{\theta -\gamma }_2^2+\frac{c_2}{2}{\theta -\xi }_2^2.& \left(35\right)\end{array}$

Starting from any initial vectors o⁽⁰⁾, γ⁽⁰⁾, ξ⁽⁰⁾, η⁽⁰⁾, μ⁽⁰⁾, each iteration j of the ADMoM (within each iteration i of the block coordinate descent) proceeds in these steps:

$\begin{array}{cc}{\theta }^{\left(i,j\right)}=\arg \underset{\theta }{\min }ℒ\left(\theta ,{\gamma }^{\left(j-1\right)},{\xi }^{\left(j-1\right)},o^{\left(i,j-1\right)},{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(36a\right)\\ {\gamma }^{\left(i,j\right)}=\arg \underset{\gamma }{\min }ℒ\left({\theta }^{\left(i,j\right)},\gamma ,{\xi }^{\left(j-1\right)},o^{\left(i,j-1\right)},{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(36b\right)\\ {\xi }^{\left(j\right)}=\arg \underset{\xi \succcurlyeq 0}{\min }ℒ\left({\theta }^{\left(i,j\right)},{\gamma }^{\left(j\right)},\xi ,o^{\left(i,j-1\right)},{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(36c\right)\\ o^{\left(i,j\right)}=\arg \underset{o}{\min }ℒ\left({\theta }^{\left(i,j\right)},\gamma ,{\xi }^{\left(j-1\right)},o,{\eta }^{\left(j-1\right)},{\mu }^{\left(j-1\right)}\right)& \left(36d\right)\\ {\eta }^{\left(j\right)}={\eta }^{\left(j-1\right)}+c_1\left({\theta }^{\left(j\right)}-{\gamma }^{\left(j\right)}\right)& \left(36e\right)\\ {\mu }^{\left(j\right)}={\mu }^{\left(j-1\right)}+c_2\left({\theta }^{\left(j\right)}-{\xi }^{\left(j\right)}\right).& \left(36f\right)\end{array}$

Since (34) is convex and satisfies the requirements for the ADMoM to be convergent, iterates {Θ^(i,j),o^(i,j)} will converge to the solution {̂Θ⁽ⁱ⁾,̂o⁽ⁱ⁾} of (31).

Problem (36a) admits the closed-form solution given by

$\begin{array}{cc}{\theta }^{\left(i,j\right)}={\left[{\left(B+{\hat{E}}^{\left(i-1\right)}\right)}^T\left(B+{\hat{E}}^{\left(i-1\right)}\right)+\left(c_1+c_2\right)I_{N_bN_g}\right]}^{-1}\times \hspace{1em}[\hspace{1em}\left(B+{\hat{E}}^{\left(i-1\right)}\right)\left(\stackrel{\_}{\varphi }+o^{\left(i,j-1\right)}\right)+c_1{\gamma }^{\left(j-1\right)}+c_2{\xi }^{\left(j-1\right)}-{\eta }^{\left(j-1\right)}-{\mu }^{\left(j-1\right)}]& \left(37\right)\end{array}$

where γ^(j) and ξ^(j) are still computed via (20) and (22). Using Lemma 2, it is possible to show that the solution of (36d) is computed via soft-thresholding as

o ^(i,j)=ℑ_λ3( φ−(B+Ê^(i-1))θ^(i,j)).  (38)

At each step of the ADMoM algorithm the soft-thresholding in (38) tags as outliers the (current) residuals φ−(B+Ê^(i-1))θ^(i,j) that exceed λ₃. The overall solver is tabulated as Algorithm 3.With arguments similar to Proposition 2, the following result can be asserted.Proposition 3. For any initialization {{circumflex over (θ)}⁽⁰⁾, ô⁽⁰⁾, Ê⁽⁰⁾)}, the iterates {{circumflex over (θ)}⁽ⁱ⁾, ô⁽ⁱ⁾, Ê⁽ⁱ⁾} in (32)-(31) converge monotonically to a stationary point of problem (30).

Simulations are now described. Consider a set of N_r=50 CRs uniformly distributed in an area of 100 m×100 m, cooperating to localize N_s=2 active PUs and estimate their PSD map. CRs and PU sources are marked with blue circles and red triangles, respectively, in FIG. 4. PUs transmit raised cosine pulses with unitary amplitude (0 dB), roll-off factor =0.5, and bandwidth W=10 MHz. They share the band B=[100,200] MHz with spectra centered at frequencies f_c=115 and 175 MHz for “PU 1” and “PU 2”, respectively. CRs adopt a path loss-only model to accomplish the sensing task.

Transmitted signals are searched over a grid of N_b=10 evenly spaced center frequencies f_c=95+vW, vε{1, . . . , 10}. Each CR computes periodogram samples at N=64 frequencies at signal-to-noise-ratio (SNR)−5 dB, and averages them across τ=100 time-slots to form ̂Φ_xr(τ,f_k), k=1, . . . , 64, as in (4).

In the first experiment, the PSD generated by PU s experiences only small-scale fading in its propagation from x_s to any location x, where it is measured in the presence of noise with variance σ_v²=0.1. To simulate small-scale fading {h_xs→x}, a 6-tap Rayleigh model with exponential power delay profile is adopted. Since the expected gain adheres to a path loss propagation law, the regression matrix is in this case perfectly known. FIG. 5A depicts the true PSD map summed across frequencies, which peaks at the active PU locations. To localize and estimate the transmit-PSDs, a grid of N_g=100 equidistant points is used. The map obtained by using the GS-Lasso-based sensing algorithm is shown in FIG. 5B, which also depicts the estimated positions of the transmitting PUs along with their transmission powers (in dB) represented by the coefficients of the normalized raised-cosines. One can readily notice that “PU 1” is perfectly localized and its transmit-PSD is estimated accurately. As for “PU 2”, its location is revealed although spurious power is also leaked to an adjacent grid point. The sparsity-promoting parameters λ₁ and λ₂ are set to λ₁=30·max{B^T̂φ} and λ₂=10·max_r{|B_xr^T̂φ|₂}, respectively. FIG. 6 corroborates the convergence of the GS-Lasso solver by showing the evolution of the elements of Θ^(j).

In FIG. 7, the transmit-PSDs undergo not only small-scale fading but also log-normal shadowing. As the CRs employ a path loss-based model, shadowing here perturbs the regression matrix. Expressed in logarithmic scale, the shadowing process has zero mean and standard deviation 6 dB. The estimated PSD maps obtained by using the “plain” GS-Lasso and the GS-TLS-based sensing algorithm are compared in FIG. 8.

FIG. 8A illustrates that the GS-Lasso is unable to localize the two PUs, as clouds of PU sources are falsely revealed around the actual locations of “PU 1” and “PU 2”. Also, it does not accurately estimate their PSDs. On the other hand, the GS-TLS algorithm reveals the exact location of both PUs, although a rather small amount of spurious power is leaked to a grid point close to “PU 1”. Note also that the transmit-powers are estimated with considerably higher accuracy. Numerical experiments have shown that only a few (5 to 10) iterations suffice for the alternating descent algorithm to converge.

The enhanced localization and power estimation capabilities impact also the subsequent CR power allocation task, which relies on the estimated coverage region of the PU-transmitters to re-use the licensed bands without causing harmful interference to any potential PU receiver. GS-Lasso will be preferable if a coarse description of the “interference-heavy” areas in terms of PU activity is desired over say accurate localization and transmit-power estimation of the PUs. Such a coarse can be useful for e.g., temporal (rather than spatio-temporal) frequency re-use purposes.

As described, abrupt local shadow fading may severely compromise the PSD estimates at CRs, and thus degrade the sensing performance. This is the case considered in FIG. 9, where an obstacle positioned in the upper-left part of the monitored area causes deep fades of the receive-power at some CRs. FIG. 10A demonstrates that “plain” GS-Lasso fails to localize the two PU transmitters. Activity is revealed around the actual location of “PU 1”, and shadowing causes the false detection of a third low-power PU in position x_g=(80,70) transmitting over the same band of “PU 2”. This false-detection event is not present in FIG. 10B, where the robust GS-TLS algorithm is used. In fact, “PU 2” is well-localized and its transmit-PSD is estimated accurately. A small amount of power is still dribbled on an adjacent grid point of “PU 2”. With λ₃=30·max{B^T̂φ}, further analysis of the data reveals that 15% of the periodogram samples, specifically those collected by the “faded CRs”, was declared unreliable, and was thus discarded

Cognitive Radio Spectrum Prediction Using Dictionary Learning

As described above, the radio frequency (RF) spectrum is a precious resource that must be utilized efficiently. Fixed spectrum allocation, which confers exclusive access rights on spectrum license holders, has resulted in significant under-utilization of the valuable spectral resource, depending on time and locations. The cognitive radio (CR) strategy aims at alleviating this inefficiency by allowing unlicensed secondary users to opportunistically transmit, provided that the transmissions do not disturb the communication of licensed primary users (PUs). To achieve the necessary protection of PU systems, CR systems provide spectrum sensing and intelligent resource allocation. Spectrum sensing is used to identify unused spectral resources in the frequency, time and space domains. The “spectrum holes” can then be exploited through agile resource allocation.

A simplifying assumption often made for spectrum sensing is that the spectrum occupancy is more or less invariant over the deployment region of the CR systems. Based on this, spectrum sensing is often performed in a collaborative fashion, where the band occupancy by a common set of PU transmitters is detected using observations fused from multiple CRs. This mitigates effectively fading and shadowing, which impede reliable detection of PU presence.

However, the assumption might not hold when PU systems employ a small RF footprint for significant spatial reuse, or the CR network grows in size and gets deployed in a broader geographical region. An instrumental concept in this case is the RF cartography, which provides a map of RF power distribution over space and reveals the spatial variation of spectrum occupancy. Such a construct is useful for optimizing CR network operations, not only in the PHY/MAC, but also in higher layers.

Techniques are described herein that allow interference power present at each CR node in the network to be acquired, and also to predict its future levels. Some example techniques may address the challenge that CRs do not have prior information on the number of PU emitters and the corresponding PU-CR channel gains, which may be needed for combining the measurements from different sensors. Moreover, CRs might not be able to report their measurements every time, due to energy-saving sleep modes or congested signaling channels. A network controller operating accordance to the techniques described herein accounts for missing observations by performing appropriate interpolation. The future spectrum state must be inferred based on past measurements. Machine learning and compressive sensing techniques are employed herein to tackle these challenges. A general dictionary learning framework is adapted here to learn the spatial and temporal patterns of the RF power distribution. To aid in spatial interpolation, the topology information of the CR network is exploited in a semi-supervised learning fashion. Batch and online algorithms are developed. The online techniques can track the slow variation of the RF power distribution, and features low computational complexity and lax memory requirement.

Spatio-temporal spectrum prediction algorithms for cognitive radios (CRs) are developed using the framework of dictionary learning and compressive sensing. The interference power levels at each CR node locations are predicted using the measurements from a subset of CR nodes without a priori knowledge on the primary transmitters. A semi-supervised dictionary learning approach for the interpolation of missing observations in space is described. The proposed techniques are also extended for spatio-temporal prediction. Batch and online alternatives are presented, where the online algorithm features low complexity and memory requirements. Numerical tests verify the performance of the proposed novel methods.

Returning to FIG. 1, consider the CR network (system 10) consisting of M nodes (sensors 12), deployed in a geographical area, over which the interference due to the incumbent PU systems vary, albeit smoothly. CRs 12 form a mesh network by identifying their neighbors, and cooperate for spectrum sensing. System 10 acquires the interference level at each CR node location, based on the measurements collected from a subset of the CR nodes per time. The missing measurements may be due to various practical limitations, such as errors and congestion in the control channel, or the fact that radios are in the sleep mode to save battery.

Suppose that there are K PU transmitters (RF-enabled devices 15) in the area with the k-th PU transmitting at power p_k(t). Let g_mk(t) denote the channel gain from the k-th PU to the m-th CR. Then, the interference power level π_m(t) perceived at the m-th CR can be modeled as

$\begin{array}{cc}{\pi }_m\left(t\right)=\sum _{k=1}^Kg_{\mathrm{mk}}\left(t\right)\mathrm{Pk}\left(t\right),m\in ℳ:=\left\{1,2,\dots ,M\right\}& \left(1.1\right)\end{array}$

Upon defining vectors π(t):=[π₁(t), . . . , π_M(t)]^τ (.^τ denotes transposition) and P(t):=[p₁(t), . . . , p_k(t)]^τ as well as matrix G(t) whose (m,k)-entry is g_mk(t), the matrix-vector counterpart of (1) can be expressed as

π(t)=G(t)P(t).  (1.2)

At each time t, a subset M^obs(t)⊂ of CRs observe the interference power levels. The measurements from these CRs can be stacked in vector y^obs(t)εR^|Mobs(t)| given as

y ^obs(t)=0(t)π(t)+z(t)  (1.3)

where (t)εR^|Mobs(t)| is the measurement noise vector, and (t) is a matrix consisting of the m-th row of an M×M identity matrix, where mεM^obs(t).

The problem of estimating G(t) and p(t), given the past and the current measurements y^obs(τ) for τ=1, 2, . . . , t, when there were no missing observations was tackled using dictionary learning in [13]. Here, the goal is to predict the missing interference levels π_m(t) for π_m(t) for mεM^miss(t):=M\M^obs(t), given {y^obs(τ)}_τ=1^t.

Techniques for spatial spectrum prediction using semi-supervised dictionary learning are described. Prompted by (1.2), suppose that π can be represented as a linear combination of a small number of bases (atoms) taken from a dictionary. Let εR^M×Q denote a dictionary with Q atoms. Then, the preceding assumption amounts to

π=Ds  (1.4)

where vector εR^Q is sparse. Fourier bases or the wavelet bases are some of the dictionaries often used for a variety of natural or man-made signals. Based on this model, the techniques described below leverage compressive sensing and machine learning to predict the unobserved interference levels in space. In the following, a two-stage approach is first considered, in which the dictionary learned in the training phase is used for the desired prediction task in the operational phase. An alternative algorithm capable of performing the dictionary learning and spatial prediction simultaneously will be presented subsequently.

An example two-phase batch implementation is described. Instead of using off-the-shelf bases such as Fourier or the wavelet bases, D can be directly learned from training data, which can be collected through a sounding procedure. In the presence of missing entries in the data, it is helpful to augment this learning process with additional structural information. In this work, the network topology information, which is typically maintained for various network control tasks such as routing, is leveraged in the framework of semi-supervised learning.

Let Aε{1,0}^M×M denote the adjacency matrix of the CR network topology. Thus, the (m,m′)-th entry a_m,m′ of for m,m′ε is 1 if nodes m and m′ are neighbors, and 0 otherwise. Then, the Laplacian matrix is defined as :=diag(1)-, where 1 is the all-one vector, and diag(v) is a diagonal matrix with the entries of vector v on its diagonal.

Training Phase:

In the training phase, given a training set {y_n^ovs}_n=1^N, which may contain missing entries, the goal is to form an estimate {circumflex over (D)} of D such that y_n^obs≈O_n{circumflex over (D)}s_n for n=1, 2, . . . , N, where coefficients s_n are sparse, and matrix O_n discards the missing entries. Specifically, the following optimization problem is solved to obtain {circumflex over (D)}.

$\begin{array}{cc}\hat{D}:=\arg {\min }_{D\in ,\left\{s_n\right\}}\sum _{n=1}^Nf_n\left(s_n,D\right)\mathrm{where}& \left(1.5\right)\\ f_n\left(s,D\right):=\frac{1}{2}{y_n^{\mathrm{obs}}-O_n\mathrm{Ds}}_2^2+{\lambda }_s{s}_1+\frac{1}{2}{\lambda }_Ls^{\tau }D^{\tau }\mathrm{LDs}& \left(1.6\right)\\ :=\left\{\left[d_1,\dots ,d_Q\right]\epsilon {ℝ}^{M\times Q}:{d_q}_2^2\le 1,q=1,\dots Q\right\}& \left(1.7\right)\end{array}$

Here, the first term in (1.6) promotes fitness of the reconstruction to the training datum in a least-squares (LS) sense, and the l₁-norm-based regularization term encourages sparsity in s with λ_s>0 playing the role of a tuning parameter. The third term in (1.6) can be re-written as

s ^τ D ^τ LDs=Σ _m=1 ^MΣ_m′=1^Ma_m,m,(π_m−π_m′)²  (1.8)

indicating that it encourages the interference levels experienced at neighboring nodes to be similar, with λ_L>0 being a tuning parameter.

To appreciate the role of this Laplacian matrix-based regularization, suppose that a CR never reports an observation during the entire training period. Then, without the last term, the corresponding row in cannot be estimated, making it impossible to predict the interference level at this CR's location. The presence of the Laplacian term allows one to estimate the missing entry relying on neighbors' measurements.

Operational Phase: Once 15 has been obtained as in (5), the operational phase predicts the interference levels. First, a sparse coding step is performed at each time t to estimate the sparse coefficient (t) corresponding to the measurement y^obs(t); that is,

$\begin{array}{cc}\hat{s}\left(t\right):=\arg \underset{s}{\min }\frac{1}{2}{y^{\mathrm{obs}}\left(t\right)-O\left(t\right)\hat{D}s}_2^2+{\lambda }_2{s}_1+\frac{1}{2}{\lambda }_Ls^{\tau }{\hat{D}}^{\tau }L\hat{D}s& \left(1.9\right)\end{array}$

Then, the desired interference levels in π(t) that include the missing entries can be recovered by {circumflex over (π)}(t):{circumflex over ({circumflex over (=)}(t):={circumflex over (D)}ŝ(t).

Implementation:

Problem (1.9) is convex and there are various specialized algorithms available for solving the problems of this sort extremely fast. On the other hand, (5) is nonconvex, and it is difficult to obtain globally optimal solutions. However, the problem is convex with respect to and {s_n} individually. Thus, to find a locally optimal solution, a block-coordinate descent (BCD) algorithm can be employed until convergence.

Specifically, at the k-th iteration, updates are done as:

$\begin{array}{cc}\left\{{\hat{s}}_n^{\left(k\right)}\right\}:\arg \underset{\left\{s_n\right\}}{\min }\sum _{n=1}^Nf_n\left(s_n,{\hat{D}}^{\left(k-1\right)}\right)& \left(1.10\right)\\ {\hat{D}}^{\left(k\right)}:=\arg \underset{D\epsilon }{\min }\sum _{n=1}^Nf_n\left({\hat{s}}_n^{\left(k\right)},D\right)& \left(1.11\right)\end{array}$

where {circumflex over (D)}^(k) and {ŝ_n^(k)} are the k-th iterates. Note that (1.10) can be solved separately per n=1, 2, . . . , N using the same solver as the one for (1.9). To solve (1.11), a BCD algorithm can be once again employed over the columns of D. Define {tilde over (L)}_n:=O_n^τO_n+λ_LL, and let s_n,q denote the q-th entry of vector s_n. Then, the overall dictionary training algorithm is presented in Table 1.

TABLE 1
Input: training set {ynobs}n=1N, {On}, initial dictionary D0, λs, L and λL
Output: {circumflex over (D)} := [{circumflex over (d)}1, {circumflex over (d)}2, . . . , {circumflex over (d)}Q]
1:  Set D = D0
2:  Repeat
    Perform sparse coding with fixed {circumflex over (D)}.
3:  For n = 1, 2, . . . , N
4:  ŝn = arg mins fn (s, {circumflex over (D)})
5.  Next n
    Perform dictionary update with fixed {ŝn}
6.  Repeat
7.  For q = 1, 2, . . . , Q
8:  dq ={circumflex over (d)}q + (Σn =1N  {tilde over (L)}nsn,q2)−1·
    [Σn =1N sn,q (OnTynobs − {tilde over (L)}n{circumflex over (D)}sn)]
9:  {circumflex over (d)}q = dq/ max{∥ dq∥2,1}
10: Next q
11: Until convergence
12: Until convergence

Example online techniques are now described. In order to track time-varying statistics of the interference patterns, an online algorithm can be derived, in which the dictionary training and spatial interference prediction are performed jointly at the same time. Compared to the batch training discussed, the online algorithm can perform the computation recursively, resulting in significant savings in complexity and memory.

Specifically, the following formulation is adopted, which weights recent observations more heavily.

$\begin{array}{cc}\hat{D}\left(t\right),\left\{\hat{s}\left(t\right)\right\}=\arg \underset{D\epsilon ,\left\{s\left(t\right)\right\}}{\min }\sum _{\tau =1}^t{\beta }^{t-\tau }\left(\frac{1}{2}{y^{\mathrm{obs}}\left(\tau \right)-O\left(\tau \right)\mathrm{Ds}\left(\tau \right)}_2^2+{\lambda }_ss\left(\tau \right)1+\frac{1}{2}{\lambda }_Ls^{\tau }\left(\tau \right)D^{\tau }\mathrm{LDs}\left(\tau \right)\right)& \left(1.12\right)\end{array}$

where βε(0,1] is a forgetting factor. Instead of solving problem (1.12) in a batch fashion for the entire time horizon τ=1, 2, . . . , t whenever a new observation y^obss (t) arrives at each time t, an online approach updates only the “current” coefficient vector ŝ(t), while the past ones ŝ(t−1), . . . , ŝ(1) are held fixed. Nevertheless, it can be shown under mild conditions that {circumflex over (D)} so obtained converges as t→∞ to the same {circumflex over (D)} as would be obtained from a batch approach.

Although the dictionary update depends on the entire observation history, a recursive computation can avoid storing the past observations and calculations. For this, it is useful to maintain the following quantities:

$\begin{array}{cc}A\left(t\right):=\sum _{r=1}^t{\beta }^{t-\tau }\hat{s}\left(t\right){\hat{s}}^{\tau }\left(\tau \right)=\beta A\left(t-1\right)+\hat{s}\left(t\right){\hat{s}}^{\tau }\left(t\right)& \left(1.13\right)\\ \begin{array}{c}{A}_m\left(t\right):=\sum _{\tau =1}^t{\beta }^{t-\tau }{}_{\left\{m\in {ℳ}^{\mathrm{obs}}\left(t\right)\right\}}\hat{s}\left(t\right){\hat{s}}^{\tau }\left(t\right)\\ =\beta A_m\left(t-1\right)+{}_{\left\{m\in {ℳ}^{\mathrm{obs}}\left(t\right)\right\}}\hat{s}\left(t\right){\hat{s}}^{\tau }\left(t\right),m\in ℳ\end{array}& \left(1.14\right)\\ \begin{array}{c}B\left(t\right):=\sum _{\tau =1}^t{\beta }^{t-\tau }O^{\tau }\left(\tau \right)y^{\mathrm{obs}}\left(\tau \right){\hat{s}}^{\tau }\left(\tau \right)\\ =\beta B\left(t-1\right)+O^{\tau }\left(t\right)y^{\mathrm{obs}}\left(t\right){\hat{s}}^{\tau }\left(t\right)\end{array}& \left(1.15\right)\end{array}$

where ∥_{·} is an indicator function equal to 1 if the condition inside the braces are satisfied, and 0 otherwise.

Similar to the batch case described above, the dictionary update amounts to solving (12) for {circumflex over (D)}(t) with {ŝ(τ)}_τ=1^t fixed. Let ŝ_j(τ) denote the j-th entry of ŝ(τ), and A_m,jq(t) and A_jq(t) the (j,q)-th entry of matrices A_m(t) and A(t), respectively. Also, let b_j(t) represent the j-th column of B(t). Then, upon defining

$\begin{array}{cc}{\Phi }_{j,q}\left(t\right):=\sum _{\tau =1}^T{\beta }^{t-\tau }{\hat{s}}_j\left(\tau \right){\hat{s}}_q\left(\tau \right)\left(O^{\tau }\left(\tau \right)O\left(\tau \right)+{\lambda }_LL\right)& \left(1.16\right)\\ =\mathrm{diag}\left(\left[A_{1,\mathrm{jq}}\left(t\right),A_{2,\mathrm{jq}}\left(t\right),\dots ,A_{M,\mathrm{jq}}\left(t\right)\right]\right)+{\lambda }_LA_{\mathrm{jq}}\left(t\right)L& \left(1.17\right)\end{array}$

the column-wise BCD leads to the following update for the j-th column of {circumflex over (D)}(t)

$\begin{array}{cc}{\stackrel{\_}{d}}_j:={{\Phi }_{j,j}\left(t\right)}^{-1}\left[b_j\left(t\right)-\sum _{q=1,q\ne j}^Q{\Phi }_{j,q}\left(t\right){\widehat{\stackrel{\_}{d}}}_q\left(t\right)\right]& \left(1.18\right)\\ {\hat{d}}_j\left(t\right)=\frac{{\stackrel{\_}{d}}_j}{\max \left\{{{\stackrel{\_}{d}}_j}_2,1\right\}}.& \left(1.19\right)\end{array}$

An example method for online spectrum prediction is given in Table II.

TABLE II
Input:  online observations {yobs(t)}, {O (t)},
        initial dictionary D0, λs, L, λLand β ∈ (0, 1]
Output: {{circumflex over (π)}(t)}
1:      Set D(0) = D0, A(0) = 0, Am(0) = 0 ∀m ∈  ,
        and B(0) = 0.
2:      for t = 1, 2, . . .
        Perform sparse coding
3:      s                        ^                                                                                          (                        t                        )                                                              =                                                                                            argmin                          s                                                                                                  1                          2                                                                                                                                                                                                                                                                                                                            y                                    obs                                                                                                                                          (                                    t                                    )                                                                                                  -                                                                                                      O                                                                                                              (                                      t                                      )                                                                                                                                                                                                                    D                                      ^                                                                                                                                                  (                                                                              t                                        -                                        1                                                                            )                                                                                                                                          s                                                                                                                                                                                                        2                          2                                                                    +                                                                                                                                                                  λ                          s                                                                                                                                                                                                                      s                                                                                                                                          1                                                                    +                                                                        1                          2                                                                                                  λ                                                      L                                                          S                                                                                                                                                                                                                                                D                              ^                                                                                                                                                                  (                                                          t                              -                              1                                                        )                                                                                                  L                                                                                                      D                            ^                                                                                                          (                                                          t                              -                              1                                                        )                                                                                                  s
        Perform prediction
4:      Output {circumflex over (π)}(t) = {circumflex over (D)}(t − 1)ŝ(t)
        Perform dictionary update
5:      A(t) = βA(t − 1) +  ŝ(t)   (t)
6:      Am(t) = βAm(t − 1) +   ŝ(t)   (t)
        for ∀m ∈
7:      B(t) = βB(t − 1) +   (t)yobs(t)   (t)
8:      Set [{circumflex over (d)}1(t), . . . ,{circumflex over (d)}Q(t)] = {circumflex over (D)}(t − 1)
9:      Repeat
10:     For q = 1, 2, . . . , Q
11:     Update {circumflex over (d)}q(t) as (18)-(19)
12:     Next q
13:     Until convergence
14:     Set {circumflex over (D)}(t) = [{circumflex over (d)}1(t), . . . , {circumflex over (d)}Q(t)]
15:     Next t

Techniques for spatio-temporal spectrum prediction are now described. The algorithms developed so far provide imputations for missing measurements of spatial interference distributions, given the (incomplete) measurements of the current and the past time instants. However, in order to predict the future interference status, for which not even partial measurements can be available, one has to incorporate temporal correlation structures into the model, or learn such structures from the data. In this disclosure, we leverage the data-driven dictionary learning framework to learn temporal dynamics from the data even in an online fashion.

The observations over T consecutive intervals are concatenated into a super-vector, and apply the algorithms described above. That is, define

y ^obs(t):=[y^obsτ(t), . . . , y^obsτ(t−T+1)]^τ  (1.20)

O(t):=[O^τ(t), . . . , O^τ(t−T+1)]^τ  (1.21)

L:=I _T L  (1.22)

where denotes the Kronecker product, which are used in place of y^obs(t), O(t), and L, respectively, in the method of Table II.

To perform prediction for π(t+1), after executing line 14 in Table II, compute sparse coefficient ŝ^f(t+1) for a fictitious observation y^obs,f(t+1):=[y^obsτ(t), . . . , y^obsτ(t−T+2)]^τ, assuming that the entire y^obs(t+1) is missing; i.e., O^f(t+1):=[O^τ(t), . . . O^τ(t−T+2)]^τ. Then, {circumflex over (π)}(t+1) can be obtained as

{circumflex over (π)}(t+1)={circumflex over (D)}(t)[1:M,:]ŝ^f(t+1)  (1.23)

where {circumflex over (D)}(t)[1:M,:] denotes the first M rows of {circumflex over (D)}(t).

The performance of the proposed techniques was evaluated via numerical tests. A CR network consisting of M=20 nodes with the topology depicted in FIG. 11 is considered, where the circles denote the CR nodes and the lines represent the connections established among neighbors. The interference power distribution due to K=3 PU transmitters is also depicted in FIG. 1, where the emitter locations are clearly revealed. The pathloss was computed as

${\left(\frac{d}{d_0}\right)}^{\alpha },$

where d was the distance, d₀=0.01 and α=2.5. The number of atoms of the dictionary was set to Q=50.

First, the two-phase algorithm was tested. Each of the PUs turned on with a 30% chance, and transmitted at a power level p_k(t) chosen from a uniform distribution with support. Each CR made a measurement with a 70% chance. The measurements were corrupted by additive noise generated from a zero-mean Gaussian distribution with variance 10⁻⁵, which was then clipped to ensure non-negativity. The measurements were also normalized by the maximum amplitude observed in the training set. No shadowing or small-scale fading was considered, signifying a quasi-static scenario. N=300 samples were used to train the dictionary, and then another 300 samples were supplied for the operational phase. The values of λ_s and λ_L were set to 0.1 and 0.005, respectively. To compensate for the bias inherent in Lasso-type estimators, de-biasing was performed in the operational phase; that is, after performing sparse coding in (1.9) to obtain ŝ(t), (1.9) was re-solved without the l₁ regularization term only for the non-zero entries in ŝ(t). FIG. 12 shows the true interference level and the reconstructed one at CR m=12 in the thick blue and the thin red curves, respectively. The missing (true) levels are denoted by the circles, whose interpolations are marked by the crosses. It can be seen that the missing entries are accurately recovered through the proposed method.

To test the online algorithm, the Rayleigh fading coefficient h_mk(t) for the channel from PU k to CR m at time t was generated using a first-order autoregressive model

h _mk(t)=αh_mk(t−1)+√{square root over (1−α²)}w_mk(t)  (1.24)

where α=0.9995 was used, and w_mk(t) was circularly symmetric zero-mean complex Gaussian noise with variance 1. The overall channel gain g_mk(t) was formed by multiplying the pathloss with |h_mk(t)|². The forgetting factor was β=0.95, and λ_s=0.25 and λ_L=0.005 were used. The transmit-power of all PUs was fixed to 150 so that the tracking performance could be clearly visible. FIG. 13 shows the interference level for CR 7, where it is evident that the online algorithm tracks the slow variation of interference levels due to channel fading. Also, it is noted that the initial transient for the online learning is quite short.

The normalized root mean square error (RMSE) for the missing observations, averaged over 20 CRs, is depicted in FIG. 14 for various values of λ_s and λ_L. The online algorithm was again used. Two sets of curves are presented, corresponding to the chance of missing observations P_miss equal to 0.3 and 0.5. It can be deduced that the Laplacian-based regularization becomes more important when a larger fraction of observations are unavailable.

In order to test the temporal prediction, certain traffic patterns were assumed. That is, at each time interval t, PU 1 tossed a coin and transmitted with probability 0.1. If PU 1 did transmit, PU 2 transmitted in the next time slot, followed by PU 3's transmission in the third time slot. Likewise, at each time t, PU 3 started transmission with probability 0.15, followed by PU 2 in the second time slot, and PU 1 in the third. FIG. 15 shows the result of the one-time slot-ahead prediction of the interference power at CR 6 using the online algorithm with T=4. No missing measurements were assumed. It can be seen that the traffic patterns are successfully acquired by dictionary learning to predict future interference levels.

Spectrum prediction algorithms for CR networks have been described. Using a dictionary learning framework, the techniques can predict the interference power experienced at each CR node based on the current and the past measurements collected from a subset of nodes in the network. Exploiting the fact that the spatial variation of interference is smooth, a regularization term based on the CR network topology was also incorporated. Batch and online algorithms were derived, where the online alternative possessed a tracking capability at lower complexity and memory requirements. Temporal prediction was also discussed. Numerical tests verified the efficacy of the techniques.

FIG. 16 shows a detailed example of various devices that may be configured to execute program code to practice some embodiments in accordance with the current disclosure. For example, device 500 may be a CR 12, a FC 16, a workstation, a computing center, a cluster of servers or other example embodiments of a computing environment, centrally located or distributed, capable of executing the techniques described herein.

In this example, a computer 500 includes a processor 510 that is operable to execute program instructions or software, causing the computer to perform various methods or tasks. Processor 510 is coupled via bus 520 to a memory 530, which is used to store information such as program instructions and other data while the computer is in operation. A storage device 540, such as a hard disk drive, nonvolatile memory, or other non-transient storage device stores information such as program instructions, data files of the multidimensional data and the reduced data set, and other information. The computer also includes various input-output elements 550, including parallel or serial ports, USB, Firewire or IEEE 1394, Ethernet, and other such ports to connect the computer to external device such a printer, video camera, surveillance equipment or the like. Other input-output elements include wireless communication interfaces such as Bluetooth, Wi-Fi, and cellular data networks.

The computer itself may be a traditional personal computer, a rack-mount or business computer or server as shown in FIG. 16, or any other type of computerized system. The computer in a further example may include fewer than all elements listed above, such as a thin client or mobile device having only some of the shown elements. In another example, the computer is distributed among multiple computer systems, such as a distributed server that has many computers working together to provide various functions.

The techniques described herein may be implemented in hardware, software, firmware, or any combination thereof. Various features described as modules, units or components may be implemented together in an integrated logic device or separately as discrete but interoperable logic devices or other hardware devices. In some cases, various features of electronic circuitry may be implemented as one or more integrated circuit devices, such as an integrated circuit chip or chipset.

If implemented in hardware, this disclosure may be directed to an apparatus such a processor or an integrated circuit device, such as an integrated circuit chip or chipset. Alternatively or additionally, if implemented in software or firmware, the techniques may be realized at least in part by a computer readable data storage medium comprising instructions that, when executed, cause one or more processors to perform one or more of the methods described above. For example, the computer-readable data storage medium may store such instructions for execution by a processor. Any combination of one or more computer-readable medium(s) may be utilized.

A computer-readable medium may form part of a computer program product, which may include packaging materials. A computer-readable medium may comprise a computer data storage medium such as random access memory (RAM), read-only memory (ROM), non-volatile random access memory (NVRAM), electrically erasable programmable read-only memory (EEPROM), flash memory, magnetic or optical data storage media, and the like. In general, a computer-readable storage medium may be any tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device. Additional examples of computer readable medium include computer-readable storage devices, computer-readable memory, and tangible computer-readable medium. In some examples, an article of manufacture may comprise one or more computer-readable storage media.

In some examples, the computer-readable storage media may comprise non-transitory media. The term “non-transitory” may indicate that the storage medium is not embodied in a carrier wave or a propagated signal. In certain examples, a non-transitory storage medium may store data that can, over time, change (e.g., in RAM or cache).

The code or instructions may be software and/or firmware executed by processing circuitry including one or more processors, such as one or more digital signal processors (DSPs), general purpose microprocessors, application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other equivalent integrated or discrete logic circuitry. Accordingly, the term “processor,” as used herein may refer to any of the foregoing structure or any other processing circuitry suitable for implementation of the techniques described herein. In addition, in some aspects, functionality described in this disclosure may be provided within software modules or hardware modules.

As described herein, spatio-temporal and dynamic re-use of the licensed bands calls for collaborative CR network sensing algorithms able to portray the ambient power spectral density at arbitrary locations in space, frequency, and time. The present paper addressed this ambitious task through a parsimonious model of the PSD in frequency and space, which reduces the sensing task to estimating a sparse vector of unknown parameters. An estimator of the model parameters was developed based on the GS-Lasso, and a low-complexity solver based on the ADMoM was presented. The location and transmit-PSD information conveyed by the estimated model parameters, complemented with either a path loss-based or more elaborated propagation models was shown to allow CRs to accurately reconstruct the PSD atlas of the primary system. To cope with uncertainty in the regression matrix, a provably convergent sensing algorithm was introduced which combines the merits of the TLS framework with the hierarchical sparsity inherent to the network-level sensing problem. To account for outliers, a robust algorithm able to discern and reject unreliable PSD data was also developed. The novel robust GS-TLS approach capitalizes on the sparsity of the unknown parameters and the outliers, and offers systematic estimation of the spectrum holes jointly in space, frequency, and time while taking into account channel uncertainties and unmodeled errors.

Further exemplary details are described in: Dall'Anese, “Group sparse Lasso for cognitive network sensing robust to model uncertainties and outliers,” Physical Communication, November, 2011; Dall'Anese, “Group sparse Lasso for cognitive network sensing robust to model uncertainties and outliers,” IEEE 12^th International Workshop on Signal Processing Advances in Wireless Communications, Jun. 26, 2011; Dall'Anese, “Distributed Cognitive Spectrum Sensing via Group Sparse Total Least-Squares,” 4^th IEEE International Workshop on Computational Advances in Multi-sensor Adaptive Processing,” Dec. 13, 2011; and Kim, “Joint Link Learning and Cognitive Radio Sensing,” in Proc. of the 45th Asilomar Conf. on Signals, Systems, and Computers, Pacific Grove, Calif., Nov. 6-9, 2011, contents of all of which are incorporated herein by reference.

The following clauses illustrate additional examples:

Clause 1. A method comprising: sensing local radio-frequency (RF) interference spectrum at each of a plurality of sensors positioned at a plurality of locations within a geographic region; and computing a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location, wherein the basis expansion model is computed as a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band, and wherein computing the basis expansion model comprises applying group sparse (GS) total least-squares (TLS) method to compute the coefficients.

Clause 2. The method of clause 1, wherein computing estimates for coefficients of the basis expansion model comprises applying an alternating direction method of multipliers (ADMoM) to compute the coefficients.

Clause 3. The method of clause 1, wherein computing a basis expansion model comprises computing the coefficients of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso).

Clause 4. A system comprising: a plurality of sensors to sense a local radio-frequency (RF) interference spectrum at each of a plurality of locations within a geographic region; and a processor that computes a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location, wherein the processor computes the basis expansion model is computed as a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band, and wherein the processor applies group sparse (GS) total least-squares (TLS) method to compute the coefficients.

Clause 5. The system of clause 4, wherein the processor applies an alternating direction method of multipliers (ADMoM) to compute the coefficients.

Clause 6. The system of claim 4, wherein the processor computes the coefficients of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso) to identify the coefficients for computation.

Clause 7. A method comprising: sensing, over a period of time, local radio-frequency (RF) interference spectrum at a subset of a plurality of sensors positioned at a plurality of locations within a geographic region; and applying dictionary learning to compute predicted interference power levels at the plurality of sensors for a current time based on current and past measurements collected from the subset of sensor in the network.

Clause 8. The method of clause 7, further comprising constructing a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location based on the predicted interference power at each of the sensors.

Various embodiments of the invention have been described. These and other embodiments are within the scope of the following claims.

Claims

1. A method comprising: sensing local radio-frequency (RF) interference spectrum at each of a plurality of sensors positioned at a plurality of locations within a geographic region;computing a basis expansion model from the sensed RF interference spectrum, wherein the basis expansion model comprises a set of reference basis functions that represent a distribution of RF power emitted by RF-enabled devices throughout the geographic region as a function of frequency and location;computing coefficients for the reference basis functions of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso); andconstructing, in accordance with the basis expansion model, a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region.

2. The method of claim 1, wherein computing the coefficients for the reference basis functions of the basis expansion model using a GS-Lasso comprises applying a GS-Lasso operator to select which coefficients of the basis expansion model are required for computation.

3. The method of claim 1, wherein computing the coefficients for the reference basis functions of the basis expansion model using a GS-Lasso comprises applying an alternating direction method of multipliers (ADMoM) to compute the coefficients.

4. The method of claim 1, wherein computing the coefficients for the reference basis functions of the basis expansion model using a GS-Lasso comprises applying a group sparse (GS) total least-squares (TLS) method to identify and remove outlier measurements from the sensed interference spectrum when computing the coefficients.

5. The method of claim 1, wherein computing the basis expansion model comprises computing the basis expansion model as a parametric model in which each of the coefficients for the reference basis functions is a scalar.

6. The method of claim 1, further comprising computing the reference basis functions to represent the aggregate distribution of RF power across the geographic region corresponding to frequencies spanned by the reference basis functions.

7. The method of claim 1, wherein computing the basis expansion model comprises computing a parsimonious model that accounts for scarce presence of transmitting RF-enabled devices by using only a small subset of the coefficients for the model as non-zero as identified by the GS-Lasso.

8. The method of claim 1, wherein computing the basis expansion model comprises performing a most parsimonious sparse signal expansion using an overcomplete basis to construct the PSD maps.

9. The method of claim 1, further comprising: communicating observation data indicative of the sensed RF interference spectrum from each of the sensors to a centralized computer; andcomputing the PSD map with the centralized computer based on the observation data.

10. The method of claim 1, further comprising: exchanging, between the sensors, observation data indicative of the sensed RF interference spectrum; andcomputing the PSD map with the sensors based on the observation data.

11. The method of claim 1, further comprising: processing the PSD map to identify a location within the geographic region where at least one frequency is unoccupied.

12. A system comprising: a plurality of sensors to sense a local radio-frequency (RF) interference spectrum at each of a plurality of locations within a geographic region; anda processor that computes a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location,wherein the processor computes the basis expansion model as a plurality of reference basis functions each having a corresponding coefficient, the reference basis functions representing a power emitted by an RF-enabled device on a corresponding frequency band, andwherein the processor computes the coefficients of the basis expansion model using a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso).

13. The system of claim 12, wherein the processor applies an alternating direction method of multipliers (ADMoM) to compute the coefficients.

14. The system of claim 12, wherein the processor applies a group sparse (GS) total least-squares (TLS) method to identify and remove outlier measurements from the sensed interference spectrum when computing the coefficients.

15. The system of claim 12, wherein the processor computes the reference basis functions to represent RF signals transmitted by RF-enabled devices within the geographic region

16. The system of claim 12, wherein the processor computes the basis expansion model as a parametric model in which each of the coefficients for the reference basis functions are scalars.

17. The system of claim 12, wherein the processor computes the reference basis functions to represent the aggregate distribution of RF power across the geographic region corresponding to frequencies spanned by the basis functions.

18. The system of claim 12, wherein the processor computes the basis expansion model by computing a parsimonious model that accounts for scarce presence of transmitting RF-enabled devices by using only a small subset of the coefficients for the model as non-zero

19. The system of claim 12, wherein the processor computes the basis expansion model by performing a most parsimonious sparse signal expansion using an overcomplete basis to construct the PSD maps

20. The system of claim 12, wherein the processor receives observation data indicative of the sensed RF interference spectrum from each of the sensors, and computes the PSD map with the processor based on the observation data.

21. The system of claim 12, wherein the sensors exchange observation data indicative of the sensed RF interference spectrum, andwherein the processor comprises one of a plurality of processors located at the sensors to compute the basis expansion model and construct the PSD map based on the observation data.

22. A mobile device comprising: a sensor to sense a local radio-frequency (RF) interference spectrum at a locations within a geographic region; anda processor that computes a basis expansion model from the sensed RF interference spectrum at each of the sensors to construct a power spectral density (PSD) map representative of the distribution of RF power throughout the geographic region as a function of frequency and location,wherein the processor computes the basis expansion model to include a plurality of functions having a corresponding coefficient, each of the functions representing a power emitted by an RF-enabled device on a corresponding frequency band, andwherein the processor applies a group sparse least-absolute-shrinkage-and-selection operator (GS-Lasso) to coefficients of the basis expansion model only for frequencies an locations within the geographic region having an active transmitter.