Three-Event Energy Detection with Adaptive Threshold for Spectrum Sensing in Cognitive Radio Systems

Abstract

Implementation of dynamic spectrum access (DSA) in cognitive radio (CR) systems requires the unlicensed secondary users (SU) to implement spectrum sensing to monitor the activity of the licensed primary users (PU). Energy detection (ED) is one of the most widely used methods for spectrum sensing in CR systems, and in this paper we present a novel ED algorithm with an adaptive sensing threshold. The three-event ED (3EED) algorithm for spectrum sensing is considered for which an accurate approximation of the optimal decision threshold that minimizes the decision error probability (DEP) is found using Newton’s method with forced convergence in one iteration. The proposed algorithm is analyzed and illustrated with numerical results obtained from simulations that closely match the theoretical results and show that it outperforms the conventional ED (CED) algorithm for spectrum sensing.

1. Introduction

The expansion of wireless communication systems in all sectors of modern society over the past decade has prompted an increased demand in spectrum resources. In this context, the traditional static allocation of the radio frequencies leads to inefficient use of the spectrum [1,2], and cognitive radio (CR) systems have emerged as a meaningful alternative for improving the efficiency of spectrum usage through dynamic spectrum access (DSA) [3]. DSA allows the access of secondary users (SU) to licensed frequency bands when these are not actively used by licensed primary users (PU), and relies on spectrum sensing, which is performed by the SU to detect the presence of active PU.

Among the methods used for spectrum sensing, energy detection (ED) [4,5,6,7,8] is the most commonly used one, due to its simplicity as well as its almost universal applicability. Alternative spectrum sensing approaches use eigenvalue-based algorithms [9], covariance-based detection methods [10,11], cyclostationary feature detection [12,13], or compressive sensing [14,15]. A thorough and very recent survey on most known spectrum sensing methods can be found in [16].

An important drawback of the ED method is implied by its sensitivity to noise uncertainty [17], which has led to improved ED algorithms that outperform the classical energy detection (CED) method [4]. These include the modified ED method for spectrum sensing in [18], which uses the average value of the ED test statistic over multiple sensing events, as well as the three-event ED (3EED) algorithm in [19], which decides that a PU is active if the ED test statistic exceeds the sensing threshold in three consecutive sensing events that include the current sensing event along with the sensing events immediately before and after it. To improve the performance of ED in the presence of Laplacian noise, [20] proposes a spectrum sensing algorithm in which the square of the received signal amplitude in the ED test statistic is replaced by an exponent that ranges between 0 and 2. Furthermore, adaptive ED approaches have been established for dynamic scenarios, which adjust the duration of the sensing window based on the PU on/off activity [21] or change the sensing threshold in response to changes in the noise characteristics [22,23].

In this paper, we present a new adaptive ED algorithm for spectrum sensing, which is based on the 3EED algorithm [19], but in which the sensing threshold is adapted similar to [22], to optimize the decision error probability (DEP), which is a weighted sum of the probabilities of missed detection and false alarm [22,23]. We note that obtaining an accurate analytical expression for the optimal threshold is elusive as the expressions involved in finding the sensing threshold do not have closed-form expressions and require approximations. In order to underline the novelty of the current work, we will enumerate its contributions as compared to the previous similar and related works. In [22], the authors proposed a threshold adaptation method by minimizing the DEP under the Gaussian approximation, but for the simple CED algorithm. Hence, an exact expression of the threshold could be obtained directly as a solution of the optimization equation, using the properties of the Q-function. However, for more efficient algorithms than CED, regularly having more complex expressions for DEP, the optimization equation is not analytically solvable. Therefore, in this paper, we aim at extending the method from [22] to a more complex and efficient ED algorithm, such as 3EED [19]. Under the Gaussian approximation, the DEP for any ED algorithm can be written as an expression based on several Q-function terms. First, we have to prove that DEP is a convex function in the threshold value and then, we propose the use of a numerical method, such as Newton’s method, to iteratively determine the root of the analytically unsolvable optimization equation. However, the main drawback of iterative methods is the increased operating time, which is a critical issue for the spectrum sensing algorithms. In order to overcome this problem, we propose in this paper a transformation of the optimization function, such that the Newton’s method for the transformed function converges faster. By analyzing the monotonicity of the optimization function, we determine a transform that linearizes the function around the root. Thus, we manage to reduce the number of iterations to the minimum possible, i.e., one iteration. As a prior validation step for the current research, in [24], we tested only empirically the performance of the adaptive threshold 3EED algorithm, where we derived the threshold value by using a brute-force search. Finally, in the current paper, we derive an analytical approximate expression of the adaptive threshold for 3EED by using a fast convergence numerical method. Moreover, we consider that this method can be generalized for most if not all ED-based spectrum sensing algorithms.

The rest of the paper is structured as follows: in Section 2, we briefly present the CED and 3EED algorithms for spectrum sensing, followed by a discussion on optimizing the sensing threshold in Section 3. The proposed 3EED algorithm with an adaptive threshold is formally stated in Section 4. In Section 5, we illustrate the performance of the proposed algorithm with numerical results obtained from simulation and compare it to that of the adaptive CED in [22], and we conclude the paper with final remarks in Section 6.

2. Spectrum Sensing by Three-Event Energy Detection (3EED)

Let $y\left(n\right)$ be the signal at the SU receiver, which consists of an active PU signal with average power ${\sigma }_s^2$ corrupted by additive white Gaussian noise (AWGN) variable with variance ${\sigma }_n^2$ when the channel is busy, or just of the AWGN when the channel is idle (PU is not active). We denote by $E_i$ the value of the received signal energy estimated during the i-th sensing slot that consists of N samples, $\lambda$ the ED sensing threshold, and $q_i$ the binary decision variable {0,1} for the i-th spectrum sensing slot. With these notations, in CED correct detection of the PU signal (active/inactive) implies setting $q_i=1$, that is, the channel is “busy”, if $E_i>\lambda$ when $y\left(n\right)$ is implied by an active PU signal corrupted by noise, and $q_i=0$, that is, the channel is “idle”, if $E_i\le \lambda$ when $y\left(n\right)$ corresponds to a noise only term. We note that, when the number of samples N is sufficiently large, the tail region of the probability density function of the received signal energy $E_i$ for values above the sensing threshold is approximated by a normal distribution [25] such that the standard Q-function [26] can be used to approximate the probabilities of detection and false alarm for CED as [18]

$$P_{fa}^{CED}\left(\lambda ,{\sigma }_n^2\right)=\mathrm{Prob}[E_i>\lambda |\mathrm{channel} \mathrm{is} “\mathrm{idle}”]=Q\left(\frac{\frac{\lambda }{N}-{\sigma }_n^2}{{\sigma }_n^2\sqrt{\frac{2}{N}}}\right)$$

(1)

and

$$P_d^{CED}\left(\lambda ,{\sigma }_s^2,{\sigma }_n^2\right)=\mathrm{Prob}[E_i>\lambda |\mathrm{channel} \mathrm{is} “\mathrm{busy}”]=Q\left[\frac{\frac{\lambda }{N}-{\sigma }_s^2-{\sigma }_n^2}{\left({\sigma }_s^2+{\sigma }_n^2\right)\sqrt{\frac{2}{N}}}\right].$$

(2)

As the name suggests, the 3EED algorithm [19] relies on ED in up to three consecutive sensing slots to decide if the signal $y\left(n\right)$ at the SU receiver is implied by an active PU signal or not. Specifically, similar to CED, if the energy in the current spectrum sensing slot i exceeds the sensing threshold, $E_i>\lambda$, then the 3EED algorithm decides that a PU is active and sets $q_i=1$. However, if $E_i\le \lambda$, the algorithm checks the energy level estimated in the previous slot, $E_{i-1}$, to decide that a PU signal is active and set $q_i=1$ if $E_{i-1}>\lambda$, or to move on to the next time slot $i+1$ to estimate $E_{i+1}$ and make the final decision that PU signal is active setting $q_i=1$ if $E_{i+1}>\lambda$, or inactive setting $q_i=0$ if $E_{i+1}\le \lambda$. Thus, if the PU signal is identified as active in any of the three consecutive sensing slots, i, $i-1$, and $i+1$, the algorithm returns $q_i=1$ corresponding to a “busy” channel, while if the PU signal is determined to be inactive in all three consecutive time slots the algorithm returns $q_i=0$ corresponding to an “idle” channel. The formal statement of the 3EED algorithm is given in Algorithm 1.

Algorithm 1 The 3EED Algorithm for Spectrum Sensing

Input: sensing threshold $\lambda$

Output: output result

1: for each spectrum sensing slot i do 2:     Estimate received signal energy $E_i$ and save its value 3:     if $E_i>\lambda$ then 4:         Set $q_i=1$⟶ PU active/Channel “busy” 5:     else 6:         Read $E_{i-1}$ (saved in slot i-1) 7:         if $E_{i-1}>\lambda$ then 8:            Set $q_i=1$⟶ PU active/Channel “busy” 9:         else 10:            Estimate $E_{i+1}$ 11:            if $E_{i+1}>\lambda$ then 12:                Set $q_i=1$⟶ PU active/Channel “busy” 13:            else 14:                Set $q_i=0$⟶ PU not active/Channel “idle” 15: return$q_i$

The probabilities of false alarm and detection for the 3EED algorithm for spectrum sensing are expressed in terms of the probabilities of false alarm and detection for the CED algorithm, $P_{fa}^{CED}$ and $P_d^{CED}$, respectively as [19]:

$$\begin{array}{ccc}\hfill P_{fa}^{3EED}& =& P_{fa}^{CED}+\left(1-P_{fa}^{CED}\right)·\hfill \\ & & \left(\frac{T-B-1}{T-B}{P}_{fa}^{CED}+\frac{1}{T-B}{P}_d^{CED}\right)·\hfill \\ & & \left[1+\frac{T-B-1}{T-B}\left(1-P_{fa}^{CED}\right)+\frac{1}{T-B}\left(1-P_d^{CED}\right)\right]\hfill \end{array}$$

(3)

and

$$\begin{array}{ccc}\hfill P_d^{3EED}& =& P_d^{CED}+\left(1-P_d^{CED}\right)·\hfill \\ & & \left(\frac{1}{B}{P}_{fa}^{CED}+\frac{B-1}{B}{P}_d^{CED}\right)·\hfill \\ & & \left[1+\frac{1}{B}\left(1-P_{fa}^{CED}\right)+\frac{B-1}{B}\left(1-P_d^{CED}\right)\right],\hfill \end{array}$$

(4)

where T represents the total duration of the transmission cycle that consists of B consecutive slots during which the channel is busy followed by a number of $T-B$ slots in which the channel is idle. We note that the ratio $\alpha =B/T$ represents the spectrum utilization ratio of the PU (that is the fraction of time the channel is actively used by the PU).

The decision threshold $\lambda$ is determined based on a desired performance level that is specified in terms of a target probability of false alarm value, $P_{fa}$, as it is usually the case with constant false alarm rate (CFAR) detectors [19],

$$\lambda =\left[Q^{-1}\left(1+\sqrt[3]{P_{fa}-1}\right)\sqrt{2N}+N\right]{\sigma }_n^2.$$

(5)

We note that the threshold must be adapted when changes in the noise variance occur, in order to keep the probability of false alarm at the desired value. We also note that, by considering only the probability of false alarm in setting the sensing threshold, the SU is favored, since by setting the value of $P_{fa}$ low the SU has more chances of utilizing the spectrum. However, in the case of a mis-detection, when the PU is active but the channel is incorrectly detected as “idle”, the SU transmission will interfere with the active PU signal and will negatively affect the PU performance. Thus, in order for spectrum sensing to be beneficial to both PU and SU, threshold setting should consider both the probability of false alarm and the probability of mis-detection.

3. Optimizing the 3EED Sensing Threshold to Minimize the Decision Error Probability

To optimize the decision threshold of the 3EED algorithm for spectrum sensing we use the DEP as the performance metric, which is formally defined as [22]:

$$P_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right)=\left(1-\alpha \right)P_{fa}^{3EED}+\alpha \left(1-P_d^{3EED}\right).$$

(6)

Analyzing expression (6) we note that the DEP is a function of the sensing threshold $\lambda$, the PU spectrum utilization ratio $\alpha$, the average power of the PU signal ${\sigma }_s^2$, and the noise variance ${\sigma }_n^2$. Furthermore, both terms in the DEP expression (6) are related to decision errors in the spectrum sensing process:

• The first term in (6) is implied by the probability of false alarm and corresponds to the case when a given sensing slot is found busy without an active PU transmission in this slot. False alarms occur only when the channel is idle and the term is weighted by the $1-\alpha$ factor corresponding to the fraction of time when the PU is not active.
• The second term in (6) corresponds to the probability of mis-detection and is associated with the case when a given sensing slot is found idle while an active PU transmission is actually in progress in this slot. Mis-detections occur when the channel is busy and the term is weighted by the $\alpha$ factor corresponding to the fraction of time when the PU is active.

Thus, the DEP provides a trade-off between the false alarm and mis-detection performance, and allows optimization of the decision threshold for the 3EED algorithm by setting it to minimize the DEP (6), that is

$${\lambda }_{opt}=arg\underset{\lambda }{min}{P}_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right).$$

(7)

Introducing (3) and (4) in (6), we rewrite the expression of DEP as

$$\begin{array}{ccc}\hfill P_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right)& =& \left(1-\alpha \right)\left\{1+{\left[P_{fa}^{CED}\left(\lambda ,{\sigma }_n^2\right)-1\right]}^3\right\}\hfill \\ & +& \alpha {\left[1-P_d^{CED}\left(\lambda ,{\sigma }_s^2,{\sigma }_n^2\right)\right]}^3,\hfill \end{array}$$

(8)

For the sake of simplicity, let us introduce two variables denoted as $a(\lambda ,{\sigma }_n^2)$ and $b(\lambda ,{\sigma }_s^2,{\sigma }_n^2)$, which are the arguments of the Q-function in (1) and (2), respectively, and expressed as

$$a(\lambda ,{\sigma }_n^2)=\frac{\lambda }{\sqrt{2N}{\sigma }_n^2}-\sqrt{\frac{N}{2}}$$

(9)

and

$$b(\lambda ,{\sigma }_s^2,{\sigma }_n^2)=\frac{\lambda }{\sqrt{2N}\left({\sigma }_s^2+{\sigma }_n^2\right)}-\sqrt{\frac{N}{2}}.$$

(10)

Using the expressions (9) and (10), we can rewrite the DEP function from (8) as following:

$$\begin{array}{ccc}\hfill P_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right)& =& \left(1-\alpha \right)\left\{1+{\left[Q\left(a\right)-1\right]}^3\right\}+\alpha {\left[1-Q\left(b\right)\right]}^3,\hfill \end{array}$$

(11)

We note that finding an accurate closed-form expression for ${\lambda }_{opt}$ is elusive as the expression of the DEP metric to be optimized has already been approximated when expressions of $P_{fa}^{3EED}$ and $P_d^{3EED}$ in terms of the Q-functions have been used to obtain (11), and resorting to further approximations of the Q-function in terms of simpler elementary functions will affect the accuracy of the optimal threshold value obtained. Thus, we focus on numerical optimization of the DEP metric (11), which takes advantage of its convexity properties.

Specifically, we notice that the DEP function has four variables, i.e., $\lambda$, $\alpha$, ${\sigma }_s^2$, and ${\sigma }_n^2$. However, the optimization aims only at the threshold variable, $\lambda$. Therefore, the other three variables are assumed to be constant for this optimization. In order to perform the threshold optimization, the DEP function from (6) and (11) has to be a convex function in $\lambda$. Therefore, we provide the next theorem.

Theorem 1.

The DEP function given by (11) is a convex function in λ.

Proof. 

Let us introduce the following notations for the functional terms from (11):

$$\begin{array}{ccc}\hfill T_1\left(a\right)& =& 1+{\left[Q\left(a\right)-1\right]}^3\hfill \\ \hfill T_2\left(b\right)& =& {\left[1-Q\left(b\right)\right]}^3\hfill \end{array}$$

(12)

Analyzing the second derivative of both functions from (12) with respect to a and b, respectively, one can notice that both have a unique and identical inflection point (i.e., the solution of $T_1^{″}\left(a\right)=0$ and $T_2^{″}\left(b\right)=0$). The value of this common inflection point is the solution of the following transcendental equation for $x=a$ or $x=b$:

$$\begin{array}{ccc}\hfill 2Q^{\prime }\left(x\right)& =& x\left[Q\left(x\right)-1\right]\hfill \end{array}$$

(13)

An approximate value for the solution of the Equation (13) can be obtained using numerical methods or graphically, i.e., $x_0=0.76527\dots$. Considering this approximate value of $x_0$ and using the relations in (9) and (10), we can determine the values of the inflection point values of ${\lambda }_{inflec}^{\left(1\right)}$ and ${\lambda }_{inflec}^{\left(2\right)}$ for the functions $T_1\left(a\right)$ and $T_2\left(b\right)$, respectively, from (12):

$$\begin{array}{ccc}\hfill {\lambda }_{inflec}^{\left(1\right)}& =& \sqrt{2N}{\sigma }_n^2\left(x_0+\sqrt{\frac{N}{2}}\right)\hfill \\ \hfill {\lambda }_{inflec}^{\left(2\right)}& =& \sqrt{2N}\left({\sigma }_s^2+{\sigma }_n^2\right)\left(x_0+\sqrt{\frac{N}{2}}\right)\hfill \end{array}$$

(14)

In order to evaluate the convexity of the functions from (12), we have to estimate the values in the inflection points, given by (14), of the first derivative of these functions with respect to $\lambda$, respectively, denoted as $T_1^{\prime }\left(\lambda ={\lambda }_{inflec}^{\left(1\right)}\right)$ and $T_2^{\prime }\left(\lambda ={\lambda }_{inflec}^{\left(2\right)}\right)$. Analyzing these first derivative functions, we noticed that $T_1^{\prime }\left({\lambda }_{inflec}^{\left(1\right)}\right)<0$ and $T_2^{\prime }\left({\lambda }_{inflec}^{\left(2\right)}\right)>0$. Therefore, it results that $T_1\left(\lambda \right)$ is convex in $\lambda$ for $\lambda \ge {\lambda }_{inflec}^{\left(1\right)}$ and $T_2\left(\lambda \right)$ is convex in $\lambda$ for $\lambda \le {\lambda }_{inflec}^{\left(2\right)}$.

Also, we note that the DEP expression from (11), with the notations from (12), is a linear combination of the two convex functions $T_1\left(\lambda \right)$ and $T_2\left(\lambda \right)$ for the intersection set ${\lambda }_{inflec}^{\left(1\right)}\le \lambda \le {\lambda }_{inflec}^{\left(2\right)}$:

$$\begin{array}{ccc}\hfill P_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right)& =& \left(1-\alpha \right)T_1\left(\lambda ,\alpha ,{\sigma }_n^2\right)+\alpha T_2\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right),\hfill \end{array}$$

(15)

We also know that the spectrum utilization ratio of the PU is a positive number, i.e., $\alpha \in \left[0,1\right]$ and therefore, both coefficients of the linear expression from (15) are positive numbers.

Finally, Theorem $2.10$ from page 52 in [27] states that a linear combination of two convex functions on a set, with non-negative coefficients, is also a convex function. Hence, we can conclude that the DEP expression from (11) is a convex function in $\lambda$ on the set ${\lambda }_{inflec}^{\left(1\right)}\le \lambda \le {\lambda }_{inflec}^{\left(2\right)}$. □

Now, because the DEP function is convex as stated by Theorem 1, we can determine the optimum threshold from (7) by solving the following equation to find the critical points of the DEP (6):

$$\frac{\partial P_e\left(\lambda ,\alpha ,{\sigma }_s^2,{\sigma }_n^2\right)}{\partial \lambda }=0$$

(16)

Next, using the DEP expression from (11) we obtain

$$\begin{array}{ccc}\hfill \frac{\partial P_e}{\partial \lambda }& =& 3\left(1-\alpha \right){\left[Q\left(a\right)-1\right]}^2\frac{dQ\left(a\right)}{da}\frac{\partial a}{\partial \lambda }\hfill \\ & -& 3\alpha {\left[Q\left(b\right)-1\right]}^2\frac{dQ\left(b\right)}{db}\frac{\partial b}{\partial \lambda }\hfill \end{array}$$

(17)

where $a(\lambda ,{\sigma }_n^2)$ and $b(\lambda ,{\sigma }_s^2,{\sigma }_n^2)$ were previously introduced in (9) and (10).

Using the expression of the first derivative of the Q-function

$$\frac{dQ\left(x\right)}{dx}=-\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}$$

(18)

the threshold optimization Equation (16) becomes

$$\frac{Q\left(a\right)-1}{Q\left(b\right)-1}=A(\alpha ,\gamma )e^{-(b^2-a^2)/4}$$

(19)

where $\gamma ={\sigma }_s^2/{\sigma }_n^2$ is the signal-to-noise ratio (SNR) and

$$A(\alpha ,\gamma )=\sqrt{\frac{\alpha }{1-\alpha }}·\sqrt{\frac{1}{1+\gamma }}.$$

(20)

We note that Equation (19) is transcendental, which makes it difficult to find a closed-form solution for the optimal threshold ${\lambda }_{opt}$, and the alternative is to obtain a numerical solution for it as discussed in the following section.

4. Numerical Approach to Finding the Optimal 3EED Sensing Threshold

The numerical approach used to find ${\lambda }_{opt}$, the optimal threshold for the 3EED spectrum sensing, is based on the application of Newton’s iterative method [28] for Equation (19) rewritten as

$$\underset{f\left(\lambda \right)}{\underbrace{\frac{Q\left(a\right)-1}{Q\left(b\right)-1}-A(\alpha ,\gamma )e^{-(b^2-a^2)/4}}}=0.$$

(21)

As the rate of convergence for Newton’s method is quadratic, with proper initialization, a solution to (21) that approximates well the optimal sensing threshold ${\lambda }_{opt}$ can be found in only a few iterations of the type [28]

$${\lambda }_{n+1}={\lambda }_n-\frac{f\left({\lambda }_n\right)}{f^{\prime }\left({\lambda }_n\right)}.$$

(22)

Furthermore, the proposed approach uses the Q-function approximation [29]

$$Q\left(x\right)\cong \frac{\left(1-e^{-1.4x}\right)e^{-x^2/2}}{1.135\sqrt{2\pi }x},x\ge 0$$

(23)

to rewrite the $Q\left(a\right)$ and $Q\left(b\right)$ terms in (21).

We have to note that the initialization of Newton’s iterative method, used to solve the optimization equation from (21), and its convergence are strongly related to the convexity of the DEP function. According to the DEP function’s convexity demonstrated in Theorem 1, we know that the best initialization values of $\lambda$ should be the values of the inflection points, i.e., either ${\lambda }_{inflec}^{\left(1\right)}$ or ${\lambda }_{inflec}^{\left(2\right)}$ given in (14). However, for the analytical derivation of the optimum value of $\lambda$, the use of ${\lambda }_{inflec}^{\left(1\right)}$ or ${\lambda }_{inflec}^{\left(2\right)}$ in the expressions of the optimization function from (21) and its derivative function would have complex expressions. Instead, we prefer to use the values of $\lambda$ for which the argument variables a and b in the optimization Equation (21) take null values, i.e., $a=0$ or $b=0$. Hence, all terms multiplied by a or b will be eliminated from the final expressions. The $\lambda$ values that are the solutions for the equations $a(\lambda ,{\sigma }_n^2)=0$ and $b(\lambda ,{\sigma }_s^2,{\sigma }_n^2)=0$, denoted as ${\lambda }_{a,0}$ and ${\lambda }_{b,0}$, respectively, result from (9) and (10) as

$${\lambda }_{a,0}=N{\sigma }_n^2$$

(24)

and

$${\lambda }_{b,0}=N{\sigma }_n^2\left(1+\gamma \right)$$

(25)

It is important to note that the difference between the values of ${\lambda }_{a,0}$ and ${\lambda }_{b,0}$ and the values of ${\lambda }_{inflec}^{\left(1\right)}$ and ${\lambda }_{inflec}^{\left(2\right)}$ is almost neglectable. Using the expressions from (14), (24) and (25), it can be easily demonstrated that the relative difference between these values is given by:

$$\begin{array}{c}\hfill \frac{{\lambda }_{inflec}^{\left(2\right)}-{\lambda }_{b,0}}{{\lambda }_{inflec}^{\left(2\right)}}=\frac{{\lambda }_{inflec}^{\left(1\right)}-{\lambda }_{a,0}}{{\lambda }_{inflec}^{\left(1\right)}}=\frac{x_0}{x_0+\sqrt{\frac{N}{2}}}\end{array}$$

(26)

Considering that $x_0=0.76527\dots$ and N takes values of the order of tens of thousands samples, then it results in a value for the relative difference in (26) that is less than $1\%$.

Hence, considering all the above, we will use the values of ${\lambda }_{a,0}$ and ${\lambda }_{b,0}$ as initial values for Newton’s method.

Since (23) is valid only for positive arguments of the Q-function and the arguments $a(\lambda ,{\sigma }_n^2)$ and $b(\lambda ,{\sigma }_s^2,{\sigma }_n^2)$ in (21) can take both positive and negative values, as can be observed from (9) and (10), the approximation is applied in conjunction with the symmetry property of the Q-function, $Q(-x)=1-Q\left(x\right)$. Specifically, the following three distinct cases will be analyzed separately to provide a numerical solution for Equation (21): $\left(i\right)b<a\le 0$, $\left(ii\right)a>0,b\le 0$, and $\left(iii\right)a>b\ge 0$.

We have to note that the domain of $\lambda$ values is separated into these three sub-domains ($\left(i\right)$, $\left(ii\right)$, and $\left(iii\right)$), where the value of ${\lambda }_{a,0}$ delimits cases $\left(i\right)$ and $\left(ii\right)$ and ${\lambda }_{b,0}$ delimits cases $\left(ii\right)$ and $\left(iii\right)$. Therefore, ${\lambda }_{a,0}$ will be used as the initial value for the iterative method in case $\left(i\right)$ and ${\lambda }_{b,0}$ will be used as the initial value for the iterative method in case $\left(iii\right)$. Instead, for case $\left(ii\right)$, we can use either ${\lambda }_{a,0}$ or ${\lambda }_{b,0}$ as the initial value, because the minimum DEP can be reached from either side of the interval with a similar accuracy.

4.1. Case $\left(i\right)b<a\le 0$

Using the symmetry property of the Q-function, Equation (21) can be rewritten as:

$$f\left(\lambda \right)=\frac{Q(-a)}{Q(-b)}-A(\alpha ,\gamma )e^{-(b^2-a^2)/4}=0,$$

(27)

which, upon using the approximation (23) becomes:

$$\underset{f_i\left(\lambda \right)}{\underbrace{\frac{1-e^{1.4a}}{a}\frac{b}{1-e^{1.4b}}-A(\alpha ,\gamma )e^{-3\left(b^2-a^2\right)/4}}}=0$$

(28)

Denoting the function in the left-hand side of (28) by $f_i\left(\lambda \right)$, and initializing Newton’s iteration with the value ${\lambda }_{a,0}$ from (24) we obtain

$$f_i\left({\lambda }_{a,0}\right)=\frac{1.4b\left({\lambda }_{a,0}\right)}{e^{1.4b\left({\lambda }_{a,0}\right)}-1}-A(\alpha ,\gamma )e^{-0.75b^2\left({\lambda }_{a,0}\right)}$$

(29)

$$\begin{array}{ccc}\hfill f_i^{\prime }\left({\lambda }_{a,0}\right)& =& \frac{0.98}{\sqrt{2N}{\sigma }_n^2}\frac{b\left({\lambda }_{a,0}\right)}{e^{1.4b\left({\lambda }_{a,0}\right)}-1}-\frac{1.4}{\sqrt{2N}{\sigma }_n^2\left(1+\gamma \right)}·\left[\frac{1-e^{1.4b\left({\lambda }_{a,0}\right)}+1.4b\left({\lambda }_{a,0}\right)e^{1.4b\left({\lambda }_{a,0}\right)}}{{\left(e^{1.4b\left({\lambda }_{a,0}\right)}-1\right)}^2}\right]\hfill \\ & +& \frac{3A(\alpha ,\gamma )}{2}{e}^{-3b^2\left({\lambda }_{a,0}\right)/4}b\left({\lambda }_{a,0}\right)\frac{1}{\sqrt{2N}{\sigma }_n^2\left(1+\gamma \right)}\hfill \end{array}$$

(30)

The derivative $f_i^{\prime }\left({\lambda }_{a,0}\right)$ is given by the expression (30), which, together with (29), can be used to start Newton’s iterations (22) to obtain a numerical solution for the optimal sensing threshold ${\lambda }_{opt}^{\left(i\right)}$.

4.2. Case $\left(ii\right)a>0,b\le 0$

Following a similar approach as in the previous case but in which the symmetry property of the Q-function is used only for the $Q\left(b\right)$ term, Equation (21) is rewritten in this case as:

$$f\left(\lambda \right)=\frac{1-Q\left(a\right)}{Q(-b)}-A(\alpha ,\gamma )e^{-(b^2-a^2)/4}=0,$$

(31)

for which approximation (23) implies that

$$\begin{array}{c}\hfill \stackrel{f_{ii}\left(\lambda \right)}{\overbrace{\frac{b}{e^{1.4b}-1}·\left[1.135\sqrt{2\pi }-\frac{\left(1-e^{-1.4a}\right)e^{-a^2/2}}{a}\right]-A(\alpha ,\gamma )e^{-(3b^2-a^2)/4}}}=0\end{array}$$

(32)

Denoting the function in the left-hand side of (32) by $f_{ii}\left(\lambda \right)$, and initializing Newton’s iteration with the value ${\lambda }_{b,0}$ from (25) we obtain

$$\begin{array}{ccc}\hfill f_{ii}\left({\lambda }_{b,0}\right)& =& \frac{1.135\sqrt{2\pi }}{1.4}-\left[1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right]\frac{e^{-0.5a^2\left({\lambda }_{b,0}\right)}}{1.4a\left({\lambda }_{b,0}\right)}-A(\alpha ,\gamma )e^{0.25a^2\left({\lambda }_{b,0}\right)}\hfill \end{array}$$

(33)

The derivative $f_{ii}^{\prime }\left({\lambda }_{b,0}\right)$ is given by expression (34), which, together with (33), can now be used to start Newton’s iterations (22) to obtain a numerical solution for the optimal sensing threshold ${\lambda }_{opt}^{\left(ii\right)}$.

$$\begin{array}{ccc}\hfill f_{ii}^{\prime }\left({\lambda }_{b,0}\right)& =& -\frac{1}{2\sqrt{2N}{\sigma }_n^2\left(1+\gamma \right)}·\left[1.135\sqrt{2\pi }-\frac{\left(1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right)e^{-\frac{a^2\left({\lambda }_{b,0}\right)}{2}}}{a\left({\lambda }_{b,0}\right)}\right]\hfill \\ & +& \frac{\left(1-a^2\left({\lambda }_{b,0}\right)\right)\left(1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right)e^{-\frac{a^2\left({\lambda }_{b,0}\right)}{2}}}{1.4\sqrt{2N}{\sigma }_n^2{a}^2\left({\lambda }_{b,0}\right)}-\frac{e^{-1.4a\left({\lambda }_{b,0}\right)}{e}^{-\frac{a^2\left({\lambda }_{b,0}\right)}{2}}}{\sqrt{2N}{\sigma }_n^{2}a\left({\lambda }_{b,0}\right)}\hfill \\ & -& \frac{A(\alpha ,\gamma )a\left({\lambda }_{b,0}\right)}{2\sqrt{2N}{\sigma }_n^2}{e}^{\frac{a^2\left({\lambda }_{b,0}\right)}{4}}\hfill \end{array}$$

(34)

In case $\left(ii\right)$, we have used ${\lambda }_{b,0}$ as the initial value, but we can obtain similar expressions for $f_{ii}\left(\lambda \right)$ and $f_{ii}^{\prime }\left(\lambda \right)$ with ${\lambda }_{a,0}$ as the initial value, as in (33) and (34), respectively, with similar performance results for Newton’s method.

4.3. Case $\left(iii\right)a>b\ge 0$

In this case, the approximation from (23) can be used directly in (21) to rewrite it as expression (35). Denoting the function in the left-hand side of (35) by $f_{iii}\left(\lambda \right)$ and initializing it in this case with the same value ${\lambda }_{b,0}$ in (25) we obtain the expressions for $f_{iii}\left({\lambda }_{b,0}\right)$ and $f_{iii}^{\prime }\left({\lambda }_{b,0}\right)$ as shown in (36) and (37), respectively, which can be used to start Newton’s iterations (22) to obtain a numerical solution for the optimal sensing threshold ${\lambda }_{opt}^{\left(iii\right)}$.

$$\underset{f_{iii}\left(\lambda \right)}{\underbrace{\left[1.135\sqrt{2\pi }-\frac{\left(1-e^{-1.4a}\right)e^{-a^2/2}}{a}\right]·\frac{b}{1.135\sqrt{2\pi }b-\left(1-e^{-1.4b}\right)e^{-b^2/2}}-A(\alpha ,\gamma )e^{-(b^2-a^2)/4}}}=0$$

(35)

$$f_{iii}\left({\lambda }_{b,0}\right)=\frac{1}{1.135\sqrt{2\pi }-1.4}·\left[1.135\sqrt{2\pi }-\frac{\left(1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right)e^{-a^2\left({\lambda }_{b,0}\right)/2}}{a\left({\lambda }_{b,0}\right)}\right]-A(\alpha ,\gamma )e^{a^2\left({\lambda }_{b,0}\right)/4}$$

(36)

$$\begin{array}{ccc}{f}_{iii}^{\prime }\left({\lambda }_{b,0}\right)\hfill & =\hfill & \frac{e^{-a^2\left({\lambda }_{b,0}\right)/2}}{\left(1.135\sqrt{2\pi }-1.4\right)\sqrt{2N}{\sigma }_n^{2}a\left({\lambda }_{b,0}\right)}·\left\{\frac{\left[1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right]\left[a^2\left({\lambda }_{b,0}\right)+1\right]}{a\left({\lambda }_{b,0}\right)}-1.4e^{-1.4a\left({\lambda }_{b,0}\right)}\right\}\hfill \\ & -\hfill & \frac{1.96}{2{\left(1.135\sqrt{2\pi }-1.4\right)}^2\sqrt{2N}{\sigma }_n^2\left(1+\gamma \right)}·\left[1.135\sqrt{2\pi }-\frac{\left(1-e^{-1.4a\left({\lambda }_{b,0}\right)}\right)e^{-a^2\left({\lambda }_{b,0}\right)/2}}{a\left({\lambda }_{b,0}\right)}\right]\hfill \\ & -\hfill & \frac{A(\alpha ,\gamma )a\left({\lambda }_{b,0}\right)e^{a^2\left({\lambda }_{b,0}\right)/4}}{2\sqrt{2N}{\sigma }_n^2}\hfill \end{array}$$

(37)

4.4. One-Step Numerical Solution and the 3EED Algorithm with Adaptive Threshold

The rate of convergence for Newton’s iterations and the number of steps n needed to find a numerical solution ${\lambda }_{opt}^{\left(j\right)}$ to the optimal sensing threshold for each of the three cases outlined in the previous sections depend on the monotonicity of the corresponding $f_j\left(\lambda \right)$ function,

$${\lambda }_{opt}^{\left(j\right)}\cong {\lambda }_{n+1}={\lambda }_n-\frac{f_j\left({\lambda }_n\right)}{f_j^{\prime }\left({\lambda }_n\right)},j\in \{i,ii,iii\},$$

(38)

as well as on the initialization

$${\lambda }_0^{\left(j\right)}=\left\{\begin{array}{cc}\hfill {\lambda }_{a,0}& \mathrm{for} \mathrm{case}j=i\hfill \\ \hfill {\lambda }_{b,0}& \mathrm{for} \mathrm{cases}j=ii,iii\hfill \end{array}\right.$$

(39)

In order to improve convergence and minimize the overhead associated with finding the optimal sensing threshold, we propose to use the alternative function

$$F_j\left(\lambda \right)=ln\left[1+|f_j\left(\lambda \right)|\right]$$

(40)

having the derivative expressed as

$$F_j^{\prime }\left(\lambda \right)=\frac{f_j\left(\lambda \right)·f_j^{\prime }\left(\lambda \right)}{|f_j\left(\lambda \right)|+f_j^2\left(\lambda \right)},$$

(41)

which yields a good approximation for the optimal sensing threshold in only one iteration

$${\lambda }_{opt}^{\left(j\right)}\cong {\lambda }_1={\lambda }_0^{\left(j\right)}-\frac{F_j\left({\lambda }_0^{\left(j\right)}\right)}{F_j^{\prime }\left({\lambda }_0^{\left(j\right)}\right)},j\in \{i,ii,iii\}.$$

(42)

This approach was prompted by empirical observations that are discussed in Section 5.1 and takes advantage of the monotonicity of $f_j\left(\lambda \right)$, which, along with proper choice of the initial value ${\lambda }_0^{\left(j\right)}$, amount essentially to a linearization of the function that yields the optimal threshold around its root.

Noting that the functions $f_j\left(\lambda \right)$, $j\in \{i,ii,iii\}$, that approximate $f\left(\lambda \right)$ in (21) are all monotonically decreasing, and that they overlap at the edges of the definition intervals for a and b, that is

$$\begin{array}{ccc}\hfill f_i\left({\lambda }_{a,0}\right)& =& f_{ii}\left({\lambda }_{a,0}\right),\mathrm{and}\hfill \\ \hfill f_{ii}\left({\lambda }_{b,0}\right)& =& f_{iii}\left({\lambda }_{b,0}\right),\hfill \end{array}$$

(43)

the proposed approach identifies first which approximation case j is applicable ($\{i,ii,oriii\}$), and determines the one-step numerical solution for the optimal sensing threshold ${\lambda }_{opt}$ using (42). Once the sensing threshold is found, the 3EED Algorithm 1 is applied with the value ${\lambda }_{opt}$ to make a decision on the PU activity and channel availability. The approach is formally stated as Algorithm 2.

Algorithm 2 The 3EED algorithm with adaptive threshold

Input: N, $\alpha$, ${\sigma }_n^2$, ${\sigma }_s^2$

Output: output result

1: Use (24) and (25) to obtain ${\lambda }_{a,0}$ and ${\lambda }_{b,0}$. 2: Use (32) and (33) to obtain $f_{ii}\left({\lambda }_{a,0}\right)$ and $f_{ii}\left({\lambda }_{b,0}\right)$ 3: if$f_{ii}\left({\lambda }_{a,0}\right)<0$ and $f_{ii}\left({\lambda }_{b,0}\right)<0$ then 4:     $j=i$ in (42) and ${\lambda }_{opt}={\lambda }_{opt}^{\left(i\right)}$ 5: else if$f_{ii}\left({\lambda }_{a,0}\right)>0$ and $f_{ii}\left({\lambda }_{b,0}\right)<0$ then 6:     $j=ii$ in (42) and ${\lambda }_{opt}={\lambda }_{opt}^{\left(ii\right)}$ 7: else if$f_{ii}\left({\lambda }_{a,0}\right)>0$ and $f_{ii}\left({\lambda }_{b,0}\right)>0$ then 8:     $j=iii$ in (42) and ${\lambda }_{opt}={\lambda }_{opt}^{\left(iii\right)}$ 9: Apply Algorithm 1 with sensing threshold ${\lambda }_{opt}$ 10: return$q_i$ output by Algorithm 1

5. Simulations and Numerical Results

In this section, we present numerical results obtained from simulations that support the proposed approach for one-step threshold calculation and illustrate the performance of the 3EED algorithm with an adaptive threshold. The algorithm is also compared to the adaptive CED algorithm in [22] as well as with the “Brute-force” adaptive 3EED algorithm [24] in which an exhaustive search is implemented in Matlab to obtain the value of ${\lambda }_{opt}$ instead of the proposed one-step approach.

The parameter values used in the simulations are: the number of samples in a sensing slot $N=65,537$, the SNR $\gamma$ is between $-25$ dB and $-15$ dB, and $T=500$ slots. The signal transmitted by the PU is implemented using binary phase shift keying (BPSK) modulation, and 2500 sensing slots were considered in each transmission sequence in all simulations.

5.1. Sensing Threshold Calculation

Figure 1 shows the values of the optimal sensing threshold ${\lambda }_{opt}$ that minimizes the DEP, as a function of $\alpha$, the spectrum utilization ratio, for SNR values of $-25$ and $-20$ dB, corresponding to the proposed approach discussed in Section 4, the “Brute-force” approach in [24], and the adaptive CED approach in [22]. From this figure we note that the values of the sensing threshold obtained using the proposed approach closely match those obtained using the “brute-force” approach [24] for values of $\alpha$ between $0.2$ and $0.7$, in particular for $\alpha$ around $0.5$, and support the application of the numerical approach presented in Section 4 for threshold calculation with minimal overhead. We also note that, according to studies performed on GSM channels [30] or in the industrial, scientific and medical (ISM) bands [31], this is the range of values for the spectrum utilization ratio $\alpha$ that is of practical interest for SU access to licensed spectrum.

5.2. Decision Error Probability

Figure 2 shows the values of the DEP $P_e$ as a function of the SNR $\gamma$ for different values of the spectrum utilization ratio $\alpha$, for the proposed 3EED algorithm with an adaptive threshold as well as for the adaptive CED algorithm in [22]. Based on the plots shown in Figure 2 we first note that the values of the DEP obtained through Monte Carlo simulations closely match the analytical values of the DEP implied by (6) for the proposed 3EED algorithm with an adaptive threshold or by the corresponding DEP expression in [22], for all SNR and $\alpha$ values. We also note that the proposed 3EED algorithm with an adaptive threshold outperforms the adaptive CED algorithm in cases of practical interest corresponding to spectrum utilization ratios above $0.2$ and below $0.7$, with maximum gain of the proposed algorithm in terms of the DEP of about 1 dB achieved for $\alpha$ around $0.5$. Furthermore, for both the proposed 3EED algorithm with an adaptive threshold and the adaptive CED algorithm, the value of the DEP becomes less sensitive to changes in the spectrum utilization ratio $\alpha$ as the SNR values increase.

Figure 3 shows the dependence of the DEP on the spectrum utilization ratio $\alpha$ for two different SNR values $\gamma =-23$ dB and $-20$ dB. As it can be observed, the DEP decreases with increasing SNR for both the proposed 3EED algorithm with an adaptive threshold and the adaptive CED algorithm. Furthermore, the proposed algorithm outperforms the adaptive CED one for all values of the spectrum utilization ratio $\alpha$. We note that, for each value of $\alpha$, the corresponding DEP values are minimum since the optimal sensing threshold is used for that $\alpha$, as discussed in Section 3.

6. Conclusions

A new ED algorithm for spectrum sensing in CR systems is presented in the paper. The proposed algorithm employs an adaptive sensing threshold that is optimized to minimize the DEP and has minimal overhead, since the value of the optimal threshold is found through a one-step iterative method.

Numerical results obtained from simulations are presented to illustrate the performance of the proposed algorithm and to compare it to the adaptive CED algorithm. The results show that the proposed algorithm outperforms the CED algorithm for spectrum sensing, resulting in lower values for the DEP for all values of the spectral utilization ratio that are of practical interest for CR systems providing SU access to licensed spectrum.

We intend to implement and validate the proposed algorithm using SDR platforms from the USRP family. In this context, we will also focus on the optimization of the expressions used in the algorithm to estimate the decision threshold for minimizing the computational complexity and the sensing time.