# A Novel Classification and Identification Scheme of Emitter Signals Based on Ward’s Clustering and Probabilistic Neural Networks with Correlation Analysis

Authors: Xiaofeng Liao; Bo Li; Bo Yang
Journal: Computational Intelligence and Neuroscience
 (2018)

---

## Abstract

The rapid development of modern communication technology makes the identification of emitter signals more complicated. Based on Ward’s clustering and probabilistic neural networks method with correlation analysis, an ensemble identification algorithm for mixed emitter signals is proposed in this paper. The algorithm mainly consists of two parts, one is the classification of signals and the other is the identification of signals. First, self-adaptive filtering and Fourier transform are used to obtain the frequency spectrum of the signals. Then, the Ward clustering method and some clustering validity indexes are used to determine the range of the optimal number of clusters. In order to narrow this scope and find the optimal number of classifications, a sufficient number of samples are selected in the vicinity of each class center to train probabilistic neural networks, which correspond to different number of classifications. Then, the classifier of the optimal probabilistic neural network is obtained by calculating the maximum value of classification validity index. Finally, the identification accuracy of the classifier is improved effectively by using the method of Bivariable correlation analysis. Simulation results also illustrate that the proposed algorithms can accurately identify the pulse emitter signals.


---

## Body

## 1. Introduction

Under the conditions of the rapid development of modern information technology, various types of communication equipment, such as radar and radio navigation equipment, radio and television equipment, and electronic computer and peripheral equipment, have been used by military and large technology companies. Its applications have been extended from the ground, air, and sea to outer space. In order to obtain important information promptly, accurately, and effectively, the signal characteristics derived from the same types of communication sources need to be extracted and analyzed. The realization of identification from general communication signals to individual signals has a significance to win the initiative. Therefore, many researchers have devoted to the study of the identification of emitter signals. Their purpose is to further improve the accuracy of signal classification and identification.The earliest research on the identification of emitter signals began in the 1970s, and it was one of the key technologies in Electronic Warfare systems. At present, the identification of emitter is mainly based on two modeling techniques: syntactic pattern-based methods and parametric pattern-based methods. In [1], Visnevski author had proposed a syntactic model to identify multifunction radar (MFR) signals. In the process of modeling, the MFRs were considered as stochastic discrete event systems that communicated information by use of radar word level modeling, radar phrase level modeling, and radar sentence level modeling. The radar word was a fixed arrangement of finite number of pulses, the radar phrase was a series of limited number of radar words, and the radar sentence was a combination of limited number of radar phrases. His simulation experiments had shown that the designed principle was effective for identifying MFRs. Based on the syntactic model of MFRs, Alex Wang and Vikram Krishnamurthy had used the stochastic context-free grammar to describe the behaviors of the MFR system, and some good results were obtained. Although the stochastic context-free grammar was a model for capturing the essential features of the MFR dynamics [2], it had some defects in estimating the parameters of stochastic context-free grammar. Therefore, the expectation maximization algorithm had been proposed by LP Dai et al. to estimate these parameters, which can be used to further estimate the characteristic parameters of MFR [3]. From this point of view, the ultimate goal of the modeling technique based on syntax was to find the feature parameters of emitters. It highlighted the importance of identification technology based on the parametric pattern. The feature parameter matching technique was a basic method of this pattern. The main characteristic of this method was to identify the emitter signals by matching the measured signal characteristic parameter vector with the corresponding characteristic parameters in the known database (the libraries of radar types). This method depended on the feature parameter database, and it can only be applied to emitter identification problem with invariant characteristic parameters. Moreover, the libraries of types had also an inherent uncertainty resulting inevitably from data collection methods. Therefore, Jan Matuszewski et al. had proposed the knowledge-based techniques to identify emitters [4, 5]. They think that the information of known radar platforms (knowledge), including position, intent, and recent operational history, plays an important role in the identification of emitters. Their approaches had been successfully used to identify some specific emitters. In order to fully utilize this knowledge to identify emitters, Janusz Dudczyk had proposed an idea of constructing Emitter DataBase based on the entity-relationship modeling. The entity-relationship diagram was introduced to realize this idea, which had pointed out a new direction for the construction of a complete and accurate Electronic Intelligence systems [6]. Meanwhile, artificial intelligence techniques and some optimized feature selection methods were used to improve the identification accuracy of emitter signals. In [7], the authors had proposed a vector neural network with a supervised learning algorithm which worked for signal classification and emitter identification. This network took carrier frequency, pulse width, and pulse repetition interval as inputs to complete the identification. In [8], the authors had proposed an identification method of radar signals based on immune radial-basis function neural network, which can improve the convergence speed and performance of the algorithm. In [9], a multichannel recognition system with an independent distance which was defined on impartiality condition had been proposed to identify the specific emitters. By modifying the distance in each special recognition channel, radio frequency, pulse width, and pulse repetition interval of radar signals can be extracted and classified into an appropriate class. Beyond that, some methods had been proposed to reduce the identification error rate. In [10], wavelet features were used as inputs of neural networks to identify emitters. In [11–13], the support vector machine was introduced to identify emitter signals. In [14, 15], the fuzzy c-means and probabilistic neural networks (PNN) were used to identify emitter signals. However, the determination of the optimal number of clusters was a major challenge for these methods, and they cannot effectively improve the identification accuracy. Therefore, Jawad et al. designed a clustering validity function in the hidden-layer output space of PNN to find the optimal number of clusters. Their methods were successfully applied to the classification of land use [16]. But, the method of determining the range of clustering number was subjective, which may reduce the identification accuracy of the algorithm. To overcome this deficiency and further improve accuracy, a classification and identification scheme of emitter signals based on the Ward clustering method (WCM) and PNN with correlation analysis was proposed in this paper. Its advantages are presented in three aspects:(1)
The self-adaptive filtering, Ward’s clustering, and clustering validity indexes are skillfully used to determine the scope of the optimal number of clusters.(2)
The classification validity indexD is flexibly used to find the optimal PNN classifier.(3)
The probabilistic neural networks with Bivariable correlation analysis approach are proposed to improve the identification accuracy of emitter signals.The rest of this paper is organized as follows. In Section2, the classification and identification schemes including adaptive filtering, frequency spectrum, evaluation indexes, WCM and PNN classifier are introduced. The flowchart and the pseudocodes of the classification algorithms are designed in Section 3. The flowchart and the pseudocodes of the identification algorithm are given in Section 4. The identification experiments are also carried out in this section. The comparison discussions of different schemes are proposed in Section 5. Some innovations and applicable conditions of the proposed method are summarized in Section 6.


## 2. Classification Model of Emitter Signals

### 2.1. Self-Adaptive Filtering

In the field of engineering technology, the signalxn received at time n usually contains two parts. One is the useful signal sn, which is what we need, and it enables us to understand the properties of the object to be studied. The other is the interference signal x1n, which is what we do not need, and it prevents us from understanding the properties of the object to be studied. The actual signal will be obtained once the two parts are combined together. That is, xn=sn+x1n.Weakening the interference signalx1n and maintaining or enhancing the useful signal sn are an important purpose of signal processing. The usual method is to use a frequency function Hf to multiply the frequency spectrum Xf of the signal xn. This process is called filtering. Its essence is to weaken the interference signal and highlight the useful signal. Currently, the most widely used filters include Kalman filtering, Wiener filtering, median filtering, sequential statistical filtering, wavelet transform, self-adaptive filtering, etc. In terms of adaptability and filtering performance, one of the best filtering methods is the self-adaptive filtering, which is developed on the basis of Kalman filtering, Wiener filtering, and linear filtering. The most important feature of the self-adaptive filtering is that it can track the time-varying characteristics of input signals and eliminate the unknown interference contained in the signals. Self-adaptive filtering based on the least mean square algorithm is proposed by Widrow and Hoff, and it has been widely used in many fields because of its simplicity, robustness, and easy implementation. The principle diagram of self-adaptive filtering technology is shown in Figure 1.Figure 1
Flowchart of the self-adaptive filtering technology.Figure1 presents the schematic diagram of noise elimination for self-adaptive filter. The actual signal xn contains interference signal x1n generated from signal channel 1. In order to eliminate it, the noise signal x0n which is independent of sn but related to x1n must be sampled from noise source through the signal channel 2. The main function of the self-adaptive filter is to process x0n so that the output yn approximates to x1n. Under the condition of convergence of filtering algorithm, the output en of the system approximates to sn when yn approaches x1n. The iterative formulas of self-adaptive filtering algorithm based on the least mean square are defined as follows [17]:(1)en=dn−ωTnXn,ωn+1=ωn+2μenXn,where dn is the desired signal; Xn=xn,xn−1,⋯,xn−M+1T is the input signal vector at time n; M is the length of filter; and μ is the fixed step size and satisfies 0<μ<1/MPin, where Pin is the input power of filter. ωn=w0n,w1n,⋯,wM−1nT is the weight vector of M order adaptive filter at time n.ωn=0 at initial time. yn=ωTnXn represents the actual output signal of filter. In noise elimination applications, xn is usually used as the desired signal dn and x0n is used as the input signal of the filter to eliminate x1n. After many iterations, the difference between dn and yn is the estimate of signal sn. This algorithm has the advantages of small amount of computation, easy implementation, and stable performance, but its convergence speed is relatively slow. Therefore, the authors proposed a variable step adaptive filtering algorithm based on bell-shaped function [18]. The variable step size is given as follows:(2)μn=μmax1−e−aenb,where μmax∈0,1/MPin is the maximum step size that can maintain the convergence of the adaptive filtering algorithm. Experiments show that this algorithm can effectively improve the convergence speed and reduce the steady-state error when a=0.08 and b=4. In subsection 3, the method is used to eliminate the interference signal of emitters.


### 2.2. The Spectrum of Signals

In order to classify and identify signals, it is necessary to analyze the frequency spectrum and energy spectrum of the signals. From the electrical knowledge,P=V2/R, where V represents voltage and R represents resistance. If the resistance R=1 and V is replaced by signal xt, the instantaneous energy is x2t. Thus, the total energy of the signals can be expressed as ∫−∞+∞x2tdt. According to Parseval’s theorem, it can be obtained by the following equation:(3)∫−∞+∞xt2dt=∫−∞+∞Xf2df,where(4)Xf=∫−∞+∞xte−2πiftdt,is the Fourier transform for signal xt and f represents the frequency of signal xt.Xf is called the amplitude spectrum and argXf is called the phase spectrum of xt. The frequency spectrum of xt consists of amplitude spectrum and phase spectrum. Xf2 is called energy spectrum density. Equation (3) indicates that the energy of xt is closely related to Xf2, and it can be obtained by calculating its integral on −∞,+∞. Therefore, we can obtain the frequency distribution and energy distribution of each signal by analyzing the spectrum of the signal. Then, the key amplitude and the frequency of energy distribution can be further obtained, which provide conditions for directly observing the similarities and differences of paired signals.


### 2.3. The Ward Clustering Method

As a hierarchical agglomerate cluster algorithm, the WCM has a wide range of applications [19–21]. First, it is started by accepting each node as a separate cluster. Then, the clusters with minimum distance between themselves are combined in pairs at each stage of the algorithm. This smallest distance is called the Ward distance and defined as follows:(5)drs=nr⋅nsnr+nsx¯r−x¯s2,where r and s represent the two distinct clusters, nr and ns represent the number of data points of two clusters, respectively, x¯r and x¯s represent the center of the corresponding cluster, and ||⋅|| is Euclidean norm. The centers and cardinal numbers of the new cluster are updated according to the following equations:(6)x¯r′=nr⋅x¯r+ns⋅x¯snr+ns,(7)nr′=nr+ns.Ward’s clustering algorithm has the following steps:Step 1. Each sample point is treated as a cluster. At this time, the sum of squares of deviations for each cluster is equal to 0.
Step 2. Two arbitrary clusters are merged, and the sum of squares of deviations is calculated from Equations (5)–(7). If we assume that there are N clusters in total, it must be calculated NN−1/2 times.
Step 3. The two clusters with the smallest squared sum of deviations are combined into one class. The method eventually aggregates all sample points into one class when the number of clusters is unknown.If the number of clusters is known, the WCM can be directly used to classify the signal data after removing the noise. Otherwise, it can be obtained by analyzing the dendrogram of clustering. This is a rather subjective approach, which is difficult to help us finding the true number of clusters for a given data set. In recent research, the clustering validity indexes, such as Calinski-Harabasz (CH) index, Gap index, Silhouette (Silh) index, and Davies–Bouldin (DB) index, have been demonstrated to be the best validation tools for determining the optimal number of clusters [22–27].


#### 2.3.1. Calinski-Harabasz Index

For a given setY=y1,y2,⋯,yN, assume that the dimension of each entity yi is v, i=1,2,⋯,N. K nonempty disjoint cluster sets S=S1,S2,⋯,SK around the centroid set C=c1,c2,⋯,cK can be obtained by minimizing the within-cluster distance WK:(8)WK=∑k=1K∑yi∈Skdyi,ck,where dyi,ck,k=1,2,⋯,K represents the squared Euclidean distance between the entity yi and the centroid ck, that is,(9)dyi,ck=∑j∈vyij−ckj2.Then, the CH index is defined as follows [24]:(10)CHK=T1−WK/K−1WK/N−K,where WK is defined as in (8), and T1 can be calculated by T1=∑i=1N∑j∈vyij−y¯j2.The CH index can reflect the compactness of the cluster by means of the overall within-cluster variance. The separation degree of the clusters can be reflected by the overall between-cluster variance. Therefore, a good clustering scheme corresponds to a higher value of CH index.


#### 2.3.2. Silhouette Index

For each entityyi, its silhouette value measures the similarity between the entity yi and the points in its own cluster, when compared to the points in other clusters. This similarity is reflected by measuring the distance between the entity yi and the points derived from different clusters. The silhouette value of the entity yi is defined as follows [24]:(11)syi=byi−ayimaxayi,byi,where ayi is the average distance from the entity yi∈Sk to all other points yj∈Sk, byi is the minimum distance from the entity yi to all other points yj, which satisfies yj∈Sl,yi∉Sl,l≠k. Therefore, −1≤syi≤1. If syi is close to zero, the entity yi could be assigned to another cluster. A negative value of syi indicates that the corresponding assignment seriously damages cluster cohesion, and the clustering result of yi is not advisable. yi is well matched to its own cluster when syi is close to 1. Finally, the validity of the whole clustering can be quantified by Silh index, and it is defined as follows:(12)SilhK=∑i=1NsyiN.The Silh index can be used with any distance metric, including the Manhattan distances and Euclidean distances.


#### 2.3.3. Davies–Bouldin Index

A good partition should have a larger intercluster separation degree and stronger within-cluster homogeneity and compactness. The DB index is proposed based on this idea [26]. More concretely, it is constructed by a ratio of within-cluster and between-cluster distances. The DB index is defined as follows:(13)DBK=∑k=1Kmaxk≠jd¯k+d¯jdkj/K,where(14)d¯k=1Sk∑i=1Skyi−ckq1/q,represents the average distance between each point yi in cluster k and the centroid of clusterk. Sk is the number of points in cluster k. If q=1, d¯k is the average Euclidean distance between the points in cluster k to the centroid of cluster k. If q=2, d¯k is the standard deviation of the distance of points in cluster k to the center of cluster k. When k in d¯k is replaced by j, d¯j can be obtained. In addition, dkj can be calculated according to the following equation:(15)dkj=∑h=1vc¯kh−c¯jhp1/p.It represents the distances between the centroids of thekth and the jth clusters. c¯kh is the hth component of the centroid of cluster k and dkj is the Minkowski metric of the centroids which characterizes clusters k and j. Specifically, p=1, dkj is the Manhattan distance between centroids p=2, dkj is the Euclidean distance between centroids.The DB index can reflect the degree of within-cluster dispersion and between-cluster separation. So, the true number of clusters may be determined according to the minimum value of the DB index.


#### 2.3.4. Gap Index

Robert Tibshirani et al. proposed the gap statistic method for estimating the number of clusters in a set of data [27]. A graph of the within-cluster dispersion versus the number of clusters k for a clustering procedure shows that the within-cluster dispersion decreases monotonically as k increases, but from some k onwards, the decrease becomes flatter obviously. Such position is called ‘elbow’, it often implies the appropriate number of clusters. The gap criterion gives an approach to estimate the number of clusters by locating this ‘elbow’. Therefore, under this criterion, the optimal number of clusters occurs at the largest gap value. The Gap index is defined as follows:(16)GapNK=ENlogUK−logUK,where N represents the number of points, K represents the number of clusters that are evaluated, UK defined in (17) represents the within-cluster dispersion degree.(17)UK=∑k=1K12NkHk,where Nk is the number of points in cluster k, Hk is the sum of the distances of any two points in the kth cluster. The expected value ENlogUK is determined by Monte Carlo sampling from a reference distribution. The Gap index can also be used for any distance metric.The WCM belongs to unsupervised categorization technique, which can help us find the centroid of each cluster. However, the classification accuracy of this method is limited, which makes the method not able to be used directly for signal recognition. Comparatively, PNN can effectively improve the accuracy of classification and identification [28, 29].


### 2.4. Probabilistic Neural Network Classifier

As a method of nonparametric Parzen windows estimation, PNN is first proposed by Specht. It is a nonlinear classification technique and essentially a parallel algorithm based on Bayesian minimum risk criterion [30]. Given a sample to be identified x, its posterior probability PSkx can be obtained by PNN classifier. However, if the probability densities of the classes to be separated are unknown, the training samples with known identity need to be used to estimate them. Finally, the trained PNN is used to determine the identity of x. A typical PNN classifier consists of an input layer, a pattern layer (hidden layer), a summation layer and a output layer. The flowchart of the PNN is shown in Figure 2.Figure 2
Flowchart of the PNN.The input layer neurons are used to receive values from training samples and send data to the neurons in the pattern layer, which is fully connected to the input layer. The number of neurons in the input layer is equal to the length of the input vector. The number of neurons in the pattern layer is the same as the number of training samples. Here, all neurons are collected into different groups, and theith neuron in group k corresponds to a Gaussian function fikx,σ,i=1,2,⋯,mk, where mk represents the number of neurons in group k, k=1,2,⋯,K. Gaussian function which is also called probability density function is defined as follows:(18)fikx,σ=12πv/2σvexp−∑j=1vxijk−xj22σ2,where v is the dimension of the input vector x=x1,x2,⋯,xv, xj is the jth component of the input vector x, xijk is the jth component of the ith neuron in class k. The so-called smoothing parameter σ∈0,1 determined experimentally by comparing their corresponding classification accuracy plays an important role in estimation error of the PNN classifier. The outputs of pattern layer are connected to the summation units depending on the class of patterns. There is one neuron for each group, and each neuron in summation layer sums the outputs derived from the pattern layer neurons as follows:(19)pkx=12πv/2σvmk∑i=1mkexp−∑j=1vxijk−xj22σ2.Finally, the output layer neuron output a number 1 and multiple numbers 0. The value of 1 corresponds to the classifier’s decision result for input vectors. More specifically, the input vectorx belongs to class k if pkx>pk′x for all k′=1,2,⋯,K and k≠k′.Hence, the main purpose of training PNN is to find the optimal estimate of probability density function according to the training samples and their labels, to ensure that the classifier works at the condition of minimum error rate and risk. When the samples to be identified are sent to the pattern layer, the output of each neuron is calculated according to the trained density function. Finally, the identified results are obtained through computations in the summation layer and output layer. Due to the following advantages, it is a wise choice to use PNN as a further classifier to classify signals [31]:(1)
It has a simple structure, and it is easy to train. In the PNN based on probability density function estimation, the weight of the neuron in pattern is directly taken from the input sample value.(2)
The training process of the network is simple, and there is no need to retrain for a long time when adding or reducing the number of groups.(3)
It is not easy to produce local optimal solution, and its precision is higher than that of other classification approaches. No matter how complex the classification problem is, as long as there are enough training samples, the optimal solution under the Bayes criterion can be obtained.


## 2.1. Self-Adaptive Filtering

In the field of engineering technology, the signalxn received at time n usually contains two parts. One is the useful signal sn, which is what we need, and it enables us to understand the properties of the object to be studied. The other is the interference signal x1n, which is what we do not need, and it prevents us from understanding the properties of the object to be studied. The actual signal will be obtained once the two parts are combined together. That is, xn=sn+x1n.Weakening the interference signalx1n and maintaining or enhancing the useful signal sn are an important purpose of signal processing. The usual method is to use a frequency function Hf to multiply the frequency spectrum Xf of the signal xn. This process is called filtering. Its essence is to weaken the interference signal and highlight the useful signal. Currently, the most widely used filters include Kalman filtering, Wiener filtering, median filtering, sequential statistical filtering, wavelet transform, self-adaptive filtering, etc. In terms of adaptability and filtering performance, one of the best filtering methods is the self-adaptive filtering, which is developed on the basis of Kalman filtering, Wiener filtering, and linear filtering. The most important feature of the self-adaptive filtering is that it can track the time-varying characteristics of input signals and eliminate the unknown interference contained in the signals. Self-adaptive filtering based on the least mean square algorithm is proposed by Widrow and Hoff, and it has been widely used in many fields because of its simplicity, robustness, and easy implementation. The principle diagram of self-adaptive filtering technology is shown in Figure 1.Figure 1
Flowchart of the self-adaptive filtering technology.Figure1 presents the schematic diagram of noise elimination for self-adaptive filter. The actual signal xn contains interference signal x1n generated from signal channel 1. In order to eliminate it, the noise signal x0n which is independent of sn but related to x1n must be sampled from noise source through the signal channel 2. The main function of the self-adaptive filter is to process x0n so that the output yn approximates to x1n. Under the condition of convergence of filtering algorithm, the output en of the system approximates to sn when yn approaches x1n. The iterative formulas of self-adaptive filtering algorithm based on the least mean square are defined as follows [17]:(1)en=dn−ωTnXn,ωn+1=ωn+2μenXn,where dn is the desired signal; Xn=xn,xn−1,⋯,xn−M+1T is the input signal vector at time n; M is the length of filter; and μ is the fixed step size and satisfies 0<μ<1/MPin, where Pin is the input power of filter. ωn=w0n,w1n,⋯,wM−1nT is the weight vector of M order adaptive filter at time n.ωn=0 at initial time. yn=ωTnXn represents the actual output signal of filter. In noise elimination applications, xn is usually used as the desired signal dn and x0n is used as the input signal of the filter to eliminate x1n. After many iterations, the difference between dn and yn is the estimate of signal sn. This algorithm has the advantages of small amount of computation, easy implementation, and stable performance, but its convergence speed is relatively slow. Therefore, the authors proposed a variable step adaptive filtering algorithm based on bell-shaped function [18]. The variable step size is given as follows:(2)μn=μmax1−e−aenb,where μmax∈0,1/MPin is the maximum step size that can maintain the convergence of the adaptive filtering algorithm. Experiments show that this algorithm can effectively improve the convergence speed and reduce the steady-state error when a=0.08 and b=4. In subsection 3, the method is used to eliminate the interference signal of emitters.


## 2.2. The Spectrum of Signals

In order to classify and identify signals, it is necessary to analyze the frequency spectrum and energy spectrum of the signals. From the electrical knowledge,P=V2/R, where V represents voltage and R represents resistance. If the resistance R=1 and V is replaced by signal xt, the instantaneous energy is x2t. Thus, the total energy of the signals can be expressed as ∫−∞+∞x2tdt. According to Parseval’s theorem, it can be obtained by the following equation:(3)∫−∞+∞xt2dt=∫−∞+∞Xf2df,where(4)Xf=∫−∞+∞xte−2πiftdt,is the Fourier transform for signal xt and f represents the frequency of signal xt.Xf is called the amplitude spectrum and argXf is called the phase spectrum of xt. The frequency spectrum of xt consists of amplitude spectrum and phase spectrum. Xf2 is called energy spectrum density. Equation (3) indicates that the energy of xt is closely related to Xf2, and it can be obtained by calculating its integral on −∞,+∞. Therefore, we can obtain the frequency distribution and energy distribution of each signal by analyzing the spectrum of the signal. Then, the key amplitude and the frequency of energy distribution can be further obtained, which provide conditions for directly observing the similarities and differences of paired signals.


## 2.3. The Ward Clustering Method

As a hierarchical agglomerate cluster algorithm, the WCM has a wide range of applications [19–21]. First, it is started by accepting each node as a separate cluster. Then, the clusters with minimum distance between themselves are combined in pairs at each stage of the algorithm. This smallest distance is called the Ward distance and defined as follows:(5)drs=nr⋅nsnr+nsx¯r−x¯s2,where r and s represent the two distinct clusters, nr and ns represent the number of data points of two clusters, respectively, x¯r and x¯s represent the center of the corresponding cluster, and ||⋅|| is Euclidean norm. The centers and cardinal numbers of the new cluster are updated according to the following equations:(6)x¯r′=nr⋅x¯r+ns⋅x¯snr+ns,(7)nr′=nr+ns.Ward’s clustering algorithm has the following steps:Step 1. Each sample point is treated as a cluster. At this time, the sum of squares of deviations for each cluster is equal to 0.
Step 2. Two arbitrary clusters are merged, and the sum of squares of deviations is calculated from Equations (5)–(7). If we assume that there are N clusters in total, it must be calculated NN−1/2 times.
Step 3. The two clusters with the smallest squared sum of deviations are combined into one class. The method eventually aggregates all sample points into one class when the number of clusters is unknown.If the number of clusters is known, the WCM can be directly used to classify the signal data after removing the noise. Otherwise, it can be obtained by analyzing the dendrogram of clustering. This is a rather subjective approach, which is difficult to help us finding the true number of clusters for a given data set. In recent research, the clustering validity indexes, such as Calinski-Harabasz (CH) index, Gap index, Silhouette (Silh) index, and Davies–Bouldin (DB) index, have been demonstrated to be the best validation tools for determining the optimal number of clusters [22–27].


### 2.3.1. Calinski-Harabasz Index

For a given setY=y1,y2,⋯,yN, assume that the dimension of each entity yi is v, i=1,2,⋯,N. K nonempty disjoint cluster sets S=S1,S2,⋯,SK around the centroid set C=c1,c2,⋯,cK can be obtained by minimizing the within-cluster distance WK:(8)WK=∑k=1K∑yi∈Skdyi,ck,where dyi,ck,k=1,2,⋯,K represents the squared Euclidean distance between the entity yi and the centroid ck, that is,(9)dyi,ck=∑j∈vyij−ckj2.Then, the CH index is defined as follows [24]:(10)CHK=T1−WK/K−1WK/N−K,where WK is defined as in (8), and T1 can be calculated by T1=∑i=1N∑j∈vyij−y¯j2.The CH index can reflect the compactness of the cluster by means of the overall within-cluster variance. The separation degree of the clusters can be reflected by the overall between-cluster variance. Therefore, a good clustering scheme corresponds to a higher value of CH index.


### 2.3.2. Silhouette Index

For each entityyi, its silhouette value measures the similarity between the entity yi and the points in its own cluster, when compared to the points in other clusters. This similarity is reflected by measuring the distance between the entity yi and the points derived from different clusters. The silhouette value of the entity yi is defined as follows [24]:(11)syi=byi−ayimaxayi,byi,where ayi is the average distance from the entity yi∈Sk to all other points yj∈Sk, byi is the minimum distance from the entity yi to all other points yj, which satisfies yj∈Sl,yi∉Sl,l≠k. Therefore, −1≤syi≤1. If syi is close to zero, the entity yi could be assigned to another cluster. A negative value of syi indicates that the corresponding assignment seriously damages cluster cohesion, and the clustering result of yi is not advisable. yi is well matched to its own cluster when syi is close to 1. Finally, the validity of the whole clustering can be quantified by Silh index, and it is defined as follows:(12)SilhK=∑i=1NsyiN.The Silh index can be used with any distance metric, including the Manhattan distances and Euclidean distances.


### 2.3.3. Davies–Bouldin Index

A good partition should have a larger intercluster separation degree and stronger within-cluster homogeneity and compactness. The DB index is proposed based on this idea [26]. More concretely, it is constructed by a ratio of within-cluster and between-cluster distances. The DB index is defined as follows:(13)DBK=∑k=1Kmaxk≠jd¯k+d¯jdkj/K,where(14)d¯k=1Sk∑i=1Skyi−ckq1/q,represents the average distance between each point yi in cluster k and the centroid of clusterk. Sk is the number of points in cluster k. If q=1, d¯k is the average Euclidean distance between the points in cluster k to the centroid of cluster k. If q=2, d¯k is the standard deviation of the distance of points in cluster k to the center of cluster k. When k in d¯k is replaced by j, d¯j can be obtained. In addition, dkj can be calculated according to the following equation:(15)dkj=∑h=1vc¯kh−c¯jhp1/p.It represents the distances between the centroids of thekth and the jth clusters. c¯kh is the hth component of the centroid of cluster k and dkj is the Minkowski metric of the centroids which characterizes clusters k and j. Specifically, p=1, dkj is the Manhattan distance between centroids p=2, dkj is the Euclidean distance between centroids.The DB index can reflect the degree of within-cluster dispersion and between-cluster separation. So, the true number of clusters may be determined according to the minimum value of the DB index.


### 2.3.4. Gap Index

Robert Tibshirani et al. proposed the gap statistic method for estimating the number of clusters in a set of data [27]. A graph of the within-cluster dispersion versus the number of clusters k for a clustering procedure shows that the within-cluster dispersion decreases monotonically as k increases, but from some k onwards, the decrease becomes flatter obviously. Such position is called ‘elbow’, it often implies the appropriate number of clusters. The gap criterion gives an approach to estimate the number of clusters by locating this ‘elbow’. Therefore, under this criterion, the optimal number of clusters occurs at the largest gap value. The Gap index is defined as follows:(16)GapNK=ENlogUK−logUK,where N represents the number of points, K represents the number of clusters that are evaluated, UK defined in (17) represents the within-cluster dispersion degree.(17)UK=∑k=1K12NkHk,where Nk is the number of points in cluster k, Hk is the sum of the distances of any two points in the kth cluster. The expected value ENlogUK is determined by Monte Carlo sampling from a reference distribution. The Gap index can also be used for any distance metric.The WCM belongs to unsupervised categorization technique, which can help us find the centroid of each cluster. However, the classification accuracy of this method is limited, which makes the method not able to be used directly for signal recognition. Comparatively, PNN can effectively improve the accuracy of classification and identification [28, 29].


## 2.3.1. Calinski-Harabasz Index

For a given setY=y1,y2,⋯,yN, assume that the dimension of each entity yi is v, i=1,2,⋯,N. K nonempty disjoint cluster sets S=S1,S2,⋯,SK around the centroid set C=c1,c2,⋯,cK can be obtained by minimizing the within-cluster distance WK:(8)WK=∑k=1K∑yi∈Skdyi,ck,where dyi,ck,k=1,2,⋯,K represents the squared Euclidean distance between the entity yi and the centroid ck, that is,(9)dyi,ck=∑j∈vyij−ckj2.Then, the CH index is defined as follows [24]:(10)CHK=T1−WK/K−1WK/N−K,where WK is defined as in (8), and T1 can be calculated by T1=∑i=1N∑j∈vyij−y¯j2.The CH index can reflect the compactness of the cluster by means of the overall within-cluster variance. The separation degree of the clusters can be reflected by the overall between-cluster variance. Therefore, a good clustering scheme corresponds to a higher value of CH index.


## 2.3.2. Silhouette Index

For each entityyi, its silhouette value measures the similarity between the entity yi and the points in its own cluster, when compared to the points in other clusters. This similarity is reflected by measuring the distance between the entity yi and the points derived from different clusters. The silhouette value of the entity yi is defined as follows [24]:(11)syi=byi−ayimaxayi,byi,where ayi is the average distance from the entity yi∈Sk to all other points yj∈Sk, byi is the minimum distance from the entity yi to all other points yj, which satisfies yj∈Sl,yi∉Sl,l≠k. Therefore, −1≤syi≤1. If syi is close to zero, the entity yi could be assigned to another cluster. A negative value of syi indicates that the corresponding assignment seriously damages cluster cohesion, and the clustering result of yi is not advisable. yi is well matched to its own cluster when syi is close to 1. Finally, the validity of the whole clustering can be quantified by Silh index, and it is defined as follows:(12)SilhK=∑i=1NsyiN.The Silh index can be used with any distance metric, including the Manhattan distances and Euclidean distances.


## 2.3.3. Davies–Bouldin Index

A good partition should have a larger intercluster separation degree and stronger within-cluster homogeneity and compactness. The DB index is proposed based on this idea [26]. More concretely, it is constructed by a ratio of within-cluster and between-cluster distances. The DB index is defined as follows:(13)DBK=∑k=1Kmaxk≠jd¯k+d¯jdkj/K,where(14)d¯k=1Sk∑i=1Skyi−ckq1/q,represents the average distance between each point yi in cluster k and the centroid of clusterk. Sk is the number of points in cluster k. If q=1, d¯k is the average Euclidean distance between the points in cluster k to the centroid of cluster k. If q=2, d¯k is the standard deviation of the distance of points in cluster k to the center of cluster k. When k in d¯k is replaced by j, d¯j can be obtained. In addition, dkj can be calculated according to the following equation:(15)dkj=∑h=1vc¯kh−c¯jhp1/p.It represents the distances between the centroids of thekth and the jth clusters. c¯kh is the hth component of the centroid of cluster k and dkj is the Minkowski metric of the centroids which characterizes clusters k and j. Specifically, p=1, dkj is the Manhattan distance between centroids p=2, dkj is the Euclidean distance between centroids.The DB index can reflect the degree of within-cluster dispersion and between-cluster separation. So, the true number of clusters may be determined according to the minimum value of the DB index.


## 2.3.4. Gap Index

Robert Tibshirani et al. proposed the gap statistic method for estimating the number of clusters in a set of data [27]. A graph of the within-cluster dispersion versus the number of clusters k for a clustering procedure shows that the within-cluster dispersion decreases monotonically as k increases, but from some k onwards, the decrease becomes flatter obviously. Such position is called ‘elbow’, it often implies the appropriate number of clusters. The gap criterion gives an approach to estimate the number of clusters by locating this ‘elbow’. Therefore, under this criterion, the optimal number of clusters occurs at the largest gap value. The Gap index is defined as follows:(16)GapNK=ENlogUK−logUK,where N represents the number of points, K represents the number of clusters that are evaluated, UK defined in (17) represents the within-cluster dispersion degree.(17)UK=∑k=1K12NkHk,where Nk is the number of points in cluster k, Hk is the sum of the distances of any two points in the kth cluster. The expected value ENlogUK is determined by Monte Carlo sampling from a reference distribution. The Gap index can also be used for any distance metric.The WCM belongs to unsupervised categorization technique, which can help us find the centroid of each cluster. However, the classification accuracy of this method is limited, which makes the method not able to be used directly for signal recognition. Comparatively, PNN can effectively improve the accuracy of classification and identification [28, 29].


## 2.4. Probabilistic Neural Network Classifier

As a method of nonparametric Parzen windows estimation, PNN is first proposed by Specht. It is a nonlinear classification technique and essentially a parallel algorithm based on Bayesian minimum risk criterion [30]. Given a sample to be identified x, its posterior probability PSkx can be obtained by PNN classifier. However, if the probability densities of the classes to be separated are unknown, the training samples with known identity need to be used to estimate them. Finally, the trained PNN is used to determine the identity of x. A typical PNN classifier consists of an input layer, a pattern layer (hidden layer), a summation layer and a output layer. The flowchart of the PNN is shown in Figure 2.Figure 2
Flowchart of the PNN.The input layer neurons are used to receive values from training samples and send data to the neurons in the pattern layer, which is fully connected to the input layer. The number of neurons in the input layer is equal to the length of the input vector. The number of neurons in the pattern layer is the same as the number of training samples. Here, all neurons are collected into different groups, and theith neuron in group k corresponds to a Gaussian function fikx,σ,i=1,2,⋯,mk, where mk represents the number of neurons in group k, k=1,2,⋯,K. Gaussian function which is also called probability density function is defined as follows:(18)fikx,σ=12πv/2σvexp−∑j=1vxijk−xj22σ2,where v is the dimension of the input vector x=x1,x2,⋯,xv, xj is the jth component of the input vector x, xijk is the jth component of the ith neuron in class k. The so-called smoothing parameter σ∈0,1 determined experimentally by comparing their corresponding classification accuracy plays an important role in estimation error of the PNN classifier. The outputs of pattern layer are connected to the summation units depending on the class of patterns. There is one neuron for each group, and each neuron in summation layer sums the outputs derived from the pattern layer neurons as follows:(19)pkx=12πv/2σvmk∑i=1mkexp−∑j=1vxijk−xj22σ2.Finally, the output layer neuron output a number 1 and multiple numbers 0. The value of 1 corresponds to the classifier’s decision result for input vectors. More specifically, the input vectorx belongs to class k if pkx>pk′x for all k′=1,2,⋯,K and k≠k′.Hence, the main purpose of training PNN is to find the optimal estimate of probability density function according to the training samples and their labels, to ensure that the classifier works at the condition of minimum error rate and risk. When the samples to be identified are sent to the pattern layer, the output of each neuron is calculated according to the trained density function. Finally, the identified results are obtained through computations in the summation layer and output layer. Due to the following advantages, it is a wise choice to use PNN as a further classifier to classify signals [31]:(1)
It has a simple structure, and it is easy to train. In the PNN based on probability density function estimation, the weight of the neuron in pattern is directly taken from the input sample value.(2)
The training process of the network is simple, and there is no need to retrain for a long time when adding or reducing the number of groups.(3)
It is not easy to produce local optimal solution, and its precision is higher than that of other classification approaches. No matter how complex the classification problem is, as long as there are enough training samples, the optimal solution under the Bayes criterion can be obtained.


## 3. Classification Scheme and Experiments of Emitter Signals

### 3.1. Flowchart and Algorithms of Classification Scheme

The flowchart of the proposed classification algorithms is shown in Figure3.Figure 3
Flowchart of signal classification.Figure3 indicates that the proposed scheme is composed of four modules, that is, data processing module, preclassification module, evaluation module, and accurate classification module. In the evaluation module, the clustering validity indexes are used to determine the range of K if it is unknown. For each K∈KMin,KMax, the classification validity index D is calculated as follows [16]:(20)DK=K∑j=1N1max1≤k≤Kpkj−N1N1K−1,where N1 is the number of input vectors, pkj is the element of the matrix Q of size K×N1 in the output of PNN’s pattern layer representing the membership of the jth input vector to the cluster k. When N1=1, pk1 is equal to pkx presented in (19). When N1>1, matrix Q=pkjK×N1 can be obtained by PNN. max1≤k≤Kpkj is the largest element of the jth column in the matrix Q. Equation (20) indicates that DK is a nonlinear function related to K, when N1≠∑j=1N1max1≤k≤Kpkj. Therefore, K corresponding to the maximum value of DK is the optimal number of clusters K∗.The pseudocodes are listed in Algorithm1 if K is known.Algorithm 1:Classification Algorithm 1.
Input: The original signal vectors yi that need to be classified.
The classification numberK.
Output: The classified label vector Cy.
Computesi by using self-adaptive filtering for yi;
Compute frequency spectrumGi of si;
Compute centerck of each class k by using WCM;
Select training samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples.
End: Classification Algorithm 1.If the classification numberK is unknown, WCM and clustering validity indexes are used to determine the range of K. The corresponding pseudocodes are listed in Algorithm 2.Algorithm 2:Classification Algorithm 2.
Input: The original signal vectors yi that need to be classified.
Output: The classified label vector Cy.
Computesi by using self-adaptive filtering for yi;
Compute frequency spectrumGi of si;
Compute Ward’s clustering dendrogram;
ComputeCH(K), Silh(K), DB(K) and Gap(K); Compute Kmin, Kmax;
ifKmin = Kmax, then
Compute centerck of each class k;
Select training samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples;
Output the classified label ofyi;
else forK = Kmin : Kmax
fork = 1 : K
Compute centerck of each class k;
Select samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples;
Compute matrixQ=pkj and DK for each K;
ComputeDmax and the optimal number of classifications K∗;
Output the classified label ofyi;
end
end
End: Classification Algorithm 2.The algorithms show that the supervised learning PNN is used to classify samples. Therefore, the training (teaching) samples must be selected first. By Ward clustering method, we have obtained a preliminary classification of all samples. That is, the identities of some samples have been determined, except for some boundary points, which need to be further determined by the trained PNN. Therefore, some labeled proximity pointsxkjk=1,2,⋯,K,j=1,2,⋯,Jk around the center ck can be selected to train PNN, where Jk represents the number of samples selected in class k and it should be preset, such as Jk=a⋅Sk, a∈0.6,0.8.


### 3.2. Classification Experiments

A signal setRs sampled from some pulse emitters is used to test the effectiveness of the proposed algorithms. Each emitter emits continuous signals in the pulse state. After a period of time, the receiver will receive multiple signals from all emitters. These signals are converted into digital signals by the analog-digital converter. The sampling frequency is 1.01 MHz. Signal samples yii=1,2,⋯,500 are randomly extracted from these digital signals. The signal set Rs=y1,y2,⋯,y500T and the dimension of yi is 1024. Considering that each yi is disturbed by signals from other emitters, the mean value of all signals is used as the noise signal Is=Is1Is2⋯Is1024, where Isn=∑i=1500yin/500, n=1,2,⋯,1024. First, the self-adaptive filtering is applied to process these signals. In this algorithm, the actual signal xn corresponds to yin and the noise signal x0n corresponds to Isn. Then, Fourier transform is used to obtain the amplitude spectrum of all processed signals. The amplitude spectrum of the thirteenth signal is shown in Figure 4.Figure 4
The double-sided amplitude spectrum of a sampling signal. (a) Before filtering; (b) after filtering.
(a)
(b)Figure4(a) shows the sampling signal contains obvious white noise, which makes the feature of the useful signal unclear. However, most of the noises are removed after using the self-adaptive filtering. Moreover, the characteristics of the signals are highlighted so that the amplitude spectrum of them can be analyzed correctly. For these transformed signals, the clustering dendrogram can be obtained by using the WCM, and it is shown in Figure 5.Figure 5
The dendrogram of the signals after using the WCM.When the signals are divided into 3, 4, 5, and 6 classes, the intercluster distances are 75.5, 57.5, 27, and 14, respectively. The increments of the distance between the clusters are 13, 30.5, and 18 in turn. If the number of elements in each class is required relatively close, the ideal number of classifications is 3, 4, and 5. However, these numbers need to be further determined by the clustering validity indexes. The number of clustersK corresponding to different evaluation indexes is shown in Figure 6.Figure 6
The number of clustersK versus the different evaluation indexes. (a) The number of clusters K versus CH values. (b) The number of clusters K versus DB values. (c) The number of clusters K versus Silh values. (d) The number of clusters K versus Gap values.
(a)
(b)
(c)
(d)Figure6 shows that the optimal number of clusters is 3 when the DB index and the Silh index are used, and it is 2 when the CH index is used. However, it becomes to 5 when the Gap index is used. Therefore, the optimal number of classifications should belong to the interval [2, 5]. In this case, the PNN needs to be used to obtain more accurate results.For eachK in this interval, seventy (a=0.7) sample points nearby each ck are selected to train PNN classifiers. Then, the optimal classifier is obtained by calculating the maximum value of DK. The results have shown that D5=1, while the other values are less than D(5). That is to say, the optimal number of classifications is 5. Since the size of matrix Rs is 500 × 1024, three columns are randomly selected as the X-axis, the Y-axis, and the Z-axis to plot the classification results diagram. Let(21)Ack=y1k1y1k2⋯y1k1024y2k1y2k2⋯y2k1024⋮⋮⋮yMkk1yMkk2⋯yMkk1024,represents a matrix consisting of all signals in class k after classifying, where Mk is the number of samples in class k, yjk is the jth sample in class k, and j=1,2,⋯,Mk. Thus, Ac=Ac1Ac2⋯Ac5T is a matrix of size 500 × 1024. If the data on columns a1,a2,a3 of Ac are selected to form a matrix Ac¯ of size 500 × 3, each row of Ac¯ is a three-dimensional array, that is, a point in the coordinate system. When a=0.5,σ=1, scatter plots of these data are shown in Figure 7. The first column of Ac¯ corresponds to the X-axis, the second column of Ac¯ corresponds to the Y-axis, and the third column corresponds to the Z-axis.Figure 7
The scatter plots of classification results. (a) The classified results whena1=245,a2=562,a3=720. (b) The classified results when a1=268,a2=369,a3=547. (c) The classified results when a1=139,a2=367,a3=741.
(a)
(b)
(c)The experimental results show that all signal samples are divided into five classes, and each class contains 100 signals. Therefore, all signals can be thought to come from five emitters, and each emitter emits 100 signals. Although only three distribution figures of all classified signals are presented in Figure7, in fact, in our experiments, we have obtained more than 1000 scatter plots which are drawn by randomly selecting three columns from matrix Ac. In these classification results, the distributions of sample sets are similar and the separations of them are obvious. Therefore, the proposed methods can effectively distinguish signals from different emitters.


## 3.1. Flowchart and Algorithms of Classification Scheme

The flowchart of the proposed classification algorithms is shown in Figure3.Figure 3
Flowchart of signal classification.Figure3 indicates that the proposed scheme is composed of four modules, that is, data processing module, preclassification module, evaluation module, and accurate classification module. In the evaluation module, the clustering validity indexes are used to determine the range of K if it is unknown. For each K∈KMin,KMax, the classification validity index D is calculated as follows [16]:(20)DK=K∑j=1N1max1≤k≤Kpkj−N1N1K−1,where N1 is the number of input vectors, pkj is the element of the matrix Q of size K×N1 in the output of PNN’s pattern layer representing the membership of the jth input vector to the cluster k. When N1=1, pk1 is equal to pkx presented in (19). When N1>1, matrix Q=pkjK×N1 can be obtained by PNN. max1≤k≤Kpkj is the largest element of the jth column in the matrix Q. Equation (20) indicates that DK is a nonlinear function related to K, when N1≠∑j=1N1max1≤k≤Kpkj. Therefore, K corresponding to the maximum value of DK is the optimal number of clusters K∗.The pseudocodes are listed in Algorithm1 if K is known.Algorithm 1:Classification Algorithm 1.
Input: The original signal vectors yi that need to be classified.
The classification numberK.
Output: The classified label vector Cy.
Computesi by using self-adaptive filtering for yi;
Compute frequency spectrumGi of si;
Compute centerck of each class k by using WCM;
Select training samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples.
End: Classification Algorithm 1.If the classification numberK is unknown, WCM and clustering validity indexes are used to determine the range of K. The corresponding pseudocodes are listed in Algorithm 2.Algorithm 2:Classification Algorithm 2.
Input: The original signal vectors yi that need to be classified.
Output: The classified label vector Cy.
Computesi by using self-adaptive filtering for yi;
Compute frequency spectrumGi of si;
Compute Ward’s clustering dendrogram;
ComputeCH(K), Silh(K), DB(K) and Gap(K); Compute Kmin, Kmax;
ifKmin = Kmax, then
Compute centerck of each class k;
Select training samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples;
Output the classified label ofyi;
else forK = Kmin : Kmax
fork = 1 : K
Compute centerck of each class k;
Select samplesxkj around ck and record their labels Lxkj;
Create the PNN classifier by usingxkj and Lxkj;
Determine the class of the remaining samples;
Compute matrixQ=pkj and DK for each K;
ComputeDmax and the optimal number of classifications K∗;
Output the classified label ofyi;
end
end
End: Classification Algorithm 2.The algorithms show that the supervised learning PNN is used to classify samples. Therefore, the training (teaching) samples must be selected first. By Ward clustering method, we have obtained a preliminary classification of all samples. That is, the identities of some samples have been determined, except for some boundary points, which need to be further determined by the trained PNN. Therefore, some labeled proximity pointsxkjk=1,2,⋯,K,j=1,2,⋯,Jk around the center ck can be selected to train PNN, where Jk represents the number of samples selected in class k and it should be preset, such as Jk=a⋅Sk, a∈0.6,0.8.


## 3.2. Classification Experiments

A signal setRs sampled from some pulse emitters is used to test the effectiveness of the proposed algorithms. Each emitter emits continuous signals in the pulse state. After a period of time, the receiver will receive multiple signals from all emitters. These signals are converted into digital signals by the analog-digital converter. The sampling frequency is 1.01 MHz. Signal samples yii=1,2,⋯,500 are randomly extracted from these digital signals. The signal set Rs=y1,y2,⋯,y500T and the dimension of yi is 1024. Considering that each yi is disturbed by signals from other emitters, the mean value of all signals is used as the noise signal Is=Is1Is2⋯Is1024, where Isn=∑i=1500yin/500, n=1,2,⋯,1024. First, the self-adaptive filtering is applied to process these signals. In this algorithm, the actual signal xn corresponds to yin and the noise signal x0n corresponds to Isn. Then, Fourier transform is used to obtain the amplitude spectrum of all processed signals. The amplitude spectrum of the thirteenth signal is shown in Figure 4.Figure 4
The double-sided amplitude spectrum of a sampling signal. (a) Before filtering; (b) after filtering.
(a)
(b)Figure4(a) shows the sampling signal contains obvious white noise, which makes the feature of the useful signal unclear. However, most of the noises are removed after using the self-adaptive filtering. Moreover, the characteristics of the signals are highlighted so that the amplitude spectrum of them can be analyzed correctly. For these transformed signals, the clustering dendrogram can be obtained by using the WCM, and it is shown in Figure 5.Figure 5
The dendrogram of the signals after using the WCM.When the signals are divided into 3, 4, 5, and 6 classes, the intercluster distances are 75.5, 57.5, 27, and 14, respectively. The increments of the distance between the clusters are 13, 30.5, and 18 in turn. If the number of elements in each class is required relatively close, the ideal number of classifications is 3, 4, and 5. However, these numbers need to be further determined by the clustering validity indexes. The number of clustersK corresponding to different evaluation indexes is shown in Figure 6.Figure 6
The number of clustersK versus the different evaluation indexes. (a) The number of clusters K versus CH values. (b) The number of clusters K versus DB values. (c) The number of clusters K versus Silh values. (d) The number of clusters K versus Gap values.
(a)
(b)
(c)
(d)Figure6 shows that the optimal number of clusters is 3 when the DB index and the Silh index are used, and it is 2 when the CH index is used. However, it becomes to 5 when the Gap index is used. Therefore, the optimal number of classifications should belong to the interval [2, 5]. In this case, the PNN needs to be used to obtain more accurate results.For eachK in this interval, seventy (a=0.7) sample points nearby each ck are selected to train PNN classifiers. Then, the optimal classifier is obtained by calculating the maximum value of DK. The results have shown that D5=1, while the other values are less than D(5). That is to say, the optimal number of classifications is 5. Since the size of matrix Rs is 500 × 1024, three columns are randomly selected as the X-axis, the Y-axis, and the Z-axis to plot the classification results diagram. Let(21)Ack=y1k1y1k2⋯y1k1024y2k1y2k2⋯y2k1024⋮⋮⋮yMkk1yMkk2⋯yMkk1024,represents a matrix consisting of all signals in class k after classifying, where Mk is the number of samples in class k, yjk is the jth sample in class k, and j=1,2,⋯,Mk. Thus, Ac=Ac1Ac2⋯Ac5T is a matrix of size 500 × 1024. If the data on columns a1,a2,a3 of Ac are selected to form a matrix Ac¯ of size 500 × 3, each row of Ac¯ is a three-dimensional array, that is, a point in the coordinate system. When a=0.5,σ=1, scatter plots of these data are shown in Figure 7. The first column of Ac¯ corresponds to the X-axis, the second column of Ac¯ corresponds to the Y-axis, and the third column corresponds to the Z-axis.Figure 7
The scatter plots of classification results. (a) The classified results whena1=245,a2=562,a3=720. (b) The classified results when a1=268,a2=369,a3=547. (c) The classified results when a1=139,a2=367,a3=741.
(a)
(b)
(c)The experimental results show that all signal samples are divided into five classes, and each class contains 100 signals. Therefore, all signals can be thought to come from five emitters, and each emitter emits 100 signals. Although only three distribution figures of all classified signals are presented in Figure7, in fact, in our experiments, we have obtained more than 1000 scatter plots which are drawn by randomly selecting three columns from matrix Ac. In these classification results, the distributions of sample sets are similar and the separations of them are obvious. Therefore, the proposed methods can effectively distinguish signals from different emitters.


## 4. Identification Scheme and Experiments of Emitter Signals

### 4.1. Flowchart and Algorithms of Identification Scheme

Besides classification, data identification is an important function of PNN. Since PNN is based on the maximum posterior probability, it will give the optimal solution under the Bayesian criterion, whether or not the samples to be identifiedIxi belong to the five determined classes that have been obtained in Section 3. So, if a sample belongs to one of them, it will accurately identify it. However, it may lose its function when the sample do not belong to them. At this time, the amplitude spectrum of all samples can be analyzed before, so that they can be identified whether they belong to the determined classes in advance. Considering that it is difficult to comparatively analyze the amplitude spectrum of the signals to be identified and each signal in every class, we adopt the curve fitting method to find the feature sequences of each class and the signals to be identified. Finally, the correlation degree of these sequences is calculated to get the preliminary identity information of the samples to be identified. This method is called Bivariable correlation analysis, and it is introduced as follows.Step 1. Simplifying amplitude spectrum. Let Li represent the length of the signal Ixi, Fsi represent the Fourier transform of Ixi, Fsi is the amplitude spectrum of Ixi, Fsi=gi1gi2⋯giLi. The sequence Fi of length N2 is extracted from Fsi which takes ΔTi as the step size, where ΔTi=Li/N2. For each training sample xkl in class k, the corresponding sequence Fkl can be obtained according to the same method, where k=1,2,⋯,K∗,l=1,2,⋯,L,i=1,2,⋯,p, and L represents the number of training samples in class k, and p represents the number of samples to be identified.Step 2. Fitting curve. The fitting curve z^i can be obtained by fitting the sequence Fi and the fitting curve y^kl in class k can be obtained by fitting the sequence Fkl.Step 3. Constructing the feature sequence. First, for different k, y^klt,t=1,2,⋯,T can be calculated after giving the upper bound T. Then, the tth signal feature of the class k can be obtained by the following equation:(22)Y^kt=1L∑l=1Ly^klt.Finally,Csk=Y^k1,Y^k2,⋯,Y^kT is the feature sequence of the class k. Similarly, the feature sequences Zi=Z^i1,Z^i2,⋯,Z^iT,i=1,2,⋯,p of the samples to be identified can be obtained.Step 4. Correlation analysis. The correlation coefficient between the Zi and the Csk can be calculated as the following equation:(23)rik=T∑t=1TY^kt⋅Z^it−∑t=1TY^kt⋅∑t=1TZ^itT∑t=1TY^k2t−∑t=1TY^kt2·T∑t=1TZ^i2t−∑t=1TZ^it2,where T∑t=1TY^kt⋅Z^it−∑t=1TY^kt⋅∑t=1TZ^it represents the sample covariance of Csk and Zi, T∑t=1TY^k2t−∑t=1TY^kt2 and T∑t=1TZ^i2t−∑t=1TZ^it2 are the sample standard variation of Csk and Zi, respectively. The performance of the correlation test indicates that the sample to be identified belongs to the corresponding class when r≥0.95. All samples that do not satisfy this condition should be removed. Finally, the remaining samples can be effectively identified by the trained PNN. The flowchart of the proposed identification scheme is shown in Figure 8.Figure 8
Flowchart of signal identification.Figure8 shows that the flowchart mainly includes two blocks. One is the preidentification module, and the other is the identification module. The role of the previous module is to eliminate the sample with a small correlation to the determined classes. When the correlations between the samples to be identified and several different classes are high, PNN will accurately identify it based on Bayesian criteria. The pseudocodes are given in Algorithm 3, which is called identification algorithm.Algorithm 3:Identification algorithm.
Input: The training samples xkl, label matrix Bxkl.
Samples to be identifiedIxi.
Output: Ixi is the element of the class k.
Normalizexkl, Ixi;
Create PNN model by usingxkl and Bxkl;
Adjust the dimension ofIxi to make it consistent with the xkl;
Determining sequencesFi and Fkl;
Compute fitting curvesz^i, y^kl;
Create feature sequencesCsk,Zi;
Computerik, and if rik≥0.95, then
Put theIxi into the trained PNN;
else
Ixi does not belong to any class;
End: Identification algorithm.Through the methods proposed in Section3, all signal samples have been classified to five classes. That is to say, the identity of each sample has been determined. Thus, the identification algorithm consists of two parts. On the one hand, some proximity points xkl around the center ck of class k are selected to train PNN, k=1,2,⋯,K∗,l=1,2,⋯,L. xkl is the input vector of the network, and its label is the output of the system. On the other hand, the trained PNN classifier is used to determine the identity of Ixi,i=1,2,⋯,p.If the dimension ofIxi is greater than xkl, the method of adding time window can be used to adjust their dimensions. Let Dx represent the dimension of xkl and DIx represent the dimension of Ixi. The time window Ts=DIx/Dx is used to reduce the value of DIx. Thus, the adjusted samples to be identified I¯xi=Ixi1,Ixi1+Ts,⋯,Ixi1+Dx−1Ts. When DIx<Dx, the same method can be used to reduce the dimension of xkl and make it consistent with DIx. At this moment, DIx and Dx need to be exchanged.


### 4.2. Identification Experiments

By adjusting some parameters of emitters, such as pulse width and output power, signals that are different fromyii=1,2,⋯,500 can be obtained. Suppose that sample Ixii=1,2,3,4 is randomly selected from these signals for identification and the dimension of Ixi is 10240. xkl represents the lth training sample in class k, Bxkl is the label matrix. First, I¯xi=Ixi1,Ixi1+10,⋯,Ixi1+10230 can be obtained according to the method of adding time window. Then, the single-sided amplitude spectrum of center ck for class k is shown in Figure 9; k=1,2,⋯,5. The double-sided amplitude spectrum of I¯xii=1,2,3,4 is shown in Figure 10. Finally, the identity of I¯xii=1,2,3,4 is preliminary determined by observing the similarity between these amplitude spectrum.Figure 9
The single-sided amplitude spectrum of classk,k=1,2,⋯,5. (a) k=1; (b) k=2; (c) k=3; (d) k=4; (e) k=5.
(a)
(b)
(c)
(d)
(e)Figure 10
The double-sided amplitude spectrum ofI¯xii=1,2,3,4. (a) i=1; (b) i=2; (c) i=3; (d) i=4.
(a)
(b)
(c)
(d)It is observed that the amplitude spectrum ofI¯x1 is similar to that of the signals in class 2 or class 3, and I¯x2 is likely to belong to class 5. But, the class of I¯x3 and I¯x4 is difficult to determine only by observation. At this time, it is necessary to carry out correlation test for these samples. The test results are shown in Table 1 when N2=300 and T=100.Table 1
Correlation test of the samples to be identified.
r
i
k
Class 1
Class 2
Class 3
Class 4
Class 5
I
¯
x
1
0.3687
0.9704
0.9945
0.5473
0.6910
I
¯
x
2
0.7991
0.7447
0.6674
0.7882
0.9941
I
¯
x
3
0.2364
0.8907
0.7982
0.5202
0.5583
I
¯
x
4
0.7864
0.3572
0.3083
0.6887
0.6841Table1 shows that r12=0.9704 and r13=0.9945, which are all greater than 0.95. So, I¯x1 might be the element of class 2 or class 3. I¯x2 should be classified into class 5. These results are consistent with the observations. I¯x3 and I¯x4 should be removed before formal classification. Therefore, it is necessary to identify I¯x1 and I¯x2 by using PNN. L points are selected as training samples in each class to train PNN, L=70. Finally, The trained PNN is used to judge the identity of I¯x1 and I¯x2. Since Ix1 corresponds to I¯x1, Ix2 corresponds to I¯x2, we can obtained the identity of Ix1 and Ix2 by identifying I¯x1 and I¯x2. The identified results are shown in Figure 11 when σ=1.Figure 11
The identified results ofIxia1Ixia2Ixia3,i=1,2. (a) a1=41,a2=256,a3=682andi=1; (b) a1=35,a2=709,a3=929andi=2.
(a)
(b)In order to obtain intuitive results, the data on columns 41, 256, 682 ofIx1 are selected to form a three-dimensional array Ix141Ix1256Ix1682. Data on the same columns of Ac are selected to form a matrix of size 500×3. All data are drawn in Figure 11(a). Obviously, Ix1 should be the element in class 3 (black solid ball) and Ix2 should be the element in class 5 (red solid ball). The same results can also be obtained when the data on columns 35, 709, and 929 of Ix1 and Ac are selected. The corresponding results are drawn in Figure 11(b).


## 4.1. Flowchart and Algorithms of Identification Scheme

Besides classification, data identification is an important function of PNN. Since PNN is based on the maximum posterior probability, it will give the optimal solution under the Bayesian criterion, whether or not the samples to be identifiedIxi belong to the five determined classes that have been obtained in Section 3. So, if a sample belongs to one of them, it will accurately identify it. However, it may lose its function when the sample do not belong to them. At this time, the amplitude spectrum of all samples can be analyzed before, so that they can be identified whether they belong to the determined classes in advance. Considering that it is difficult to comparatively analyze the amplitude spectrum of the signals to be identified and each signal in every class, we adopt the curve fitting method to find the feature sequences of each class and the signals to be identified. Finally, the correlation degree of these sequences is calculated to get the preliminary identity information of the samples to be identified. This method is called Bivariable correlation analysis, and it is introduced as follows.Step 1. Simplifying amplitude spectrum. Let Li represent the length of the signal Ixi, Fsi represent the Fourier transform of Ixi, Fsi is the amplitude spectrum of Ixi, Fsi=gi1gi2⋯giLi. The sequence Fi of length N2 is extracted from Fsi which takes ΔTi as the step size, where ΔTi=Li/N2. For each training sample xkl in class k, the corresponding sequence Fkl can be obtained according to the same method, where k=1,2,⋯,K∗,l=1,2,⋯,L,i=1,2,⋯,p, and L represents the number of training samples in class k, and p represents the number of samples to be identified.Step 2. Fitting curve. The fitting curve z^i can be obtained by fitting the sequence Fi and the fitting curve y^kl in class k can be obtained by fitting the sequence Fkl.Step 3. Constructing the feature sequence. First, for different k, y^klt,t=1,2,⋯,T can be calculated after giving the upper bound T. Then, the tth signal feature of the class k can be obtained by the following equation:(22)Y^kt=1L∑l=1Ly^klt.Finally,Csk=Y^k1,Y^k2,⋯,Y^kT is the feature sequence of the class k. Similarly, the feature sequences Zi=Z^i1,Z^i2,⋯,Z^iT,i=1,2,⋯,p of the samples to be identified can be obtained.Step 4. Correlation analysis. The correlation coefficient between the Zi and the Csk can be calculated as the following equation:(23)rik=T∑t=1TY^kt⋅Z^it−∑t=1TY^kt⋅∑t=1TZ^itT∑t=1TY^k2t−∑t=1TY^kt2·T∑t=1TZ^i2t−∑t=1TZ^it2,where T∑t=1TY^kt⋅Z^it−∑t=1TY^kt⋅∑t=1TZ^it represents the sample covariance of Csk and Zi, T∑t=1TY^k2t−∑t=1TY^kt2 and T∑t=1TZ^i2t−∑t=1TZ^it2 are the sample standard variation of Csk and Zi, respectively. The performance of the correlation test indicates that the sample to be identified belongs to the corresponding class when r≥0.95. All samples that do not satisfy this condition should be removed. Finally, the remaining samples can be effectively identified by the trained PNN. The flowchart of the proposed identification scheme is shown in Figure 8.Figure 8
Flowchart of signal identification.Figure8 shows that the flowchart mainly includes two blocks. One is the preidentification module, and the other is the identification module. The role of the previous module is to eliminate the sample with a small correlation to the determined classes. When the correlations between the samples to be identified and several different classes are high, PNN will accurately identify it based on Bayesian criteria. The pseudocodes are given in Algorithm 3, which is called identification algorithm.Algorithm 3:Identification algorithm.
Input: The training samples xkl, label matrix Bxkl.
Samples to be identifiedIxi.
Output: Ixi is the element of the class k.
Normalizexkl, Ixi;
Create PNN model by usingxkl and Bxkl;
Adjust the dimension ofIxi to make it consistent with the xkl;
Determining sequencesFi and Fkl;
Compute fitting curvesz^i, y^kl;
Create feature sequencesCsk,Zi;
Computerik, and if rik≥0.95, then
Put theIxi into the trained PNN;
else
Ixi does not belong to any class;
End: Identification algorithm.Through the methods proposed in Section3, all signal samples have been classified to five classes. That is to say, the identity of each sample has been determined. Thus, the identification algorithm consists of two parts. On the one hand, some proximity points xkl around the center ck of class k are selected to train PNN, k=1,2,⋯,K∗,l=1,2,⋯,L. xkl is the input vector of the network, and its label is the output of the system. On the other hand, the trained PNN classifier is used to determine the identity of Ixi,i=1,2,⋯,p.If the dimension ofIxi is greater than xkl, the method of adding time window can be used to adjust their dimensions. Let Dx represent the dimension of xkl and DIx represent the dimension of Ixi. The time window Ts=DIx/Dx is used to reduce the value of DIx. Thus, the adjusted samples to be identified I¯xi=Ixi1,Ixi1+Ts,⋯,Ixi1+Dx−1Ts. When DIx<Dx, the same method can be used to reduce the dimension of xkl and make it consistent with DIx. At this moment, DIx and Dx need to be exchanged.


## 4.2. Identification Experiments

By adjusting some parameters of emitters, such as pulse width and output power, signals that are different fromyii=1,2,⋯,500 can be obtained. Suppose that sample Ixii=1,2,3,4 is randomly selected from these signals for identification and the dimension of Ixi is 10240. xkl represents the lth training sample in class k, Bxkl is the label matrix. First, I¯xi=Ixi1,Ixi1+10,⋯,Ixi1+10230 can be obtained according to the method of adding time window. Then, the single-sided amplitude spectrum of center ck for class k is shown in Figure 9; k=1,2,⋯,5. The double-sided amplitude spectrum of I¯xii=1,2,3,4 is shown in Figure 10. Finally, the identity of I¯xii=1,2,3,4 is preliminary determined by observing the similarity between these amplitude spectrum.Figure 9
The single-sided amplitude spectrum of classk,k=1,2,⋯,5. (a) k=1; (b) k=2; (c) k=3; (d) k=4; (e) k=5.
(a)
(b)
(c)
(d)
(e)Figure 10
The double-sided amplitude spectrum ofI¯xii=1,2,3,4. (a) i=1; (b) i=2; (c) i=3; (d) i=4.
(a)
(b)
(c)
(d)It is observed that the amplitude spectrum ofI¯x1 is similar to that of the signals in class 2 or class 3, and I¯x2 is likely to belong to class 5. But, the class of I¯x3 and I¯x4 is difficult to determine only by observation. At this time, it is necessary to carry out correlation test for these samples. The test results are shown in Table 1 when N2=300 and T=100.Table 1
Correlation test of the samples to be identified.
r
i
k
Class 1
Class 2
Class 3
Class 4
Class 5
I
¯
x
1
0.3687
0.9704
0.9945
0.5473
0.6910
I
¯
x
2
0.7991
0.7447
0.6674
0.7882
0.9941
I
¯
x
3
0.2364
0.8907
0.7982
0.5202
0.5583
I
¯
x
4
0.7864
0.3572
0.3083
0.6887
0.6841Table1 shows that r12=0.9704 and r13=0.9945, which are all greater than 0.95. So, I¯x1 might be the element of class 2 or class 3. I¯x2 should be classified into class 5. These results are consistent with the observations. I¯x3 and I¯x4 should be removed before formal classification. Therefore, it is necessary to identify I¯x1 and I¯x2 by using PNN. L points are selected as training samples in each class to train PNN, L=70. Finally, The trained PNN is used to judge the identity of I¯x1 and I¯x2. Since Ix1 corresponds to I¯x1, Ix2 corresponds to I¯x2, we can obtained the identity of Ix1 and Ix2 by identifying I¯x1 and I¯x2. The identified results are shown in Figure 11 when σ=1.Figure 11
The identified results ofIxia1Ixia2Ixia3,i=1,2. (a) a1=41,a2=256,a3=682andi=1; (b) a1=35,a2=709,a3=929andi=2.
(a)
(b)In order to obtain intuitive results, the data on columns 41, 256, 682 ofIx1 are selected to form a three-dimensional array Ix141Ix1256Ix1682. Data on the same columns of Ac are selected to form a matrix of size 500×3. All data are drawn in Figure 11(a). Obviously, Ix1 should be the element in class 3 (black solid ball) and Ix2 should be the element in class 5 (red solid ball). The same results can also be obtained when the data on columns 35, 709, and 929 of Ix1 and Ac are selected. The corresponding results are drawn in Figure 11(b).


## 5. Comparative Experiments

The performance of the proposed algorithm can be reflected through comparative experiments. Therefore, it is necessary to compare it with some usual identification methods, such as support vector machine based schemes [32], particle swarm optimization and support vector machine based schemes [12], artificial neural networks and intelligent filter based schemes [33], PNN and simplified fuzzy adaptive resonance theory map neural networks based schemes [15], and fuzzy c-means based schemes [34]. When the class label of each signal is known, the performance of the above methods can be compared by calculating the classification accuracy and the identification accuracy. The comparison results are shown in Table 2. The samples that cannot be correctly judged are shown in Figure 12.Table 2
Accuracy comparison of different methods (K = 5).
Literature
Gryllias’s [32]
Liu’s [12]
Zarei’s [33]
Jaouher’s [15]
Cannon’s [34]
This work
Classification accuracy (%)
98.20
100.00
96.40
97.80
94.60
100.00
Identification accuracy (%)
75.00
75.00
50.00
75.00
50.00
100.00Figure 12
The samples that cannot be judged correctly by different schemes. (a) The samples that cannot be correctly classified. (b) The samples that cannot be correctly identified.
(a)
(b)Table2 shows that the classification accuracy obtained from Liu’s and our methods is 100%. The value of identification accuracy obtained from the proposed method also reaches 100%. Figure 12(a) shows that some algorithms fail to give the correct classification results for some signals. In terms of signal classification, the number of samples without correct classification for each method is 9, 0, 18, 11, 27, and 0. By comparison, the identification accuracy derived from Zarei’s and Cannon’s scheme is poor, and they fail to identify the sample 3 and 4 (Figure 12(b)). This is mainly because artificial neural networks and fuzzy c-means algorithm can give a judgment result for each signal, whether or not the signal belongs to the determined classes. Although the PNN has some similar problems, they can be skillfully solved by the method of Bivariable correlation analysis. It makes the proposed method have unique advantages in signal identification.


## 6. Concluding Remarks

It is an indisputable fact that the probabilistic neural networks can be used to classify and identify patterns. It has a wide range of application, including the identification of emitter signals. In this paper, a novel classification and identification scheme, which is designed by the WCM and the PNN with correlation analysis, has been proposed for emitter signals. The scheme starts with self-adaptive filtering processing and spectrum analysis, then the WCM and clustering validity indexes including CH, Silh, DB, and Gap are utilized to determine the range of the optimal number of clusters. For different classification numberK, 70 samples nearby each center are selected as training samples to establish PNN classifiers. Finally, the optimal PNN classifier which is used to identify signals is determined by the maximum of the classification validity index D(K). At this stage, the method of Bivariable correlation analysis is cleverly used to improve the identification accuracy of PNN classifier. Experiments show that the proposed method can obtain higher accuracy, and it is more stable than other schemes in identification problems.Finally, it should be pointed out that the scheme presented above is mainly used to classify the obtained signals, which are derived from some pulse emitters. The classification and identification of signals derived from some continuous wave emitters or the mixed emitters of these two types are the next topic to be studied. In addition, the proposed method can also be used to identify other signals, such as, biomedical signals and monitoring signals of digital virtual assets. When the data set to be classified are not numeric, they can first be converted to binary string and then converted to the required format.



---

*Source: 1458962-2018-11-05.xml*