Non-Cooperative Spectrum Sensing Strategy based on Recurrence Quantification Analysis in the Context of the Cognitive Radio

1. Introduction

The need to make better use of the radio spectrum is leading to the development of new spectrum access strategies. Among these strategies, the opportunistic spectrum access based on the Cognitive Radio concepts allows the sharing of a spectral bandwidth between two categories of users: Primary User “PU” and Secondary User “SU”. PU holds the license to exploit the bandwidth and SU is an opportunistic user willing to use the channel when the PU is idle. One of the most crucial challenge for SU is the identification of a free bandwidth by conducting a spectrum sensing [1,2]. Many reliable spectrum sensing methods have been eveloped to help the SU to limit their interference to PU’s transmission [3,4,5,6,7,8,9,10,11]. Among spectrum sensing approaches, we can mention the Waveform Detection (WFD) [12], the Cyclostationary Features based Detection (CFD)[13] and the Energy based Detection (ED) [14,15]. One of the most reliable methods, WFD, requires a prior knowledge on the PU’s signal characteristics. Based on the cyclic spectrum estimation, the CFD requires a relatively high computational cost for a high frequency resolution. The ED is the simplest detection method, but it is unable to distinguish a communication signal from an energetic noise, when the noise is not a weak sense stationary stochastic process or the Signal to Noise Ratio (SNR) is very low. Recently, spectrum sensing algorirthms, based on the promising concept of machine or deep learning, have been proposed [16,17,18,19]. However, these algorithms do not perform well in a non-cooperative context and at a low SNR (SNR ≤ -3 dB) and require a huge database to be optimized. To overcome these issues, we developped a blind strategy based on the Recurrence Quantification Analysis (RQA) of the received signal [20]. The quantification analysis of this recurrence can be used to find out some intrinsic features of a dynamic system, such as: hidden periodicities, stationarity features or linearity properties. Due to modulation standards, transmitted signals may contain hidden periodicities. Using this fact, we use Recurrence Quantification Analysis (RQA) tools to detect if the bandwidth allocated to PU is available or not. The main RQA tool used to quantify the recurrence level is the Recurrence Rate (RR) considered as the probability of having recurring states in a signal. In a recent work, we have proposed the RR based Detector (RRD) [20]. However, RRD is very sensitive to the SNR and depends on the choice of a recurrence threshold. To overcome the RRD limits, we propose in this paper an efficient algorithm called Recurrence Analysis based Detector (RAD). RAD exploits the similitude of distances among various states of the signal in a multidimensional space. This similitude of distances is evaluated by a square symmetrical matrix named the distance matrix. Using symmetrical property, we only exploit the upper triangular part of this matrix in order to reduce considerably the computational cost of the RAD. Then, we show that for a White Gaussian Noise (WGN), the coefficients of the first top diagonal of the distance matrix becomes a representative sample of all other coefficients. This is not the case for a communication signal even with a small SNR. This new approach can detect a communication signal in a very low SNR. We have analytically established the probabilities of detection $P_d$ and false alarm $P_{fa}$. Through Monte Carlo simulations, we studied the Receiver Operating Characteristic (ROC) curves of RAD. The theoretical and experimental results show the ability of RAD to detect the presence of a communication signal as soon as the SNR is greater than $-12$ dB with a very low probability of false alarm.

The rest of this paper is organized as follows: Section 2 presents the problem of spectrum sensing and our motivation for RQA. Section 3 deals with the concepts of RQA and the state of the art in the determination of embedding parameters in order to exhibit the hidden recurrences. Section 4 and Section 5 present the Recurrence Analysis based detection model and its theoretical and experimental performance. The last section contains the conclusion and perspectives.

2. Spectrum Sensing Problem

The radio spectrum is a limited natural ressource. Many techniques, such as cooperative communication systems and heterogeneous networks, have been developed to deal with this scarcity of the radio spectrum [21]. However, none of them can meet the strong demand for radio spectrum. The cognitive radio introduced by Mitola [22,23] is a promising solution allowing a dynamic access to the radio spectrum [1]. Dynamic spectrum access (DSA) is defined as a technique by which the operating spectrum of a radio network can be selected dynamically from the available spectrum [21]. The DSA allows the SU to exploit the holes in the spectrum dedicated to PU. The great challenge of DSA for SU remains the spectrum sensing stage during which the SU should detect the presence of PU on a given bandwidth.

2.1. Spectrum Sensing as a Statistical Decision

PU’s signal detection by the SU can be modeled as a binary hypothesis testing problem, given as [21]:

• Hypothesis ${\mathcal{H}}_0$: PU’s signal is absent
• Hypothesis ${\mathcal{H}}_1$: PU’s signal is present

Let $y\left(n\right)$ be the signal observed by the Secondary User (SU):

$$y\left(n\right)=hs\left(n\right)+b\left(n\right)$$

(1)

with h the channel gain and $b\left(n\right)$ the channel noise. The noise samples $b_k$ are assumed to be independent and identically distributed (i.i.d). The noise $b\left(n\right)$ is a Complex Circular Gaussian variable $\left(E\left[b\left(n\right)\right]=0;E\left[{\left|b\left(n\right)\right|}^2\right]={\sigma }_b^2\right);E\left[.\right]$ denotes the mathematical expectation.

The Test Statistic T, should be compared to a predetermined threshold $\lambda$ for decision making. In this case, the probability of false alarm $P_{fa}$ and the probability (Pr) of detection $P_d$ [24] are defined as follows:

$$\begin{array}{c}\hfill P_{fa}\triangleq Pr\{T\ge \lambda |{\mathcal{H}}_0\}\end{array}$$

(2)

$$\begin{array}{c}\hfill P_d\triangleq Pr\{T>\lambda |{\mathcal{H}}_1\}\end{array}$$

(3)

where ≜ stands for equal by definition.

For a good performance, $P_d$ should be increased as much as possible while keeping $P_{fa}$ under a small preselected value. In the non-cooperative spectrum sensing context, the detection algorithms should be able to find out intrinsic features of the communication signal to enhance a spectrum sensing task. One of the main features can be the recurrence of internal states in the communication signals. Hereinafter, we develop a robust algorithm based on RQA.

2.2. RQA Benefits

The RQA enjoys several specific properties:

1.

RQA is based on the chaos theory and normally used to extract the hidden recurring states of a dynamic system. The various part of a transmission chain such as modulation, filtering, coding, multiplexing, etc, generate hidden recurring states in the communication signals. Therefore, RQA can help detecting the presence of PU’s signal on a desired bandwidth.

2.

In a previous work [20], we showed that RQA is a promising tool for the spectrum sensing task. Indeed, in a non-cooperative context, we proposed the Recurrence Rate based detection model (RRD) and that previous algorithm was able to detect the presence of PU’s signals with SNR$\ge -5$ dB.

3.

During the detection procedure, a spectrum sensing algorithm based on RQA does not require the estimation of the noise variance as required by some spectrum sensing algorithms such as ED; which is a great advantage.

4.

RQA can help detecting a communication signal in a very low SNR; and contrary to ED, Recurrence Analysis can distinguish a noisy communication signal from a high energy noise.

5.

RQA does not need a high computational cost like CFD. In a non-cooperative spectrum sensing, RQA is more robust compared to widely used ED or CFD.

3. Recurrence Quantification Analysis

Recurrence is a fundamental feature of dynamic systems, which can be exploited to study the behavior of these systems and to discover their intrinsic properties [20,25,26]. To analyze hidden recurrences in dynamic and non-linear systems, an important concept in the chaos theory, RQA can be used [27,28,29,30]. During the evolution of a dynamic system, some internal states can be quasi periodically repeated. In the phase space, the successive states of a dynamic system, form trajectories that characterize the temporal evolution of the system. Knowing the temporal evolution equation, any state of the system can be precisely determined at any time. Unfortunately in real world situations, this equation is unknown; instead, we dispose only a sequence of scalar measurements as time series [31,32]. From these time series, we should reconstruct the phase space. The gold standard method for the phase space reconstruction is the delay-coordinate embedding [31]. The objective of a delay-coordinate embedding method is to use the delayed versions of $y\left(n\right)$ to form a multi-dimensional observable ${\mathbf{x}}_k$ called state vector or state of $y\left(n\right)$. The state vector ${\mathbf{x}}_k$ at moment k is defined as follows:

$${\mathbf{x}}_k={[y_k,y_{k+\tau },\dots ,y_{k+(m-1)\tau }]}^T$$

(4)

where $y_k$ denotes the sample of y at the moment k, $\tau$ is the time delay and m, the embedding dimension. The main challenge for the delay-coordinate embedding method is the reliable estimation of $\tau$ and m.

3.1. Estimation of the embedding parameters

3.1.1. Time Delay

Theoretically, the time delay $\tau$ can be almost arbitrarily chosen if the observation $y\left(n\right)$ is noise free with an infinite number of samples. However, these conditions cannot be satisfied in real applications because of noise which can generate statistical dependence among the state vectors ${\mathbf{x}}_k$; Therefore, the time delay $\tau$ has to be wisely chosen in order to reduce this statistical dependence [20,31,33]. To determine the optimal time delay ${\tau }_0$, one can use the autocorrelation function, Average Mutual Information (AMI) or phase portrait approach. The most appropriate method is the AMI because it measures the general dependence between two random variables [34]. Therefore, it could provide a better criterion for the optimal time delay ${\tau }_0$. The concept of AMI consists, first of all, in estimating the mutual information $I\left(\tau \right)$ between $y\left(n\right)$ and its delay version $y(n-\tau )$ by varying the value of $\tau$ from 0 to N ; N denotes the number of samples contained in $y\left(n\right)$. After that, the optimal ${\tau }_0$ is chosen as the first $\tau$ value that minimizes $I\left(\tau \right)$ [20,35,36]. Based on equation (1), we can conclude that $I\left(\tau \right)$ and indirectly ${\tau }_0$, depends on three parameters: the SNR value, the sampling rate and the modulation scheme. Consequently, use AMI method to determine ${\tau }_0$ is not so suitable in the context of non-cooperative spectrum sensing algorithm.

3.1.2. Embedding Dimension

The optimal embedding dimension $m_0$ for an observed signal $y\left(n\right)$ is the minimum dimension for the state vectors ${\mathbf{x}}_k$ to give a reliable reconstruction of $y\left(n\right)$ phase space [30,32]. From the literature, many approaches to estimate $m_0$ have been developed [33,37,38]. The most used approach is based on the False Nearest Neighbors (FNN) method [38]. According to the principle of FNN, any two true neighboring points in the $m-dimensional$ reconstructed phase space must remain neighbors in the $(m+1)-dimensional$ reconstructed phase space. Otherwise, they are called false neighbors. A perfect embedding means all neighboring points should be true neighbors [20,37,38]. Inspired by the FNN-based algorithm proposed in[37], we identify $m_0$ by using the distance ratio, $a(i,m)$, defined as follows [37]:

$$a(k,m)=\frac{{∥{\mathbf{x}}_k^{(m+1)}-{\mathbf{x}}_l^{(m+1)}∥}_{\infty }}{{∥{\mathbf{x}}_k^{\left(m\right)}-{\mathbf{x}}_l^{\left(m\right)}∥}_{\infty }};$$

(5)

where $k,l=\in \{1,2,\cdots ,N-m\tau \}$; ${\mathbf{x}}_k^{\left(m\right)}$ are the state vectors from the $m-dimensional$ phase space; and ${∥.∥}_{\infty }$ is $L-infinity$ norm. The major drawback of using the distance ratio $a(k,m)$ is its sensitivity to the reference state vector ${\mathbf{x}}_k$. To overcome this drawback, we consider the average value $\overline{a}(m,\tau )$ instead of $a(k,m)$:

$$\overline{a}(m,\tau )=\frac{1}{N-m\tau }\sum _{k=1}^{N-m\tau }a(k,m)$$

(6)

To investigate the dependence of $\overline{a}(m,\tau )$ on m, the coefficient of proximity $r(m,\tau )$ is defined as follows:

$$r(m,\tau )\triangleq \frac{\overline{a}(m+1,\tau )}{\overline{a}(m,\tau )}$$

(7)

When m becomes greater than a limit value $m_{lim}$, the proximity coefficient $r(m,\tau )$ converges to one. Hence, $m_0=m_{lim}+1$ becomes the minimum embedding dimension [37]. Based on, equation (7), we can observe that $r(m,\tau )$ depends of the time delay value $\tau$. This assertion is corroborated by simulation results. Figure 1 illustrates, for example, the evolution of $r(m,\tau )$ for two extreme values of $\tau$; for $\tau =2$, the optimal value of $m_0=6$; whereas $\tau =100$, we obtain $m_0=9$. As m depends on the choice of $\tau$, we propose hereinafter an optimization strategy to find simultaneously the optimal values of m and $\tau$.

3.1.3. Optimal Values of m and $\tau$

$r(m,\tau )$ is a bivariate function whose expression leads to

$$\underset{(m,\tau )\to (m_0,{\tau }_0)}{lim}r(m,\tau )=1$$

(8)

So that, we define the cost function $f(m,\tau )$:

$$f(m,\tau )=\overline{a}(m,\tau )-\overline{a}(m+1,\tau )$$

(9)

The optimal values $m_0$ and ${\tau }_0$ become:

$$(m_0,{\tau }_0)=\underset{(m,\tau )\in {\mathbb{N}}^{*}}{argmin}\left\{f(m,\tau )\right\}$$

(10)

Numerical results show that values of $(m_0,{\tau }_0)$ obtained from (10) become independent from modulation schemes, SNR and the number of samples per symbol $N_s$. Figure 2 illustrates the cost function f$(m,\tau )$ with respect to $(m,\tau )$ for a noisy 16-QAM signal. We can notice that:

$$\begin{array}{c}\hfill \forall \tau \ge 6,m\ge 1,\mathbf{f}(m,\tau )\to 0\end{array}$$

(11)

The number K of state vectors is given by:

$$K=N-(m-1)\tau \ge K_{min}$$

(12)

Based on Takens’ theorem1 and the minimum number $K_{min}$ of state vectors for reliable detection, we conclude that $m_0$ and ${\tau }_0$ should be chosen in the following ranges:

$$\begin{array}{c}\hfill \left(m_0,{\tau }_0\right)\in \left[3;\frac{N-K_{min}}{{\tau }_0}\right]\times \left[6;\frac{N-K_{min}}{2}\right]\end{array}$$

(13)

However, we notice that in practice $K_{min}\ll N$. Consequently, we can:

$$\begin{array}{c}\hfill \left(m_0,{\tau }_0\right)\in \left[3;\frac{N}{{\tau }_0}\right]\times \left[6;\frac{N}{2}\right]\end{array}$$

(14)

The recurring states in the phase space can be highlighted by the Recurrence Plot.

3.2. Recurrence Plot

The Recurrence Plot (RP) illustrates recurrences contained in a signal. RP is based on the recurrence matrix ${\mathbf{R}}_{i,j}^{(\epsilon ,m)}$ [25,31]:

$${\mathbf{R}}_{i,j}^{(\epsilon ,m)}=\Theta \left\{\epsilon -{∥{\mathbf{x}}_i-{\mathbf{x}}_j∥}_2\right\}$$

(15)

where $i,j\in \{1,\cdots ,K\}$ and $K=N-(m-1)\tau$ denotes the number of reconstructed state vectors ${\mathbf{x}}_k$; $\epsilon$ is the recurrence threshold, $\Theta (.)$ represents Heavisides’ step function, and ${∥.∥}_2$ is the $L_2$ or Frobenius norm. $d_{ij}={∥{\mathbf{x}}_i-{\mathbf{x}}_j∥}_2$ are the coefficients of the distance matrix M.

According to equation (15), the states ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ are recurring states iff $∥{\mathbf{x}}_i-{\mathbf{x}}_j∥<\epsilon$. In RP, a recurrence is represented by a black dot. If the parameters m, $\tau$ and $\epsilon$ are optimal, the RP presents some intrinsic patterns of the system. For example, Figure 3, Figure 4 and Figure 5 represent respectively the RP of 16-QAM signal, GMSK signal and a random zero mean White Gaussian Noise (WGN) with a variance ${\sigma }_b^2$. We can notice that the patterns in the RP of 16-QAM signal are different from those of GMSK signal, whereas the RP of the WGN has no particular pattern. Heuristically, we can set the recurrence threshold is $\epsilon =0.5{\sigma }_y$ where ${\sigma }_y^2$ denotes the variance of the observed signal $y\left(n\right)$; the embedding parameters are established to $m=16$ and $\tau =6$.

As the visual analysis of RP is not objective, Zbilut and Webber introduced a procedure to quantify RP structures [40,41]. In the literature, one have five classical tools to perform the RQA. Some RQA tools are based on the recurrence density, while others use the line structures of the RP [31]. The classical RQA tools are developed and extended in [25,31]. The main challenge with classical or extended RQA tools, is the choice of an optimal recurrence threshold $\epsilon$ of equation (15). A small value of $\epsilon$ does not reveal noticeable occurrences, while a large value of $\epsilon$ may lead to the appearance of neighbors for most of existing points and cause false occurrences [25,42]. The choice of $\epsilon$ is a delicate issue affecting the reliability of the spectrum sensing based on classical or extended RQA tools [20]. To overcome this issue, we propose in the following sections a new algorithm called the Recurrence Analysis based Detector (RAD) which only exploits the distance matrix. This distance matrix does not depend on the recurrence threshold $\epsilon$; and nor does proposed detection model RAD.

4. Recurrence Analysis Based Detector

Recurrence Rate (RR) is an essential tool of RQA. Our previous works based on the RQA detection model [20] show that RR-based detection model (RRD) suffers from major shortcomings, such as:

1.

It cannot detect the presence of a communication signal when SNR $\le -5$ dB.

2.

It is very sensitive to the recurrence threshold $\epsilon$.

3.

The computational cost is relatively high

4.

The performance of RRD is sensitive to the type of modulations for a communication signal.

To overcome the above limitations, we develop hereinafter a new detection model, the Recurrence Analysis Detection model (RAD), which is able to operate in very low SNR conditions. In order to reduce drastically the computational cost and avoid the delicate issue of the recurrence threshold $\epsilon$, the RAD only uses the distances $d_{ij}={∥{\mathbf{x}}_i-{\mathbf{x}}_j∥}_2$ of the upper triangular part of the distance matrix.

4.1. Detection Model

Usually, the RQA is performed on the entire distance matrix $D={\left(d_{ij}\right)}_{1\le i,j\le K}$ of different state vectors. As D is a symmetrical matrix, in order to reduce the computational cost, we use the upper triangular part of D, without the main diagonal defined by its general coefficients:

$$D_{ij}=d_{ij};\forall 1\le i<j\le K$$

(16)

We denote by ${\mathbf{u}}_d$ the coefficients of the first upper diagonal of D and by ${\mathbf{u}}_{\Delta }$, the other coefficients of the upper triangular part of D, without the main diagonal. ${\mathbf{u}}_d$ and ${\mathbf{u}}_{\Delta }$ verify two great properties:

1.

For a WGN, ${\mathbf{u}}_d$ and ${\mathbf{u}}_{\Delta }$ have the same Probability Density Function (PDF); which is not the case for a noisy communication signal. Hence, to detect the presence of a communication signal, we check if ${\mathbf{u}}_d$ is representative of ${\mathbf{u}}_{\Delta }$. To this purpose, we use a statistical test of conformity to evaluate this representativeness [43].

2.

${\mathbf{u}}_d$ of a WGN or a communication signal have the same PDF. This remark allows us to design a detector free of noise variance estimation.

To provide a better understanding of this notion of representativeness, we present the histogram of ${\mathbf{u}}_d$ and ${\mathbf{u}}_{\Delta }$ of WGN and communication signal respectively in Figure 6 and Figure 7, and we define the confidence interval ${\Phi }_{\Delta }$ as follows:

$${\Phi }_{\Delta }=\left[{\overline{\mathbf{u}}}_{\Delta }-\lambda ;{\overline{\mathbf{u}}}_{\Delta }+\lambda \right]$$

(17)

where ${\overline{\mathbf{u}}}_{\Delta }$ is the average value of ${\mathbf{u}}_{\Delta }$ and $\lambda$ is the predetermined detection threshold for RAD. Further details about $\lambda$ can be founded in subSection 4.2.4.

On the basis of the histogram of distance and by using Kolmogov-Smirnov test [44], for a WGN only, we note that ${\mathbf{u}}_d$ and ${\mathbf{u}}_{\Delta }$ can be approximated by the same PDF (see Figure 6). Consequently, by using the statistical test of conformity based on the estimation average value, we show that ${\mathbf{u}}_d$ is representative of ${\mathbf{u}}_{\Delta }$. Indeed, ${\overline{\mathbf{u}}}_d\in {\Phi }_{\Delta }$; ${\overline{\mathbf{u}}}_d$ denotes the average value of ${\mathbf{u}}_d$. Contrariwise in Figure 7, which gives the histogram of distance for a communication signal buried in the WGN with a SNR$=0$ dB, we can notice that ${\overline{\mathbf{u}}}_d\notin {\Phi }_{\Delta }$; so ${\mathbf{u}}_d$ is not representative of ${\mathbf{u}}_{\Delta }$.

Using the statistical test of conformity, we propose the following test statistic T:

$$T=|{\overline{\mathbf{u}}}_{\Delta }-{\overline{\mathbf{u}}}_d|\stackrel{{\mathcal{H}}_1}{\underset{{\mathcal{H}}_0}{\gtrless }}\lambda$$

(18)

In the course of our work, we will establish the analytical expressions of the probability of detection $P_d$ and the probability of false alarm $P_{fa}$ in order to compare the theoretical and experimental results.

4.2. Analytical Expression of the Probability of False Alarm

The probability of false alarm $P_{fa}$ can be expressed from the PDF of T under ${\mathcal{H}}_0$. In order to find this PDF, one should evaluate the PDF of ${\overline{\mathbf{u}}}_{\Delta }$ and ${\overline{\mathbf{u}}}_d$, see equation (18).

4.2.1. The PDF of ${\overline{\mathbf{u}}}_{\Delta }$ under hypothesis ${\mathcal{H}}_0$

Under ${\mathcal{H}}_0$, the average value ${\overline{\mathbf{u}}}_{\Delta }$ of ${\mathbf{u}}_{\Delta }$ is defined as follows:

$${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}=\frac{1}{K_{\Delta }}\sum _{k=1}^{K_{\Delta }}{u}_{\Delta /{\mathcal{H}}_0}^{\left(k\right)}$$

(19)

where $K_{\Delta }=(K-1)(K/2-1)$ and $u_{\Delta /{\mathcal{H}}_0}^{\left(k\right)}\in \{d_{ij};1\le i\le K-2;i+2\le j\le K\}$.

Based on the central limit theorem for independent random variables and for a large K, ${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}$ is asymptotically normally distributed with mean ${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}$ and variance ${\sigma }_{\Delta /{\mathcal{H}}_0}^2$:

$${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}\mapsto \mathcal{N}\left({\mu }_{\Delta /{\mathcal{H}}_0},{\sigma }_{\Delta /{\mathcal{H}}_0}^2\right)$$

(20)

where

$$\begin{array}{c}\hfill {\mu }_{\Delta /{\mathcal{H}}_0}=\mathbb{E}\left[{\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}\right]\end{array}$$

(21)

$$\begin{array}{c}\hfill {\sigma }_{\Delta }^2=var\left[{\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}\right]\end{array}$$

(22)

To calculate ${\mu }_{\Delta /{\mathcal{H}}_0}$ and ${\sigma }_{\Delta /{\mathcal{H}}_0}^2$, we can use the PDF $f_{{\mathbf{u}}_{\Delta /{\mathcal{H}}_0}}$ (u) of ${\mathbf{u}}_{\Delta /{\mathcal{H}}_0}$. An outcome $u_{\Delta /{\mathcal{H}}_0}$ of ${\mathbf{u}}_{\Delta /{\mathcal{H}}_0}$ is given by:

$$u_{\Delta /{\mathcal{H}}_0}=d_{ij}=\sqrt{\sum _{k=0}^{m-1}{\left(b_{i+k\tau }-b_{j+k\tau }\right)}^2}$$

(23)

As $b\left(n\right)\mapsto \mathcal{N}\left(0,{\sigma }_b^2\right)$,then $b_{i+k\tau }-b_{j+k\tau }\mapsto \mathcal{N}\left(0,2{\sigma }_b^2\right)$ and $v_k=\left(\frac{b_{i+k\tau }-b_{j+k\tau }}{{\sigma }_b\sqrt{2}}\right)\mapsto \mathcal{N}\left(0,1\right)$.By setting: $Z=\sqrt{{\sum }_{k=0}^{m-1}{v}_k^2}$ , we can conclude that Z follows a Chi distribution $\chi \left(m\right)$ with m degrees of freedom [45,46]:

$$Z=\frac{{\mathbf{u}}_{\Delta /{\mathcal{H}}_0}}{{\sigma }_b\sqrt{2}}\mapsto \chi \left(m\right)$$

(24)

The expectation value $\mathbb{E}\left[Z\right]$ and variance $var\left[Z\right]$ of Z are given as follows [45,46]:

$$\begin{array}{c}\hfill \mathbb{E}\left[Z\right]={\mu }_Z=\sqrt{2}\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\end{array}$$

(25)

$$\begin{array}{c}\hfill var\left[Z\right]=m-{\mu }_Z^2\end{array}$$

(26)

where $\Gamma \left(m\right)={\int }_0^{+\infty }{x}^{m-1}{e}^{-x}dx$ is the Gamma function.

From equation (24), we show that the PDF $f_{{\mathbf{u}}_{\Delta /{\mathcal{H}}_0}}$ of ${\mathbf{u}}_{\Delta /{\mathcal{H}}_0}$ is defined as follows [47]:

$$f_{{\mathbf{u}}_{\Delta /{\mathcal{H}}_0}}\left(u\right)=\frac{u^{m-1}}{2^{m-1}{\sigma }_b^m\Gamma \left(\frac{m}{2}\right)}exp\left[-\frac{u^2}{4{\sigma }_b^2}\right]$$

(27)

The mean ${\mu }_{\Delta /{\mathcal{H}}_0}$ and variance ${\sigma }_{\Delta /{\mathcal{H}}_0}^2$ of ${\overline{X}}_{\Delta /{\mathcal{H}}_0}$ can be determined by using equation (24) as follows:

$$\begin{array}{c}\hfill {\mu }_{\Delta /{\mathcal{H}}_0}={\sigma }_b\sqrt{2}\mathbb{E}\left[Z\right]\end{array}$$

(28)

$$\begin{array}{c}\hfill {\sigma }_{\Delta /{\mathcal{H}}_0}^2=\frac{2{\sigma }_b^2}{K_{\Delta }}var\left[Z\right]\end{array}$$

(29)

Finally, we have ${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}\mapsto \mathcal{N}\left({\mu }_{\Delta /{\mathcal{H}}_0},{\sigma }_{\Delta /{\mathcal{H}}_0}^2\right)$ with :

$$\begin{array}{c}\hfill {\mu }_{\Delta /{\mathcal{H}}_0}=2{\sigma }_b\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\end{array}$$

(30)

$$\begin{array}{c}\hfill {\sigma }_{\Delta /{\mathcal{H}}_0}^2=\frac{2{\sigma }_b^2}{K_{\Delta }}\left[m-2{\left(\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\right)}^2\right]\end{array}$$

(31)

4.2.2. The PDF of ${\overline{\mathbf{u}}}_d$ under hypothesis ${\mathcal{H}}_0$

Under ${\mathcal{H}}_0$, the average value ${\overline{\mathbf{u}}}_d$ of ${\mathbf{u}}_d$ is defined as follows:

$${\overline{\mathbf{u}}}_{d/{\mathcal{H}}_0}=\frac{1}{K_d}\sum _{k=1}^{K_d}{u}_{d_k/{\mathcal{H}}_0}^{\left(k\right)}$$

(32)

where $u_{d/{\mathcal{H}}_0}^{\left(k\right)}\in \left\{d_{i,i+1}\right\};1\le i\le K_d=K-1$.

Because of the structure of the upper triangular distance matrix, the coefficients of the first top diagonal $u_{d/{\mathcal{H}}_0}^{\left(k\right)}$ are strongly decorrelated. In addition, $K_d$ large number. Thus, we can approximate ${\overline{\mathbf{u}}}_d$ by a Gaussian variable by using the central limit theorem:

$${\overline{\mathbf{u}}}_{d/{\mathcal{H}}_0}\mapsto \mathcal{N}\left({\mu }_{d/{\mathcal{H}}_0},{\sigma }_{d/{\mathcal{H}}_0}^2\right)$$

(33)

Using a similar approach to the calculation of the expectation value ${\mu }_{\Delta /{\mathcal{H}}_0}$ and variance ${\sigma }_{\Delta /{\mathcal{H}}_0}^2$ of ${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}$, we obtain:

$$\begin{array}{c}\hfill {\mu }_{d/{\mathcal{H}}_0}=2{\sigma }_b\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\end{array}$$

(34)

$$\begin{array}{c}\hfill {\sigma }_{d/{\mathcal{H}}_0}^2=\frac{2{\sigma }_b^2}{K_d}\left[m-2{\left(\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\right)}^2\right]\end{array}$$

(35)

4.2.3. The PDF of T under ${\mathcal{H}}_0$

${\overline{\mathbf{u}}}_{d/{\mathcal{H}}_0}$ and ${\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}$ are Gaussian random variables. Based on equations (20), (30), (31), (33), (34) and (35) we end up with:

$$\begin{array}{c}\hfill {\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}-{\overline{\mathbf{u}}}_{d/{\mathcal{H}}_0}\mapsto \mathcal{N}\left(0,{\sigma }_0^2\right)\end{array}$$

(36)

where

$$\begin{array}{c}\hfill {\sigma }_0^2={\sigma }_{\Delta /{\mathcal{H}}_0}^2+{\sigma }_{d/{\mathcal{H}}_0}^2-2Cov\left({\overline{\mathbf{u}}}_{d/{\mathcal{H}}_0},{\overline{\mathbf{u}}}_{\Delta /{\mathcal{H}}_0}\right)\end{array}$$

(37)

Applying the random variable transformation theorem [47,48] on T from equation (18), we end up with $t\ge 0$:

$$\begin{array}{c}\hfill f_{T/{\mathcal{H}}_0}\left(t\right)=\frac{2}{{\sigma }_0\sqrt{2\pi }}exp\left[-\frac{t^2}{2{\sigma }_0^2}\right]\end{array}$$

(38)

4.2.4. Probability of False Alarm and Detection Threshold

Based on equation (38), the $P_{fa}$ is expressed as follows:

$$\begin{array}{c}\hfill P_{fa}={\int }_{\lambda }^{\infty }{f}_{T/{\mathcal{H}}_0}\left(t\right)dt=erfc\left(\frac{\lambda }{{\sigma }_0\sqrt{2}}\right)\end{array}$$

(39)

$erfc\left(x\right)$ denotes the Complementary Error Function2, and the detection threshold $\lambda$ is:

$$\begin{array}{c}\hfill \lambda ={\sigma }_0\sqrt{2}erf{c}^{-1}\left(P_{fa}\right)\end{array}$$

(40)

4.3. Analytical Expression of the Probability of Detection

Taking into account the presence of the communication signal $s\left(n\right)$ and keeping the same approach as under $H_0$, we demonstrate that:

$$\begin{array}{c}\hfill {\overline{\mathbf{u}}}_{\Delta /H_1}\mapsto \mathcal{N}\left({\mu }_{\Delta /H_1},{\sigma }_{\Delta /H_1}^2\right)\end{array}$$

(41)

$$\begin{array}{c}\hfill {\overline{\mathbf{u}}}_{d/H_1}\mapsto \mathcal{N}\left({\mu }_{d/H_1},{\sigma }_{d/H_1}^2\right)\end{array}$$

(42)

where:

$$\begin{array}{c}\hfill {\mu }_{\Delta /H_1}=2\sqrt{{\sigma }_s^2+{\sigma }_b^2}\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\end{array}$$

(43)

$$\begin{array}{c}\hfill {\sigma }_{\Delta /H_1}^2=\frac{2({\sigma }_s^2+{\sigma }_b^2)}{K_{\Delta }}\left[m-2{\left(\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\right)}^2\right]\end{array}$$

(44)

$$\begin{array}{c}\hfill {\mu }_{d/H_1}=2{\sigma }_b\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\end{array}$$

(45)

$$\begin{array}{c}\hfill {\sigma }_{d/H_1}^2=\frac{2{\sigma }_b^2}{K_d}\left[m-2{\left(\frac{\Gamma \left(\frac{m+1}{2}\right)}{\Gamma \left(\frac{m}{2}\right)}\right)}^2\right]\end{array}$$

(46)

We deduce the PDF of T under the hypothesis $H_1$ as follows:

$$\begin{array}{c}\hfill f_{T/H_1}\left(t\right)=\frac{1}{{\sigma }_1\sqrt{2\pi }}\left[e^{-{\left(\frac{t-{\mu }_1}{{\sigma }_1\sqrt{2}}\right)}^2}+e^{-{\left(\frac{t+{\mu }_1}{{\sigma }_1\sqrt{2}}\right)}^2}\right]\end{array}$$

(47)

with $t\ge 0$. ${\mu }_1$ and ${\sigma }_1^2$ are given by:

$$\begin{array}{c}\hfill {\mu }_1={\mu }_{\Delta /H_1}-{\mu }_{d/H_1},\end{array}$$

(48)

$$\begin{array}{c}\hfill {\sigma }_1^2={\sigma }_{\Delta /H_1}^2+{\sigma }_{d/H_1}^2-2Cov\left({\overline{\mathbf{u}}}_{\Delta /H_1},{\overline{\mathbf{u}}}_{d/H_1}\right)\end{array}$$

(49)

Consequently, the Probability of Detection $P_d$ becomes:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_d& ={\int }_{\lambda }^{+\infty }{f}_{T/H_1}\left(t\right)dt\hfill \\ & =1-\frac{1}{2}\left[erfc\left(\frac{-\lambda -{\mu }_1}{{\sigma }_1\sqrt{2}}\right)-erfc\left(\frac{\lambda -{\mu }_1}{{\sigma }_1\sqrt{2}}\right)\right]\hfill \end{array}\end{array}$$

(50)

5. Simulations Results

To evaluate the efficiency and the robustness of the proposed detection method, we generate the Receiver Operating Characteristic (ROC) curves using Monte Carlo simulations. The parameters defined in Table 1 are used with different kind of communication signals such as 64-QAM, 16-QAM, BPSK and 4-ASK signals.

In order to compare the theoretical and experimental performance of RAD, we have generated the performance curve $P_d$ versus SNR. Figure 8 shows these performance curves. The theoretical curves were generated by equation (50) and the experimental curves according to the Monte Carlo simulations. It can be clearly observed that the theoretical results match with simulated ones. From theoretical and experimental curves, the reliable detection is possible as soon as $SNR\ge -12$ dB with $P_d\ge 0.95$ and $P_{fa}=0.05$.

The Receiver Operating Characteristic curves for different $SNR$ values can be viewed on Figure 9. Here, the detector RAD is applied on a 16-QAM signal and the observation time is $31.25$ ms. For SNR $=-14$ dB and $P_{fa}=0.1$, RAD detects the presence of a communication signal with $P_d=0.76$; but for SNR $=-13$ dB, $P_d=0.9$. The detector proves itself powerful as soon as the SNR $\ge -12$ dB, since the detection probability $P_d=0.95$ with a very low value of $P_{fa}=0.05$; see Figure 9.

Another advantage of the proposed detector is its robustness against any type of the classical modulations such as QAM, PSK and ASK; see Figure 10, where the ROC curves for 4-ASK, BPSK and 64-QAM signal are almost identical.

The performance of RAD increases with the number of samples N. Figure 11 depicts this performance evolution based on N in a Gaussian channel where SNR $=-12$ dB. We notice that for $P_{fa}=0.1$, the RAD detects the communication signal with $P_d=0.5$ when $N=500$ samples, with $P_d=0.8$ when $N=1000$ samples and with $P_d=0.97$ when $N=2000$ samples.

By comparing RAD to two other blind spectrum sensing algorithms: Energy Detection (ED) [49,50] and Cyclostationary Features Detector (CFD) [51,52] and according to the results presented in Figure 12 and Figure 13, we notice the superiority of RAD to ED and CFD in a Gaussian channel. For example, for the $SNR=-12$ dB and $P_{fa}=0.05$, the RAD detects the signal with $P_d=0.95$, whereas the probability of detection for the ED is only $P_d=0.67$ (see Figure 12). Figure 13 shows that the RAD is able to detect the communication signal in very weak low SNR conditions. For $-12$ dB ⩽ SNR $⩽-10$ dB, RAD detects the presence of a communication signal with $P_d⩾0.95$ for $P_{fa}=0.05$. ED achieves RAD performance only when SNR $⩾-10$ dB and CFD works correctly when SNR$⩾-7$ dB with $P_d⩾0.97$.

After testing the performance of the detector in a Gaussian channel, we are now interested in the behaviour of the RAD in the multi-path channel. As a model, we use the model D of Rayleigh channel defined in [53] with the following parameters:

Table 2. Rayleigh Channel Model Features.

Model D of Rayleigh Channel
Entity	Parameters	Value
Features	Number of Path	6
	Doppler Frequency $f_{d_{max}}$	1,2 kHz

Table 2. Rayleigh Channel Model Features.

Model D of Rayleigh Channel
Entity	Parameters	Value
Features	Number of Path	6
	Doppler Frequency $f_{d_{max}}$	1,2 kHz

Table 3. Delay and Gain values for 6 paths [53].

Path Number	1	2	3	4	5	6
Delay (ns)	0	300	8900	12900	17100	20000
Gain (dB)	0	-2.5	-12.8	-10.0	-25.2	-16.0

Table 3. Delay and Gain values for 6 paths [53].

Path Number	1	2	3	4	5	6
Delay (ns)	0	300	8900	12900	17100	20000
Gain (dB)	0	-2.5	-12.8	-10.0	-25.2	-16.0

Figure 14 summarizes the performance of RAD in a noisy Rayleigh channel. The RAD performance remains almost unchanged. It detects a communication signal with $P_d\ge 0.92$ with $P_{fa}=0.05$ as soon as SNR $\ge -12$ dB.

6. Complexity Analysis of Recurrence Analysis Based Detector

Theoretical and experimental analyses show the superiority of RAD compare to Energy based Detector (ED) and Cyclostationary Feature based Detector (CFD) in low SNR scenario detection process. The complexity of the algorithms is measured through the number of complex multiplications that the algorithms has to perform for the calculation of the test statistics [54]. In this section, we provide the complexity analyses of ED, CFD and RAD.

6.1. Complexity Analysis of Energy Based Detector

Let $y\left(n\right)$ be the observed siganl with N samples $y_k$. The energy ${\mathcal{E}}_y$ of $y\left(n\right)$ is given by:

$${\mathcal{E}}_y=\frac{1}{N}\sum _{k=1}^N{y}_k^2$$

(51)

The ED detection process can be summarize through the Figure 15.

The complexity computational can be evaluated according to equation (51). N multiplication operations are required to perform $y_k^2$. Consequently, the computation complexity $C_{ED}$ becomes[54]:

$$C_{ED}=N$$

(52)

6.2. Complexity Analysis of Cyclostationary Feature Based Detector

In a blind context, the CFD is based on the reliable estimation of the cyclic spectrum [13,55,56]. The crest factor $F_c$ of the cyclic spectrum can be used as a decision statistic [13]:

$$F_c\left(\alpha \right)=\frac{{max}_{\alpha }D\left(\alpha \right)}{\sqrt{\frac{1}{2N+1}{\sum }_{k=1}^{2N+1}{D}^2\left(\alpha \right)}}\stackrel{{\mathcal{H}}_1}{\underset{{\mathcal{H}}_0}{\gtrless }}{\lambda }_{CFD}$$

(53)

where $\alpha$ stands for the cyclic frequency and $D\left(\alpha \right)$ denotes the Cyclic Domain Profile (CDP). The computation of $F_c$ requires 6 steps [13]. First, a Hamming sliding window $w\left(n\right)$ is used to apodize the observed signal $y\left(n\right)$. Then, we apply a first Fast Fourier Transform (FFT) and calculate the complex demodulates of the apodized signal. After that, we compute the spectral correlation of the different complex demodulates and apply a second FFT. At the end, we obtain the cyclic spectrum and we compute the crest factor $F_c$ of the cyclic spectrum.

Let N be the number of samples contained in the observed signal $y\left(n\right)$ and $N_p$ the number of samples contained in each apodized block of $y\left(n\right)$. $L=\frac{N_p}{4}$ denotes the decimation factor and $P=\frac{N}{L}$ is the number of apodized blocks from $y\left(n\right)$. The apodization is carried out with a Hamming window [13,55,57]. Figure 16 summarizes the essential steps of $F_c$ computation.

The apodization equation is defined as follows:

$$r_k=w_k{y}_k;k=1..N_p$$

(54)

where $w_k$ denotes the Hamming window coefficients, $y_k$ are the samples of the observed signal $y\left(n\right)$ and $r_k$ becomes the coefficients resulting from the apodization.

The apodization step requires $P{N}_p$ multiplication operations. The second step in CFD detection process is the first FFT applied on apodized blocks. The complexity of FFT computation is known as $N_p{log}_2{N}_p$. Consequently, the P apodized blocks require $P{N}_{p}lo{g}_2{N}_p$ multiplication operations. The complex demodulates computation requires $P{N}_p^2$ multiplication operations and the step of complex demodulates multiplication requires $P^2{N}_p^2$ multiplication operations.

The second FFT with P data points to obtain the estimation of cyclic spectrum requires $N_p^{2}P{log}_{2}P$ multiplication operation and the calculation of crest factor $F_c$ alone requires $2N+1$ multiplication operations.

Finally, the algorithmic complexity of CFD, $C_{CFD}$, is:

$$C_{CFD}=34N^2+64NL+10N+2$$

(55)

6.3. Algorithmic Complexity of Recurrence Analysis Based Detector

To make decision, RAD uses statistic test defined in equation (18). RAD detection process can be summarized through the Figure 17.

The algorithmic complexity of RAD concern essentially the distance matrix computation and the calculation of average values of distance matrix coefficients. From the observed signal $y\left(n\right)$ containing N samples, we obtain $K=N-(m-1)\tau$ state vectors. Each state vector contains m coordinates. The distance $d_{ij}=∥{\mathbf{x}}_{-}{\mathbf{x}}_j∥$ form the state vector ${\mathbf{x}}_i$ to other state vectors ${\mathbf{x}}_j$ requires $mK$ additions operations. Because we exploit only the upper triangular part of distance matrix with of the main diagonal coefficients, we use $\frac{m}{2}K(K-2)$ addition operations. The computation of the average value of the first upper diagonal requires K elementary operations and the average value of other coefficients requires $\frac{{(K-2)}^2}{2}$ elementary operations. The statistic test of RAD requires one addition operation. Consequently, the computation complexity $C_{RAD}$ of RAD becomes:

$$\begin{array}{c}{C}_{RAD}=2N^2+4N(\tau -m\tau -1)+\hfill \\ \hfill \left[2{\left(m\tau \right)}^2+2{\tau }^2-4m{\tau }^2-4\tau +2m\tau +6\right]\end{array}$$

(56)

From the equations (52), (52) and (56), we generated in Figure 18, the curves of evolution of the algorithmic complexity of the detectors ED, CFD and RAD, based on the number of samples contained in the observed signal $y\left(n\right)$. These results in Figure 18 show that ED the simplest algorithm, whereas CFD is the most complicated one. We also notice that RAD algorithm is less complicated than CFD algorithm.

7. Conclusions

This paper deals with the problematic of non-cooperative spectrum sensing in very low SNR conditions. Many algorithms have been recently developed to overcome the scarcity of radio spectrum. However, most of them suffer from the noise uncertainty and do not work correctly in a very low SNR conditions. In this paper, we use the promising approach of the Recurrence Quantification Analysis (RQA) to propose a robust detection model, named the Recurrence Analysis based Detection (RAD). RAD benefits from the exploitation of the similitude among the different state vectors. Indeed, our analyzes reveal that for a White Gaussian Noise, the coefficients contained on the first upper diagonal are representative of other coefficients of distance matrix; which is not the case for a communication signal. Thus, by applying a conformity test between the coefficients of the first upper diagonal and other coefficients of the distance matrix, the presence or absence of a communication signal can be revealed. RAD presents five major advantages: It is more robust compared to Energy based Detector (ED) and Cyclostationary Feature Detector (CFD), that are widely used in the non-cooperative spectrum sensing context; it does not suffer from the noise variance estimation because the estimation of the noise variance is not required during the spectrum sensing process; it is able to detect the communication signal in a very low SNR condition; contrary to ED, RAD is able to distinguish a noisy communication signal and a high energy noise; RAD does not need a high computational cost like CFD. Our present work also presents two major contributions. First, we have determined, for Digital communication signals, the optimal values of the time delay $\tau$ and embedding dimension m needed for the phase space reconstruction. Second, we have established the analytical expression of detection threshold $\lambda$, the probability of detection $P_d$ and the probability of false alarm $P_{fa}$ of detection model based on RQA. From ROC curves, we can notice without ambiguity that the RAD is more robust than ED and CFD algorithms. Our current simulations show that RAD is able to detect the communication signal for SNR ≥ -12 dB. In addition to facilitating the blind detection of communication signals, RQA could be used to accurately estimate the characteristic frequencies of the signal of interest. In our future work, we will optimize the performance of the RAD detector based on this reliable estimation of the characteristic frequencies of the signal of interest.