# A Cooperative Spectrum Sensing Method Based on Empirical Mode Decomposition and Information Geometry in Complex Electromagnetic Environment

Authors: Yonghua Wang; Shunchao Zhang; Yongwei Zhang; Pin Wan; Jiangfan Li; Nan Li
Journal: Complexity
 (2019)

---

## Abstract

In a complex electromagnetic environment, there are cases where the noise is uncertain and difficult to estimate, which poses a great challenge to spectrum sensing systems. This paper proposes a cooperative spectrum sensing method based on empirical mode decomposition and information geometry. The method mainly includes two modules, a signal feature extraction module and a spectrum sensing module based on K-medoids. In the signal feature extraction module, firstly, the empirical modal decomposition algorithm is used to denoise the signals collected by the secondary users, so as to reduce the influence of the noise on the subsequent spectrum sensing process. Further, the spectrum sensing problem is considered as a signal detection problem. To analyze the problem more intuitively and simply, the signal after empirical mode decomposition is mapped into the statistical manifold by using the information geometry theory, so that the signal detection problem is transformed into geometric problems. Then, the corresponding geometric tools are used to extract signal features as statistical features. In the spectrum sensing module, the K-medoids clustering algorithm is used for training. A classifier can be obtained after a successful training, thereby avoiding the complex threshold derivation in traditional spectrum sensing methods. In the experimental part, we verified the proposed method and analyzed the experimental results, which show that the proposed method can improve the spectrum sensing performance.


---

## Body

## 1. Introduction

Spectrum sensing is a key step in the cognitive radio (CR) technology [1, 2]. The main task is to quickly and accurately detect whether the primary user (PU) is using the spectrum. Traditional spectrum sensing methods include energy detection [3], matched filter detection [4, 5], and cyclostationary feature detection [6]. These methods all have shortcomings. For example, in the case of uncertain noise, the detection performance of energy detection degrades sharply. Matching filter detection needs to know the PU signal and noise information in advance.In complex electromagnetic environments, noise uncertainty greatly impairs the performance of spectrum sensing [7, 8]. Therefore, it is necessary to reduce the noise and redundant information in the signals collected by the secondary users (SUs). The main challenge is that most of the signals collected by the SUs are nonstationary and nonlinear. To cope with that, traditional noise reduction methods are proposed, including short time Fourier transform, Wigner Ville distribution, and wavelet transform [9, 10]. However, these methods are all based on Fourier transforms, limited by the principle of uncertainty. Boudraa et al. proposed an empirical mode decomposition (EMD) algorithm, which can decompose the noise-infected signal into multiple intrinsic mode function (IMF) components [11]. An EMD denoising method based on continuous mean square error criterion is proposed [12]. Through analyzing the energy abrupt point between IMF components, the low-order IMF component before the mutation point is filtered out to achieve noise reduction.With the rapid development of information geometry, the concept of statistical manifolds can be applied to the signal detection problems. These problems can be transformed into geometric problems on manifolds and then analyzed using geometric tools, so as to indirectly solve the signal detection problems. Liu et al. used information geometry for radar signal detection and proposed the matrix constant false alarm probability and a geodesic distance detector [13]. Lu et al. applied the information geometry theory to spectrum sensing and obtained the close-form expression of the decision threshold using the matching methods, but the algorithm has higher complexity [14]. Chen et al. added the manifold geometry measurement and obtained the decision threshold through simulation [15]. However, all these methods need to derive and calculate the decision threshold, which is not only complicated but also inaccurate.The application of machine learning in various fields has also attracted great interest of many researchers [16, 17]. Thilina et al. used K-means, the Gaussian mixture model (GMM) in unsupervised learning and neural network (NN), support vector machine (SVM) in supervised learning to study spectrum sensing [18]. A spectrum sensing method based on signal energy is proposed, which uses K-means to classify these energy features. At the same time, the spectrum sensing performance of the algorithm under different channel models is further analyzed [19]. A spectrum sensing method based on dominant features and maximum and minimum eigenvalues is proposed, and K-means and GMM are selected to form the unsupervised learning framework [20]. Clustering algorithm and eigenvalue-based algorithms are also proposed. In the signal feature extraction phase, a feature extraction method is proposed to increase the number of cooperative SUs logically. Then, K-means and K-medoids are selected as the learning framework [21].Based on the previous studies, this paper proposes a cooperative spectrum sensing method based on empirical mode decomposition and information geometry (EMDIGK), to further improve the spectrum sensing performance in complex electromagnetic environments. In the feature extraction phase, the EMD algorithm is first used to denoise the signals collected by the SUs, as to obtain the signal characteristics more accurately. We introduce the order decomposition and recombination (ODAR), interval decomposition and recombination (IDAR) to logically increase the number of SUs. Then, the covariance matrix of the split and recombined signal matrix are calculated separately, and then the covariance matrix is mapped onto the manifolds using information geometry theory, and the geodesic distance on the manifold is calculated and used as the signal feature. Then, the K-medoids clustering algorithm is used in the training and spectrum sensing phase. The EMDIGK method proposed in this paper does not need to obtain the signal information of the communication PU in advance, so it is easy to be obtained for the labeled training data. In the experimental part, we verified the efficacy of EMDIGK through simulation. The experimental results show that EMDIGK successfully improves the spectrum sensing performance in complex electromagnetic environment.For ease of reference, the symbols and notations used in this paper are summarized in Table1.Table 1
Symbols and notations.
Symbols Notations wn Gaussian noise in the environment sn The signal transmitted by the PU xn The signal received by the SU N The number of sampling points M The number of SU H 1, H0 PU signal exists, PU signal does not exist P f, Pd False alarm probability, detection probability, xi The signal collected by thei-th SU X aM×N signal matrix IMFj(n) Thej-th intrinsic mode function ximax(n) The maximum envelope ofxi ximin(n) The minimum envelope ofxi m(n) Average value ofximax(n) and ximin(n) h(n) Component ofxi r(n) Residual ytn Reconstructed signal t Critical point x~i x i after EMD processing X~ A newM×N matrix consisting of x~in q The split parameter s The length of the split signal vector after splitting Y O D A R, YIDAR A matrix after ODAR or IDAR R O, RI Covariance matrix ofYODAR or YIDAR x x is an n dimensional sample θ θ is an m dimensional parameter vector Ζ A random variableΖ Θ Θ is the probability distribution space R Covariance matrix corresponding toX R w, Rs+Rw Covariance matrix corresponding toX under H0 and H1 θ A point on the manifold θt Arbitrary curve between the two points ηu Theu-th eigenvalues of the matrix R1-1/2R2R1-1/2 P The number of environmental noise signal matrices R p O, RpI Thep-th noise signal covariance matrix after ODAR and IDAR τ Iteration step size v Iteration step R - O, R-I Riemann Mean ofRO and RI d 1, d2 Distance between the sensing signal and reference pointR-O and R-I on the manifold Dj Thej-th Geodesic distance feature vector (GDFV) D- Training set Ck The set of training feature vectors belonging to classk ψk Center point ofCk D^ The feature extracted under the channel of the unknown PU state


## 2. The Basic System Model of Cooperative Spectrum Sensing

In the cognitive radio network (CRN), from the perspective of a single SU, the existence of the PU signal can be defined as a binary hypothesis as follows, whereH0 indicates that the PU signal does not exist and H1 indicates that the PU signal exists.(1)xn=wn,H0sn+wn,H1n=1,2,…,Nsn represents the signal transmitted by the PU and wn represents the Gaussian noise in the environment, and N denotes the number of sampling points. Based on (1), the detection probability (Pd) and the false alarm probability (Pf) of the system can be defined as (2) and (3), respectively. (2)Pf=PH1∣H0(3)Pd=PH1∣H1For a proposed algorithm, ifPf remains unchanged, a larger Pd indicates better performance. Similarly, a smaller Pf is desired when Pd is fixed.In CR system, spectrum sensing is usually performed in a complex environment, and then single SU sensing is often affected by multipath loss, shadows, and hidden terminals, resulting in performance degradation of the entire system [22]. In CRN, the fusion center (FC) collects the information collected by all SUs involved in cooperative sensing and then makes a final judgment based on this information [23–26]. The cooperative spectrum sensing system model is shown in Figure 1.Figure 1
The basic system model of cooperative spectrum sensing.In a CRN, we suppose that there areM SUs participating in cooperative spectrum sensing, then the signals collected by M SUs form a signal matrix X=[x1,x2,…,xM]T, where xi=[xi1,xi2,…,xi(N)] represents the signal acquired by the i-th SU. Therefore, X is a M×N matrix.(4)X=x1,x2,…,xMT=x11x12⋯x1Nx21x22⋯x2N⋮⋮⋱⋮xM1xM2⋯xMN


## 3. Feature Extraction Based on EMD and Information Geometry

### 3.1. Feature Extraction Model

In a complex environment, it is first necessary to estimate the noise after EMD denoising to establish a representative reference point. The process of feature extraction based on empirical modal decomposition and information geometry is shown in Figure2. To make the reference point more representative, we first collect enough noise signals, and then perform the EMD noise reduction, DAR, and covariance conversion (as shown in the dashed box in Figure 2). Then the Riemann mean solution method is used to solve the Riemann mean of these covariance matrices and is used as the reference point. Similarly, the signal matrix to be collected is also transformed by EMD, split recombination, and covariance matrix. Finally, the geodesic distances of these covariance matrices from the reference point are calculated and used as statistical features of the signal.Figure 2
The process of feature extraction based on EMD and information geometry.


### 3.2. Empirical Modal Decomposition

Assuming the signal acquired by thei-th SU is xi=[xi1,xi2,…,xi(N)], and the aim is to reduce the noise and redundant information in the collected signal, thereby improving the overall sensing system performance. Firstly, the EMD is used to reduce the noise of the signal collected by the SU. EMD can adaptively decompose any complex signal into a series of IMFs, shown as follows.(5)xin=∑j=1JIMFjn+rJnwhererJ represents the residual. The specific steps of the EMD decomposition algorithm are listed as follows.Step 1.
Find all local maxima and local minima of signalxi.Step 2.
The cubic spline interpolation method is used to fit all local maxima and local minima, respectively, to construct the maximum envelopeximax(n) and the minimum envelope ximin(n) and then calculate the average of ximax(n) and ximin(n). Therefore, the average value m1(n) can be obtained as(6)m1n=ximaxn+ximinn2Step 3.
By subtractingm1(n) from the original signal xi, the first component h1(n) can be obtained as(7)h1n=xin-m1n
Check whetherh1(n) satisfies the conditions of the IMF. If so, continue to Step 4; otherwise, redo Steps 1 and 2 for h1(n) to obtain the envelope value m11(n) and then obtain h11n as(8)h11n=h1n-m11n
Execute sequentially until thek-th step h1kn satisfies the conditions of the IMF, then(9)IMF1=h1knStep 4.
Subtractexi from the IMF1 component, resulting in the first residual r1(n) as(10)r1n=xin-IMF1n
Treatr1(n) as the original signal, and repeat Steps 1~4 above, to obtain r2(n), and so on, until rJ(n) becomes a monotonic function or constant.According to the characteristics of EMD, after the original signal is decomposed by EMD, several IMF components and one residual are obtained. The main component of the low-order IMF component is the high-frequency part of the signal, which contains sharp signals and noise. The high-order IMF component has less noise components, mainly the low-frequency part of the signal. According to this feature, several IMF components in the low frequency band can be selected to reconstruct the original signal to achieve the purpose of noise reduction. The reconstructed signal is shown in (11).(11)ytn=∑j=tJIMFjn+rJnIn the process of selecting the IMF component of the low frequency band, there must be a certain criticalIMFt component, so that the IMF component after the component is useful for the main part. Therefore, the ultimate goal of EMD decomposition is to find the IMFt component of the criticality. A continuous mean square error criterion is proposed in [10], shown as follows.(12)CMSEy~t,y~t+1=1N∑nNy~tn+y~t+1n2=1N∑n=1NIMFtn2t=argmin⁡CMSEy~t,y~t+11≤t≤J-1Therefore, the signal after the EMD processing can be obtained as(13)x~in=∑j=tJIMFjn+rJnAfter all the SUs collected signalxi after the EMD noise reduction, we can get a new matrix X~=[x~1,x~2,…,x~M], where x~i=x~i1,x~i2,…,x~iN represents the signal acquired by the i-th SU after EMD. X~ is an M×N matrix.(14)X~=x~11x~12⋯x~1Nx~21x~22⋯x~2N⋮⋮⋱⋮x~M1x~M2⋯x~MN


### 3.3. Decomposition and Recombination

After the EMD processing, the signal matrix after EMD denoising is processed using ODAR and IDAR, respectively [27], as to logically increase the number of cooperative SUs and further acquire the characteristic information of the signal. Let YODAR and YIDAR denote the signal matrix after ODAR and IDAR, respectively. Both matrices are qM×s. YODAR and YIDAR are defined as follows.(15)YODAR=x~11⋮x~1q⋮x~im⋮x~Mq=x~11x~12⋯x~1s⋮x~1q-1s+1x~1q-1s+1⋯x~1qs⋮x~im-1s+1x~im-1s+2⋯x~ims⋮x~Mq-1s+1x~Mq-1s+2⋯x~Mqs(16)YIDAR=x~11⋮x~1q⋮x~im⋮x~Mq=x~11x~1q+1⋯x~1s-1q+1⋮x~1qx~1q+q⋯x~1s-1q+q⋮x~imx~iq+m⋯x~is-1q+m⋮x~Mqx~Mq+q⋯x~Ms-1q+qwhereqq>0 is the split parameter and s is the length of the split signal vector after splitting. The covariance matrices RO and RI are then calculated using YODAR and YIDAR, shown in (17)RO=1sYODARYODART(18)RI=1sYIDARYIDART


### 3.4. Information Geometry Theory

In information geometry theory, consider a set of probability density functionspx∣θ, where x is an n dimensional sample associated with a random variable Ζ, x∈Ζ⊆Cn. θ is an m dimensional parameter vector, θ∈Θ⊆Cm. Therefore, the probability distribution space can be described by parameter set Θ. The probability distribution function family S is as shown in (19)S=px∣θ∣θ∈Θ⊆CmR is the covariance matrix composed of the signals collected by the SUs, shown as(20)R=1NXXTwhere X denotes the signal matrix described above, N denotes the number of samples, and XT denotes the transposition of X. We can parameterize the probability distribution family by R as(21)F=px∣R∣R∈CM×MwhereCM×M is the open set of the M×M dimensional vector space. According to the information geometry theory, F forms a microscopic manifold structure, which is called a statistical manifold. R is the coordinate of the manifold. Since the parameters of the manifold are covariance matrices, they can also be called matrix manifolds. Under the two hypotheses H0 and H1, we can obtain two types of covariance matrices, Rw and Rs+Rw. At the same time, Rw and Rs+Rw correspond to two types of coordinates on the manifold. As the SNR increases, the difference between the two hypothetical sensing signals becomes larger, corresponding to an increase in the distance between the two points on the matrix manifold [14]. In the following, based on this, we will extract the signal characteristics.


#### 3.4.1. Geodesic Distance

In information geometry, the distance between two probability distributions on a statistical manifold can be measured in a variety of ways. On the manifold, the geodesic distance is a widely accepted method for measuring the distance between two probability distributions. Due to the nature of the manifold curvature, the distance between the two points on the manifold depends on the choice of the curve between the two points. A curve that minimizes the distance between two points is defined as a geodesic, and the corresponding distance is called a geodesic distance.Assuming any two pointsθ1 and θ2 on the manifold, the arbitrary curve between the two points is θt(t1≤t≤t2), θt1=θ1, θt2=θ2. Then the curve distance between points θ1 and θ2 is expressed as [28]:(22)Dθ1,θ2≜∫t1t2dθdtTGθdθdtdtwhereG(θ) is the Fisher information matrix in the information geometry for statistical metrics on the manifold.(23)Gθ=E∂ln⁡px∣θ∂θi·∂ln⁡px∣θ∂θjWhen the minimum value is obtained using (22), i.e., the minimum distance between θ1 and θ2, the geodesic distance is obtained, and the corresponding curve is the geodesic.In the cognitive radio spectrum sensing technology, the sensing signal matrix is a multivariate Gaussian distribution family with the same mean but different covariance matrices [29, 30]. Assume having the covariance matrices R1 and R2, the geodesic distance [31] between them is as shown in (24)DR1,R2=12trlog2R1-1/2R2R1-1/2=12∑u=1nlog2⁡ηuwhereηu is the u-th eigenvalues of the matrix R1-1/2R2R1-1/2.


#### 3.4.2. Riemann Mean

We perform ODAR and IDAR on theP environmental noise matrices to obtain RpO(p=1,2,…,P) and RpI(p=1,2,…,P) and then find their Riemann mean. The objective function Φ is shown in (25)ΦR-O=1P∑p=1PD2RpO,R-O(26)ΦR-I=1P∑p=1PD2RpI,R-IwhereD(·,·) is the geodesic distance between two points on the manifold. The matrices R-O and R-I (which minimize the value of the objective function Φ(·)) are Riemann mean [32], which can be expressed as (27)R-O=arg⁡minRO∈Θ⁡ΦR-O(28)R-I=arg⁡minRI∈Θ⁡ΦR-IAssume that whenP=2, there are two points R1 and R2 on the manifold, and its Riemann mean R- is located at the midpoint of the geodesic of the two points R1 and R2 on the manifold. R- is calculated using (29)R-=R11/2R1-1/2R2R1-1/21/2R11/2When there areP points (P>2), the Riemann mean is difficult to calculate. In this case, [31, 33] have proposed the iterative calculation of the Riemann mean R- using the gradient descent algorithm, shown in (30)R-v+1=R-v1/2e-τ/P∑p=1Plog⁡R-v-1/2RpR-v-1/2R-v1/2,0≤τ≤1whereτ is the iteration step size and v is the iteration step. According to the above description, we can get the Riemann mean R-O of RpO(p=1,2,…,P) and R-I of RpI(p=1,2,…,P).The Riemann meansR-O and R-I are selected as reference points, and the signal matrix to be collected is also transformed by EMD, DAR, and covariance matrix to obtain RO and RI. Finally, the geodesic distances from RO to R-O and RI to R-I are calculated separately, and the distance is taken as the statistical characteristic of the signal.According to (24), we can calculate the corresponding geodesic distance as(31)d1=DRO,R-O=12trlog2RO-1/2R-ORO-1/2=12∑u=1Mqlog2⁡ηu(32)d2=DRI,R-I=12trlog2RI-1/2R-IRI-1/2=12∑u=1Mqlog2⁡ηuwhered1 and d2 represent the distance between the sensing signal and reference point R-O and R-I on the manifold, respectively. If the PU signal does not exist, then the distance from the reference points is small. On the contrary, if the PU signal exists, then the distance is large. According to such feature, we can intuitively analyze whether the PU is using the licensed spectrum. In order to facilitate the training of the clustering algorithm framework and test the spectrum sensing performance, a feature vector capable of reflecting the signal characteristics is needed. Therefore, the two geodesic d1 and d2 constitute a geodesic distance feature vector (GDFV) D=d1,d2.


## 3.1. Feature Extraction Model

In a complex environment, it is first necessary to estimate the noise after EMD denoising to establish a representative reference point. The process of feature extraction based on empirical modal decomposition and information geometry is shown in Figure2. To make the reference point more representative, we first collect enough noise signals, and then perform the EMD noise reduction, DAR, and covariance conversion (as shown in the dashed box in Figure 2). Then the Riemann mean solution method is used to solve the Riemann mean of these covariance matrices and is used as the reference point. Similarly, the signal matrix to be collected is also transformed by EMD, split recombination, and covariance matrix. Finally, the geodesic distances of these covariance matrices from the reference point are calculated and used as statistical features of the signal.Figure 2
The process of feature extraction based on EMD and information geometry.


## 3.2. Empirical Modal Decomposition

Assuming the signal acquired by thei-th SU is xi=[xi1,xi2,…,xi(N)], and the aim is to reduce the noise and redundant information in the collected signal, thereby improving the overall sensing system performance. Firstly, the EMD is used to reduce the noise of the signal collected by the SU. EMD can adaptively decompose any complex signal into a series of IMFs, shown as follows.(5)xin=∑j=1JIMFjn+rJnwhererJ represents the residual. The specific steps of the EMD decomposition algorithm are listed as follows.Step 1.
Find all local maxima and local minima of signalxi.Step 2.
The cubic spline interpolation method is used to fit all local maxima and local minima, respectively, to construct the maximum envelopeximax(n) and the minimum envelope ximin(n) and then calculate the average of ximax(n) and ximin(n). Therefore, the average value m1(n) can be obtained as(6)m1n=ximaxn+ximinn2Step 3.
By subtractingm1(n) from the original signal xi, the first component h1(n) can be obtained as(7)h1n=xin-m1n
Check whetherh1(n) satisfies the conditions of the IMF. If so, continue to Step 4; otherwise, redo Steps 1 and 2 for h1(n) to obtain the envelope value m11(n) and then obtain h11n as(8)h11n=h1n-m11n
Execute sequentially until thek-th step h1kn satisfies the conditions of the IMF, then(9)IMF1=h1knStep 4.
Subtractexi from the IMF1 component, resulting in the first residual r1(n) as(10)r1n=xin-IMF1n
Treatr1(n) as the original signal, and repeat Steps 1~4 above, to obtain r2(n), and so on, until rJ(n) becomes a monotonic function or constant.According to the characteristics of EMD, after the original signal is decomposed by EMD, several IMF components and one residual are obtained. The main component of the low-order IMF component is the high-frequency part of the signal, which contains sharp signals and noise. The high-order IMF component has less noise components, mainly the low-frequency part of the signal. According to this feature, several IMF components in the low frequency band can be selected to reconstruct the original signal to achieve the purpose of noise reduction. The reconstructed signal is shown in (11).(11)ytn=∑j=tJIMFjn+rJnIn the process of selecting the IMF component of the low frequency band, there must be a certain criticalIMFt component, so that the IMF component after the component is useful for the main part. Therefore, the ultimate goal of EMD decomposition is to find the IMFt component of the criticality. A continuous mean square error criterion is proposed in [10], shown as follows.(12)CMSEy~t,y~t+1=1N∑nNy~tn+y~t+1n2=1N∑n=1NIMFtn2t=argmin⁡CMSEy~t,y~t+11≤t≤J-1Therefore, the signal after the EMD processing can be obtained as(13)x~in=∑j=tJIMFjn+rJnAfter all the SUs collected signalxi after the EMD noise reduction, we can get a new matrix X~=[x~1,x~2,…,x~M], where x~i=x~i1,x~i2,…,x~iN represents the signal acquired by the i-th SU after EMD. X~ is an M×N matrix.(14)X~=x~11x~12⋯x~1Nx~21x~22⋯x~2N⋮⋮⋱⋮x~M1x~M2⋯x~MN


## 3.3. Decomposition and Recombination

After the EMD processing, the signal matrix after EMD denoising is processed using ODAR and IDAR, respectively [27], as to logically increase the number of cooperative SUs and further acquire the characteristic information of the signal. Let YODAR and YIDAR denote the signal matrix after ODAR and IDAR, respectively. Both matrices are qM×s. YODAR and YIDAR are defined as follows.(15)YODAR=x~11⋮x~1q⋮x~im⋮x~Mq=x~11x~12⋯x~1s⋮x~1q-1s+1x~1q-1s+1⋯x~1qs⋮x~im-1s+1x~im-1s+2⋯x~ims⋮x~Mq-1s+1x~Mq-1s+2⋯x~Mqs(16)YIDAR=x~11⋮x~1q⋮x~im⋮x~Mq=x~11x~1q+1⋯x~1s-1q+1⋮x~1qx~1q+q⋯x~1s-1q+q⋮x~imx~iq+m⋯x~is-1q+m⋮x~Mqx~Mq+q⋯x~Ms-1q+qwhereqq>0 is the split parameter and s is the length of the split signal vector after splitting. The covariance matrices RO and RI are then calculated using YODAR and YIDAR, shown in (17)RO=1sYODARYODART(18)RI=1sYIDARYIDART


## 3.4. Information Geometry Theory

In information geometry theory, consider a set of probability density functionspx∣θ, where x is an n dimensional sample associated with a random variable Ζ, x∈Ζ⊆Cn. θ is an m dimensional parameter vector, θ∈Θ⊆Cm. Therefore, the probability distribution space can be described by parameter set Θ. The probability distribution function family S is as shown in (19)S=px∣θ∣θ∈Θ⊆CmR is the covariance matrix composed of the signals collected by the SUs, shown as(20)R=1NXXTwhere X denotes the signal matrix described above, N denotes the number of samples, and XT denotes the transposition of X. We can parameterize the probability distribution family by R as(21)F=px∣R∣R∈CM×MwhereCM×M is the open set of the M×M dimensional vector space. According to the information geometry theory, F forms a microscopic manifold structure, which is called a statistical manifold. R is the coordinate of the manifold. Since the parameters of the manifold are covariance matrices, they can also be called matrix manifolds. Under the two hypotheses H0 and H1, we can obtain two types of covariance matrices, Rw and Rs+Rw. At the same time, Rw and Rs+Rw correspond to two types of coordinates on the manifold. As the SNR increases, the difference between the two hypothetical sensing signals becomes larger, corresponding to an increase in the distance between the two points on the matrix manifold [14]. In the following, based on this, we will extract the signal characteristics.


### 3.4.1. Geodesic Distance

In information geometry, the distance between two probability distributions on a statistical manifold can be measured in a variety of ways. On the manifold, the geodesic distance is a widely accepted method for measuring the distance between two probability distributions. Due to the nature of the manifold curvature, the distance between the two points on the manifold depends on the choice of the curve between the two points. A curve that minimizes the distance between two points is defined as a geodesic, and the corresponding distance is called a geodesic distance.Assuming any two pointsθ1 and θ2 on the manifold, the arbitrary curve between the two points is θt(t1≤t≤t2), θt1=θ1, θt2=θ2. Then the curve distance between points θ1 and θ2 is expressed as [28]:(22)Dθ1,θ2≜∫t1t2dθdtTGθdθdtdtwhereG(θ) is the Fisher information matrix in the information geometry for statistical metrics on the manifold.(23)Gθ=E∂ln⁡px∣θ∂θi·∂ln⁡px∣θ∂θjWhen the minimum value is obtained using (22), i.e., the minimum distance between θ1 and θ2, the geodesic distance is obtained, and the corresponding curve is the geodesic.In the cognitive radio spectrum sensing technology, the sensing signal matrix is a multivariate Gaussian distribution family with the same mean but different covariance matrices [29, 30]. Assume having the covariance matrices R1 and R2, the geodesic distance [31] between them is as shown in (24)DR1,R2=12trlog2R1-1/2R2R1-1/2=12∑u=1nlog2⁡ηuwhereηu is the u-th eigenvalues of the matrix R1-1/2R2R1-1/2.


### 3.4.2. Riemann Mean

We perform ODAR and IDAR on theP environmental noise matrices to obtain RpO(p=1,2,…,P) and RpI(p=1,2,…,P) and then find their Riemann mean. The objective function Φ is shown in (25)ΦR-O=1P∑p=1PD2RpO,R-O(26)ΦR-I=1P∑p=1PD2RpI,R-IwhereD(·,·) is the geodesic distance between two points on the manifold. The matrices R-O and R-I (which minimize the value of the objective function Φ(·)) are Riemann mean [32], which can be expressed as (27)R-O=arg⁡minRO∈Θ⁡ΦR-O(28)R-I=arg⁡minRI∈Θ⁡ΦR-IAssume that whenP=2, there are two points R1 and R2 on the manifold, and its Riemann mean R- is located at the midpoint of the geodesic of the two points R1 and R2 on the manifold. R- is calculated using (29)R-=R11/2R1-1/2R2R1-1/21/2R11/2When there areP points (P>2), the Riemann mean is difficult to calculate. In this case, [31, 33] have proposed the iterative calculation of the Riemann mean R- using the gradient descent algorithm, shown in (30)R-v+1=R-v1/2e-τ/P∑p=1Plog⁡R-v-1/2RpR-v-1/2R-v1/2,0≤τ≤1whereτ is the iteration step size and v is the iteration step. According to the above description, we can get the Riemann mean R-O of RpO(p=1,2,…,P) and R-I of RpI(p=1,2,…,P).The Riemann meansR-O and R-I are selected as reference points, and the signal matrix to be collected is also transformed by EMD, DAR, and covariance matrix to obtain RO and RI. Finally, the geodesic distances from RO to R-O and RI to R-I are calculated separately, and the distance is taken as the statistical characteristic of the signal.According to (24), we can calculate the corresponding geodesic distance as(31)d1=DRO,R-O=12trlog2RO-1/2R-ORO-1/2=12∑u=1Mqlog2⁡ηu(32)d2=DRI,R-I=12trlog2RI-1/2R-IRI-1/2=12∑u=1Mqlog2⁡ηuwhered1 and d2 represent the distance between the sensing signal and reference point R-O and R-I on the manifold, respectively. If the PU signal does not exist, then the distance from the reference points is small. On the contrary, if the PU signal exists, then the distance is large. According to such feature, we can intuitively analyze whether the PU is using the licensed spectrum. In order to facilitate the training of the clustering algorithm framework and test the spectrum sensing performance, a feature vector capable of reflecting the signal characteristics is needed. Therefore, the two geodesic d1 and d2 constitute a geodesic distance feature vector (GDFV) D=d1,d2.


## 3.4.1. Geodesic Distance

In information geometry, the distance between two probability distributions on a statistical manifold can be measured in a variety of ways. On the manifold, the geodesic distance is a widely accepted method for measuring the distance between two probability distributions. Due to the nature of the manifold curvature, the distance between the two points on the manifold depends on the choice of the curve between the two points. A curve that minimizes the distance between two points is defined as a geodesic, and the corresponding distance is called a geodesic distance.Assuming any two pointsθ1 and θ2 on the manifold, the arbitrary curve between the two points is θt(t1≤t≤t2), θt1=θ1, θt2=θ2. Then the curve distance between points θ1 and θ2 is expressed as [28]:(22)Dθ1,θ2≜∫t1t2dθdtTGθdθdtdtwhereG(θ) is the Fisher information matrix in the information geometry for statistical metrics on the manifold.(23)Gθ=E∂ln⁡px∣θ∂θi·∂ln⁡px∣θ∂θjWhen the minimum value is obtained using (22), i.e., the minimum distance between θ1 and θ2, the geodesic distance is obtained, and the corresponding curve is the geodesic.In the cognitive radio spectrum sensing technology, the sensing signal matrix is a multivariate Gaussian distribution family with the same mean but different covariance matrices [29, 30]. Assume having the covariance matrices R1 and R2, the geodesic distance [31] between them is as shown in (24)DR1,R2=12trlog2R1-1/2R2R1-1/2=12∑u=1nlog2⁡ηuwhereηu is the u-th eigenvalues of the matrix R1-1/2R2R1-1/2.


## 3.4.2. Riemann Mean

We perform ODAR and IDAR on theP environmental noise matrices to obtain RpO(p=1,2,…,P) and RpI(p=1,2,…,P) and then find their Riemann mean. The objective function Φ is shown in (25)ΦR-O=1P∑p=1PD2RpO,R-O(26)ΦR-I=1P∑p=1PD2RpI,R-IwhereD(·,·) is the geodesic distance between two points on the manifold. The matrices R-O and R-I (which minimize the value of the objective function Φ(·)) are Riemann mean [32], which can be expressed as (27)R-O=arg⁡minRO∈Θ⁡ΦR-O(28)R-I=arg⁡minRI∈Θ⁡ΦR-IAssume that whenP=2, there are two points R1 and R2 on the manifold, and its Riemann mean R- is located at the midpoint of the geodesic of the two points R1 and R2 on the manifold. R- is calculated using (29)R-=R11/2R1-1/2R2R1-1/21/2R11/2When there areP points (P>2), the Riemann mean is difficult to calculate. In this case, [31, 33] have proposed the iterative calculation of the Riemann mean R- using the gradient descent algorithm, shown in (30)R-v+1=R-v1/2e-τ/P∑p=1Plog⁡R-v-1/2RpR-v-1/2R-v1/2,0≤τ≤1whereτ is the iteration step size and v is the iteration step. According to the above description, we can get the Riemann mean R-O of RpO(p=1,2,…,P) and R-I of RpI(p=1,2,…,P).The Riemann meansR-O and R-I are selected as reference points, and the signal matrix to be collected is also transformed by EMD, DAR, and covariance matrix to obtain RO and RI. Finally, the geodesic distances from RO to R-O and RI to R-I are calculated separately, and the distance is taken as the statistical characteristic of the signal.According to (24), we can calculate the corresponding geodesic distance as(31)d1=DRO,R-O=12trlog2RO-1/2R-ORO-1/2=12∑u=1Mqlog2⁡ηu(32)d2=DRI,R-I=12trlog2RI-1/2R-IRI-1/2=12∑u=1Mqlog2⁡ηuwhered1 and d2 represent the distance between the sensing signal and reference point R-O and R-I on the manifold, respectively. If the PU signal does not exist, then the distance from the reference points is small. On the contrary, if the PU signal exists, then the distance is large. According to such feature, we can intuitively analyze whether the PU is using the licensed spectrum. In order to facilitate the training of the clustering algorithm framework and test the spectrum sensing performance, a feature vector capable of reflecting the signal characteristics is needed. Therefore, the two geodesic d1 and d2 constitute a geodesic distance feature vector (GDFV) D=d1,d2.


## 4. Cooperative Spectrum Sensing Based on K-Medoids Clustering

After extracting the features, we use the K-medoids clustering algorithm to perform spectrum sensing. The traditional spectrum sensing method needs to derive the threshold value to determine whether the PU exists. Nevertheless, these methods often have problems, such as inaccurate thresholds and difficulties in calculation. EMDIGK uses the K-medoids clustering algorithm as the learning framework. This algorithm is similar to the K-means clustering algorithm. Unlike K-means, the K-medoids clustering algorithm does not calculate the mean of the samples as the centroid point in the class. Instead, the actual sample point in the class is taken as the center point, and the center point has the smallest distance to all sample points in the class. Compared with the centroid point in K-means, the center point in K-medoids has the advantage of being less affected by extreme values [34]. The K-medoids algorithm is less sensitive to noise points, so that the outliers will not cause excessive deviation of the results of the division. The whole process is divided into two phases, the training and the spectrum sensing phase. The cooperative spectrum sensing system model based on K-medoids is shown in Figure 3.Figure 3
System model based on the K-medoids clustering algorithm.Before training, we need to prepare a training setD-:(33)D-=D1,D2,…,DLwhereDl is the feature vector extracted by the method described in the third section, and L represents the number of feature vectors in the training set D-. Let Ck denote the set of the training feature vectors belonging to class k, where k=1,2,…,K.(34)Ck=Dl∣Dl∈Clusterk∀lClass Ck has a centroid ψk. The training process is as follows.Step 1.
Input the training sample setD- and the number of clusters K.Step 2.
Initialize the centroidsψ1,ψ2,…,ψK.Step 3.
Calculate the distance of each sample to each centroid and put it into the nearest class.Step 4.
Updateψk with ψk=arg⁡minDl∈Ck⁡∑Dj∈CkDl-ψk2.Step 5.
CalculateminC1,…,CKψ1…,ψK⁡∑k=1K∑Dl∈CkDl-ψk2; if the result of the calculation is not changing, the algorithm stops; otherwise, return to Step 3.Step 6.
Outputψ1,ψ2,…,ψK.After training, according to the centroid of each class, we can construct a classifier for spectrum sensing, as shown in (35).(35)fD^=D^-ψ1mink=2,3,…,K⁡D^-ψkIn (35), D^ represents the feature extracted under the channel of the unknown PU state. If fD^>ξ, it means the PU signal exists, and thus the channel is not available. Otherwise, fD^<ξ means the PU signal does not exist, and the channel can be used. Parameter ξ is used to control the false alarm probability and the missed detection probability in the spectrum sensing process.The overall flow of the EMDIGK algorithm proposed in this paper is as shown in Algorithm1.Algorithm 1:The overall process of EMDIGK.
Algorithm: EMDIGK Input:Sensing signal matrix X Step 1: X is denoised by EMD to obtain an matrix X~ Step 2:Perform ODAR and IDAR on X~ to obtain YODAR and YIDAR, then calculate the corresponding covariance matrix. Step 3:Extract signal characteristics using information geometry theory to obtain D=d1,d2. Step 4:Construct a training vector set D- and train the K-medoids clustering algorithm to get Eq. (35) Step 5:Import the training set D^ and get test result according to Eq. (35) Step 6:Calculate Pd and Pf by Eq. (3) and Eq. (4) Output: P d and Pf


## 5. Simulation Results and Performance Analysis

In this section, we used the multicomponent signalSt=cos⁡t+cos⁡4t+0.2t2 as the experimental simulation signal. To ensure the accuracy of the experiment, we obtained 2,000 signal features, of which 1,000 were used as the training set and the other 1,000 were used as the test set. The test set is used to verify the spectrum sensing performance of the proposed method. We compare the detection performance under different sensing methods, and the experimental results suggest that the proposed EMDIGK-based method can obtain a better sensing performance.


### 5.1. The Clustering Effect of Clustering Algorithm

Figures4 and 5 show the effect of GFV before and after clustering by K-medoids at SNR=-14 dB. It can be seen that under such condition, after the K-medoids clustering algorithm, the signal characteristics in two different states can be well distinguished. The red ‘∗’s indicate that the PU exists, while the blue ‘x’s represent the PU does not. The black circle and triangle in the figure represent the centroid of the noise class and the centroid of the containing PU signal class, respectively.Figure 4
Classification of the K-medoids clustering algorithm with GFV.Figure 5
Classification of the K-medoids clustering algorithm with GFV.


### 5.2. The ROC Curve under Different SNR Conditions

Figures6 and 7 illustrate the ROC curves of different methods at two different SNR, respectively. IQMSE, IQDMM, IQMME, and IQRMET are the methods proposed in literature [14], and ED is a method using energy characteristics as a statistic.Figure 6
The ROC curve with GFV underSNR=-14 dB.Figure 7
The ROC curve with GFV underSNR=-16 dB.Tables2 and 3 show the detection probabilities of the methods under different SNR when the false alarm probability is a constant. We observe that the EMDIGK method has higher detection probability than other methods. It can be calculated that the detection probability of EMDIGK is 36.6%, 38.5%, 56.4%, 61.6%, and 136.0% higher than other methods under the condition of SNR=-14 dB and Pf=0.1. That is because the EMDIGK method reduces the noise in the communication signal, thereby reducing the impact of noise on the spectrum sensing system. At the same time, the feature of the information geometry method can be used to analyze the feature values in different states more intuitively, so that a better detection performance can be obtained.Table 2
ThePd of different methods in SNR=-14 dB.
Algorithm P f = 0.1 P f = 0.2 EMDIGK 0.97 0.98 IQMSE 0.71 0.98 IQDMM 0.70 0.77 IQRMET 0.62 0.74 IQMME 0.60 0.73 ED 0.41 0.51Table 3
ThePd of different methods in SNR=-16 dB.
Algorithm P f = 0.1 P f = 0.2 EMDIGK 0.68 0.78 IQMSE 0.42 0.60 IQDMM 0.35 0.53 IQRMET 0.31 0.49 IQMME 0.30 0.49 ED 0.27 0.41


### 5.3. The ROC Curve with Different Number of Cooperative SUs

Figure8 shows the ROC curve of the EMDIGK method with different SU numbers. The number of sampling points N=1000. When SNR=-18 dB, it shows that the performance is better with the increase of the number of cooperative SUs. This is because cooperative spectrum sensing can reduce the interference of various factors such as multipath fading and shadow in the propagation environment. Therefore, with the increase of the cooperative SU, the anti-interference ability of the system is stronger and the performance is better.Figure 8
The ROC curve with different numbers of SUs underSNR=-18 dB.


### 5.4. The ROC Curve with Different Sampling Points

Figure9 shows the ROC curve with different sampling points for the EMDIGK method under SNR=-18 dB. The number of cooperative SUs M=2. As the number of sampling points increases, the information of the sensing signal is more complete, so the extracted features are more representative. As a result, from the figure we can see that when the number of sampling points increases, the detection performance of the EMDIGK method is also improved.Figure 9
The ROC curve with different numbers of sampling pointsSNR=-16 dB.


## 5.1. The Clustering Effect of Clustering Algorithm

Figures4 and 5 show the effect of GFV before and after clustering by K-medoids at SNR=-14 dB. It can be seen that under such condition, after the K-medoids clustering algorithm, the signal characteristics in two different states can be well distinguished. The red ‘∗’s indicate that the PU exists, while the blue ‘x’s represent the PU does not. The black circle and triangle in the figure represent the centroid of the noise class and the centroid of the containing PU signal class, respectively.Figure 4
Classification of the K-medoids clustering algorithm with GFV.Figure 5
Classification of the K-medoids clustering algorithm with GFV.


## 5.2. The ROC Curve under Different SNR Conditions

Figures6 and 7 illustrate the ROC curves of different methods at two different SNR, respectively. IQMSE, IQDMM, IQMME, and IQRMET are the methods proposed in literature [14], and ED is a method using energy characteristics as a statistic.Figure 6
The ROC curve with GFV underSNR=-14 dB.Figure 7
The ROC curve with GFV underSNR=-16 dB.Tables2 and 3 show the detection probabilities of the methods under different SNR when the false alarm probability is a constant. We observe that the EMDIGK method has higher detection probability than other methods. It can be calculated that the detection probability of EMDIGK is 36.6%, 38.5%, 56.4%, 61.6%, and 136.0% higher than other methods under the condition of SNR=-14 dB and Pf=0.1. That is because the EMDIGK method reduces the noise in the communication signal, thereby reducing the impact of noise on the spectrum sensing system. At the same time, the feature of the information geometry method can be used to analyze the feature values in different states more intuitively, so that a better detection performance can be obtained.Table 2
ThePd of different methods in SNR=-14 dB.
Algorithm P f = 0.1 P f = 0.2 EMDIGK 0.97 0.98 IQMSE 0.71 0.98 IQDMM 0.70 0.77 IQRMET 0.62 0.74 IQMME 0.60 0.73 ED 0.41 0.51Table 3
ThePd of different methods in SNR=-16 dB.
Algorithm P f = 0.1 P f = 0.2 EMDIGK 0.68 0.78 IQMSE 0.42 0.60 IQDMM 0.35 0.53 IQRMET 0.31 0.49 IQMME 0.30 0.49 ED 0.27 0.41


## 5.3. The ROC Curve with Different Number of Cooperative SUs

Figure8 shows the ROC curve of the EMDIGK method with different SU numbers. The number of sampling points N=1000. When SNR=-18 dB, it shows that the performance is better with the increase of the number of cooperative SUs. This is because cooperative spectrum sensing can reduce the interference of various factors such as multipath fading and shadow in the propagation environment. Therefore, with the increase of the cooperative SU, the anti-interference ability of the system is stronger and the performance is better.Figure 8
The ROC curve with different numbers of SUs underSNR=-18 dB.


## 5.4. The ROC Curve with Different Sampling Points

Figure9 shows the ROC curve with different sampling points for the EMDIGK method under SNR=-18 dB. The number of cooperative SUs M=2. As the number of sampling points increases, the information of the sensing signal is more complete, so the extracted features are more representative. As a result, from the figure we can see that when the number of sampling points increases, the detection performance of the EMDIGK method is also improved.Figure 9
The ROC curve with different numbers of sampling pointsSNR=-16 dB.


## 6. Conclusions

This paper proposes an EMDIGK method to effectively improve the spectrum sensing performance in a complex electromagnetic environment. In terms of feature extraction, firstly, the EMD algorithm is used to denoise the signals collected by the SUs, thereby reducing the impact of noise on subsequent sensing process. Further, in order to analyze this kind of problem more intuitively and simply, we adopt the information geometry theory. The signal after EMD is mapped into the statistical manifold, thereby transforming the signal detection problem into a geometric problem. Then, the corresponding geometric tools are used to extract the signal features as statistical features. Finally, the K-medoids clustering algorithm is used as a learning framework to construct a spectrum-aware classifier to achieve spectrum sensing. In the experiments, the performance of the EMDIGK algorithm is verified and analyzed, which shows its superiority over other popular methods. In our future work, we will study the combination cyclostationary feature detection with clustering algorithms and information geometry. In future work, research could be extended in the area of how to reduce the complexity of the algorithm to achieve spectrum sensing more efficiently.



---

*Source: 5470974-2019-02-18.xml*