Specific Emitter Identification Based on Synchrosqueezing Transform for Civil Radar

Abstract

Time-frequency (TF) signal features are widely used in specific emitter identification (SEI) which commonly arises in many applications, especially for radar signals. Due to data scale and algorithm complexity, it is difficult to obtain an informative representation for SEI with existing TF features. In this paper, a feature extraction method is proposed based on synchrosqueezing transform (SST). The SST feature has an equivalent dimension to Fourier transform, and retains the most relevant information of the signal, leading to on average approximately 20 percent improvement in SEI for complex frequency modulation signals compared with existing handcrafted features. Numerous results demonstrate that the synchrosqueezing TF feature can offer a better recognition accuracy, especially for the signals with intricate time-varying information. Moreover, a linear relevance propagation algorithm is employed to attest to the SST feature importance from the perspective of deep learning.

1. Introduction

Specific emitter identification (SEI) has always been a challenging problem in signal recognition. With the support of big data and high computing capability, technologies such as deep neural networks (DNN) can perform the identification task efficiently. However, the insufficiency of the training data set for DNN limits the depth of network and hinders performance yielding inaccurate classifications. Furthermore, understanding and interpreting network learning is always a challenging problem. Therefore, an effective feature extraction strategy is needed for SEI tasks. The ideal feature, before the classifier, should have a complete but low dimensional representation of the signal. Neither time or frequency alone can adequately fulfill this requirement when dealing with non-stationary signals. As a powerful non-stationary signal processing tool, TF analysis can reveal the signal local frequency characteristics. However, the mapping from one-dimensional signal to two-dimensional representation can introduce a large dimensionality problem. Attempts to reduce the redundancy in time-frequency (TF) feature have become an important subject [1,2,3]. Towards this end, some works [4,5,6] treat the TF spectrum as images. Then, downsampling of both variables is applied, and the result is sent to the classifier for identification. It is noted, however, that fine details of the emitter signature and modulation profile can be missed by classification confusion between closely resembling emitters.

Two effective methods have been recently introduced to tackle the above problem [7,8], namely the ambiguity function representative slice (AF-RS) and the compressed sensing mask (CS-MASK). These two features not only reveal physical characteristics of measured radar signal, but also alleviate the computational costs. However, when it comes to complex frequency modulation signals, these two methods may perform unstably. For AF-RS, taking a single AF slice near the origin may only capture a part of the modulation information. For CS-MASK, it is complicated, and rather an art, to set the parameter of mask. Both methods face challenges when dealing with complex modulation signals, giving inadequate identification.

To overcome the limitations above, the synchrosqueezing transform (SST) [9], as an empirical mode decomposition-like tool, can be applied since it squeezes the energy from a wide frequency range into a narrow bandwidth near instantaneous frequency points, which can greatly reduce the scale of feature, and retains most modulation information. Hence, in this paper, a feature extraction method based on SST is proposed. This feature shows a good recognition rate, meanwhile, due to flavor and philosophy of feature extraction, it can probably facilitate the interpretation of the deep networks with linear relevance propagation (LRP) algorithm [10,11]. The data from civil aviation meteorological radar and radar signal generators are used to verify the validity of the proposed feature.

The paper is arranged as follows. Two TF feature extraction methods are discussed in the next section. Section 3 illustrates LRP algorithm. The SST method and the corresponding properties are introduced in Section 4; Section 5 proposes the SST feature extraction method. Numerous simulation results based on measured data and the interpretation of deep neural networks are given in Section 6. Finally, we conclude the paper in Section 7.

2. Time-Frequency Feature Extraction Methods

In this section, we briefly introduce the two TF feature extraction methods: Ambiguity Function Representative Slice (AF-RS) and Compressed Sensing Mask (CS-MASK) [7,8]. Ambiguity function (AF) concentrates most signal energy at and near the range-Doppler origin. Because of this property, it is possible to extract a feature with a low dimension by organizing the information near the origin in the ambiguity domain (AD).

2.1. AF-RS

The AF is the inverse Fourier transform of the instantaneous autocorrelation function of a signal over the time variable. In essence, it transforms signals to delay-frequency offset domain, whereas common methods depict signatures in joint time-frequency domain (TFD). Owing to the fact that typical targets do not rapidly transient over a short time period, most of the backscattering energy is not highly offset away from the AD origin. The AF of signal f is given by:

$$A_f(\tau ,\nu )={\int }_{\infty }^{\infty }f(t+\frac{\tau }{2})f^{*}(t-\frac{\tau }{2})e^{-j\nu \tau }dt,$$

(1)

in which $\tau$ is time delay and $\nu$ is frequency offset. It is clear from (1) that when frequency offset is 0, the AF degenerates into the integral of autocorrelation function of the signal with respect to time. According to the definition of autocorrelation function, the maximum of AF is near the origin.

For illustration, consider a single-frequency signal and an LFM signal. The AF of the former is shown in Figure 1. At zero frequency offset, the signal energy reaches the maximum value. The cross-terms arising from the two Euler representations of a real sinusoid are shown parallel to frequency axis. These terms are considered artifacts and do not aid signal identification. On the other hand, the AF of LFM signal is shown in Figure 2. It can be seen that although the energy is mainly concentrated at zero frequency offset, the time-varying frequency of the signal cannot be fully captured in only one slice of zero frequency offset.

Since the energy of AF exhibits a peak at shift frequency being zero, the slices of ambiguity function with shift frequency near zero (including zero-slice) could be considered to be the major representative feature sets of radar emitter pulse which is named representative slice. This kind of feature not only coincides with the characteristics of real radar data but also can reduce the computation complexity.

Figure 3 shows the slices of the LFM signal near the zero-frequency offset, which is from the frequency offset from 0 ($\nu =0$) to 5 ($\nu =5$):

A class-dependent method to classify radar emitter is proposed by Gillespie and Atlas [12,13], which is used to extract initial features of radar signals, followed by a kernel optimization scheme in ambiguity function plane. Paper [8] applied Direct Discriminant Ratio (DRR) as a criterion to rank the kernel, defined as:

$$\begin{array}{c}\hfill DRR(\tau ,\nu )=\sum _{i=1}^C\sum _{j=1}^C\mid {\overline{A}}^i-{\overline{A}}^j{\mid }^2,\end{array}$$

(2)

here, ${\overline{A}}^i$ means the “average” auto-ambiguity function of Class i. Frequency offset $\nu$ is usually set near zero in order to rank the points along the major direction of ambiguity function distribution. By this norm, AF-RS extracts the slices near the origin in delay-frequency offset domain as a feature set. The most representative slice is selected as the feature, resulting in a great reduction in data dimension and an obvious enhancement of identification accuracy. However, there are still some shortcomings, mainly in: (1) Only the slice at a fixed frequency is extracted as a feature, while other slices may contribute to identification in other offsets. (2) Compared with the zero-frequency-offset slice, the representative slice has more feature fusion optimization process, which probably makes it difficult to achieve real-time optimization. It must be based on recognition rate feedback, resulting in an increase in the computational complexity.

2.2. CS-MASK

A compressed sensing mask is based on that ambiguity function and Wigner-Ville distribution is a two-dimensional Fourier transform pair relationship. After transforming the signal into the ambiguity domain, the compressed sensing theory is used to optimize the sparse ambiguity function mask feature.

According to the expression of ambiguity function, the two-dimensional inverse Fourier transform of the Wigner-Ville distribution is the ambiguity function, shown below:

$$\begin{array}{c}\hfill A_f(\tau ,\nu )={\int }_{\infty }^{\infty }f(t+\frac{\tau }{2})f^{*}(t-\frac{\tau }{2})e^{-j\nu \tau }dt={\int }_{\infty }^{\infty }{\int }_{\infty }^{\infty }WV{D}_f(t,f)e^{-j(\nu \tau +\tau f)}dtdf,\end{array}$$

(3)

here $WV{D}_f(t,f)$ is Winger-Ville Distribution (WVD) of signal $f\left(t\right)$, a common quadratic time frequency representation tool, defined as:

$$\begin{array}{c}\hfill WV{D}_f(t,f)={\int }_{\infty }^{\infty }f(t+\frac{\tau }{2})f^{*}(t-\frac{\tau }{2})e^{-j2\pi f\tau }dt.\end{array}$$

(4)

Assuming WVD and ambiguity as length N points sequences, the matrix correlation of them as follows, which are related by two-dimensional Fourier transform:

$$A_{f(N\times 1)}={\mathsf{\Psi }}_{(N\times N)}·W_{D(N\times 1)},$$

(5)

where $\mathsf{\Psi }$ is 2D Fourier transform matrix in size of $N\times N$, and $W_D$ is the matrix form of Winger-Ville Distribution of the processed signal $f\left(t\right)$. To reconstruct Wigner-Ville distribution from ambiguity function with high accuracy and low data dimension, a measurement matrix $\mathsf{\Phi }$ in size of $M\times N$ ($M\ll N$) used in the compressed sensing is defined. Feature $A_{f(M\times 1)}^{Mask}$, an M sparse measurement of AF, is obtained as follows:

$$\begin{array}{c}\hfill A_{f(M\times 1)}^{Mask}={\mathsf{\Psi }}_{(N\times N)}·A_{f(N\times 1)}={\mathsf{\Phi }}_{(M\times N)}·{\mathsf{\Psi }}_{(N\times N)}·W_{D(N\times 1)}.\end{array}$$

(6)

An ideal feature should include as much information as possible with a low dimension. The sparse feature can be obtained according to compressed sensing theory, which can be optimized as:

$$A_{f(M\times 1)}^{Mask}=\underset{A}{arg}min{\parallel A_f\parallel }_1,$$

(7)

$$s.t.\sum _{k=1}^N\sum _{n=1}^M\frac{1}{MN}[{\mathcal{F}}_{2d}^{-1}\left\{A_f\right\}-W_D]\le \epsilon ,$$

(8)

where ${\mathcal{F}}_{2d}^{-1}$ is 2D Fourier inverse transform, $\epsilon$ represents a user-specified bound to constrain the error in an acceptable range. Because the 2D Fourier transform of ambiguity function is WVD, the efficacy of the feature CS-MASK can be examined by the error between its 2D Fourier inverse transform and WVD according to Equation (8).

Figure 4 represents the ambiguity function and CS-MASK feature of an LFM signal. The signal length is 256, the size of the original ambiguity function is $256\times 256$, while the size of CS-MASK extracted according to Equations (7) and (8) is $17\times 17$, meaning the feature dimension has been greatly reduced.

It can be seen from Figure 4c that although CS-MASK only takes the information in a $17\times 17$ region centered near the origin, it contains the most information about energy and frequency change. The amount of information contained in CS-MASK can be measured with its recovered WVD. Figure 5 is a comparison of WVD of the original signal and WVD of the signal reconstructed by CS-MASK.

Nevertheless, there are still two points which can be improved in this method: (1) The selection of the mask size is based on the optimization criterion, which needs to be selected by the reconstructed signal WVD, so the calculation amount is large, although the rapid optimization method in two projects is proposed [7]; (2) The size of the optimized feature is not fixed but will change with certain factors such as signal modulation and noise interference.

3. Linear Relevance Propagation

The LRP method is generally used in the middle layer visualization and feature inversion in image processing with DNN. LRP is based on the layer-wise conservation principle which states that relevance R is conserved in process of back-propagating from one layer to the next, i.e., [11],

$$\begin{array}{c}\hfill R^{\left(1\right)}=\cdots =R^{\left(l\right)}=R^{(l+1)}=\cdots =f\left(x\right).\end{array}$$

(9)

The relevance of the last layer is equivalent to the classification function $f\left(x\right)$ of the input data. In each layer, relevance R is defined as below:

$$\begin{array}{c}\hfill R^{\left(l\right)}=\sum _{n=1}^N{r}^{\left(l\right)},\end{array}$$

(10)

where $r^{\left(l\right)}$ is the relevance of the neurons in each layer. The heatmap $h_n=r_n^1$, the relevance in the first hidden layer, is obtained by the following rule:

$$\begin{array}{c}\hfill r_n^{\left(l\right)}=\sum _m\left(\alpha \frac{z_{n,m}^{-}}{{\sum }_{n^{{}^{\prime }}}{z}_{n^{{}^{\prime }},m}^{-}}+\beta \frac{z_{n,m}^{-}}{{\sum }_{n^{{}^{\prime }}}{z}_{n^{{}^{\prime }},m}^{-}}\right)r_m^{(l+1)},\end{array}$$

(11)

where $\alpha$ and $\beta$ are two fixed parameters during iterations, and $z_{n^{{}^{\prime }},m}^{+}$ and $z_{n^{{}^{\prime }},m}^{-}$ are the positive and negative part of $z_{n^{{}^{\prime }},m}$ severally. $z_{n^{{}^{\prime }},m}$ is the contribution of nth neuron at layer l to the activation of the mth neuron in the next layer $l+1$,

$$\begin{array}{c}\hfill z_{n^{{}^{\prime }},m}=a_n^{\left(l\right)}{\omega }_{n,m}^{(l,l+1)},\end{array}$$

(12)

here $a_n^{\left(l\right)}$ is the activation of neuron n in layer l, and ${\omega }_{n,m}^{(l,l+1)}$ is the weighs from nth neuron in layer l to mth neuro in layer $l+1$.

In this paper, LRP method is used to create a heatmap for the recognition result of each method that reflects which parts of the original data contribute positively to recognition of network and which contributes negatively. Combined with this heatmap, we can understand this network, perhaps even manage to interpret how it works. This research is conducted in the simulation part.

4. Synchrosqueezing Transform

SST, proposed by Daubechies et al., is based on wavelet transform and empirical mode decomposition (EMD) [9]. SST is a post TF processing tool which compresses the spectrum in TFD by squeezing the energy distributed in a frequency range into some certain instantaneous frequency points. SST has been widely applied in many fields of biomedical signals, geographic signals, sonar signals, etc. [14,15].

In this paper, Short Time Fourier Transform (STFT), a commonly used time-frequency analysis tool, is adopted as time frequency transform. STFT converts the signals into a 2D time-frequency domain by performing Fourier Transform in a fixed window traversing the time domain. Considering a signal $f\left(t\right)$, the STFT of $f\left(t\right)$ with a window $g\left(t\right)$ in Schwartz class, is defined as [16,17,18]:

$$\begin{array}{c}\hfill V_f^g(t,f)={\int }_{\mathbb{R}}f\left(\tau \right)g^{*}(\tau -t)e^{-i2\pi f\tau }d\tau ,\end{array}$$

(13)

where $g{\left(t\right)}^{*}$ is the complex conjugate of $g\left(t\right)$, and f is frequency.

Signal $f\left(t\right)$ to be considered in synchrosqueezing techniques in the STFT is defined by:

$$\begin{array}{c}\hfill f\left(t\right)=\sum _{k=1}^K{f}_k\left(t\right),\mathrm{with}{f}_k\left(t\right)=a_k\left(t\right)e^{i{\varphi }_k\left(t\right)},\end{array}$$

(14)

where K is a constant, $a_k>0$ is a continuously differentiable function, ${\varphi }_k$ is a two times continuously differentiable function satisfying ${\varphi }_k^{{}^{\prime }}\left(t\right)>0$ and ${\varphi }_{k+1}^{{}^{\prime }}\left(t\right)>{\varphi }_k^{{}^{\prime }}\left(t\right)$ for all t. In this context, ideal T-F (ITF) representations can be defined as:

$$\begin{array}{c}\hfill IF{T}_q(t,f)=\sum _{k=1}^K{a}_k{\left(t\right)}^q\delta (2\pi f-{\varphi }_k^{{}^{\prime }}\left(t\right)),\end{array}$$

(15)

where q is a positive integer depending on the chosen T-F distribution (TFD), and $\delta$ refers to Dirac Distribution. The STFT being is associated with $q=1$.

Base on this theory, the synchrosqueezing operator reassigns the STFT as follows:

$$\begin{array}{c}\hfill T_f^g(t,f)=\frac{1}{g{\left(0\right)}^{*}}{\int }_{\mathbb{R}}{V}_f^g(t,f)e^{i2\pi ft}\delta (2\pi f-{\widehat{\omega }}_f(t,f))df,g{\left(0\right)}^{*}\ne 0\end{array}$$

(16)

where ${\widehat{\omega }}_f(t,f)={\varphi }^{{}^{\prime }}\left({\widehat{\tau }}_f(t,f)\right)$. ${\widehat{\tau }}_f(t,f)$ is a time reassignment factor, defined as [17]:

$${\widehat{\tau }}_f(t,f)=\frac{1}{T_f^g(t,f)}\int {\int }_{{\mathbb{R}}^{\mathbb{2}}}\tau W_g(\tau -t,\nu -2\pi f)W_f(\tau ,\nu )\frac{d\tau d\nu }{2\pi },T_f^g(t,f)\ne 0,$$

(17)

here $W_g$ and $W_f$ are WVD of $g\left(t\right)$ and $f\left(t\right)$.

By squeezing the frequency components within a certain range of instantaneous frequency, an energy-concentrated TF representation can be obtained. Then, the reconstruction of the original signal is computed by

$$\begin{array}{c}\hfill f\left(t\right)=\frac{1}{C_h}{\int }_{-\infty }^{+\infty }{T}_f^g(t,f)\frac{df}{f},C_h\ne 0,\end{array}$$

(18)

where $C_h=Re\left(C_{\phi }{C}_{\psi }\right),C_{\phi }=e^{-i2\pi ft}{f}^2$ and $C_{\psi }=\frac{1}{2}{\int }_0^{+\infty }\overline{\widehat{\psi }\left(\xi \right)}\frac{1}{\xi }d\xi$. $\overline{\widehat{\psi }\left(\xi \right)}$ is the complex conjugate of the Fourier transform of the mother wavelet $\psi \left(t\right)=\frac{1}{\sqrt{\omega \pi }}{e}^{-\frac{t^2}{2}}{e}^{i2\pi t}$ ($\omega =2\pi f$). More details can be found in [19].

Because of the property that the energy of a signal is mainly concentrated near the instantaneous frequency points, it is only necessary to extract vital frequency, energy, and time-varying information. The effect of first-order SST is better only when the second derivative of the phase with respect to time is approximately zero or much smaller than its first derivative. However, for complex frequency modulation signals, such as secondary frequency hopping signals, the second derivative of phase is not negligible. Hence, vertical synchrosqueezing transform (VSST) is proposed [20]. In this paper, VSST is used as the TF analysis tool because the experimental datasets contain numerous complex modulation signals.

5. The Proposed SST Feature Extraction Algorithm

In this section, the SST feature extraction is proposed to achieve an accurate recognition rate. The property of VSST is elaborated and some experimental results of VSST are shown in Section 5.1. Additionally, the multipath effect, a common phenomenon in radar signal processing, usually leads to an adverse impact on feature representation. VSST can mitigate the multipath effect, which is discussed in Section 5.2 and the SST feature extraction method is presented in Section 5.3

5.1. Vertical Synchrosqueezing Transform

The latest research results [20] have achieved vertical synchrosqueezing transform (VSST) and even high-order form. For nonlinear FM signals, the concentration of the transform is further improved. D.H.Pham gives a VSST expression based on STFT:

$$\begin{array}{c}{T}_f^g(t,f)=\frac{1}{g{\left(0\right)}^{*}}{\int }_{\mathbb{R}}{V}_f^g(t,f)e^{i2\pi ft}\delta (\omega -{\widehat{\omega }}_f(t,f))df,g{\left(0\right)}^{*}\ne 0,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\tilde{\omega }}_f(t,f)& =\left\{\begin{array}{ccc}\hfill & {\widehat{\omega }}_f(t,f)+{\widehat{q}}_f(t,f)(t-{\widehat{\tau }}_f(t,f))\hfill & ,{\partial }_t\widehat{\tau }(t,f)\ne 0\hfill \\ \hfill & {\widehat{\omega }}_f(t,f)\hfill & ,\mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

Here, ${\widehat{\omega }}_f(t,f)$ and ${\widehat{\tau }}_f(t,f)$ are TF rearrangement factors, ${\tilde{\omega }}_f(t,f)$ is the second-order estimate of the instantaneous frequency of signal. Compared with the first-order SST, VSST obtains more accurate instantaneous frequency estimation by adding the adjustment factor ${\widehat{q}}_f(t,f)$, meaning the TF concentration can be further improved. The original signal can be quickly recovered using the following formula [20]:

$$\begin{array}{c}\hfill f\left(t\right)={\int }_{\mathbb{R}}{T}_f^g(t,f)df.\end{array}$$

(20)

Figure 6 shows the VSST and the inverse transform of synthetic linear frequency modulation pulse signal. The signal information is shown below:

$$\begin{array}{c}\hfill sig\left(t\right)=cos(2\pi 30t^2+2\pi 20t).\end{array}$$

(21)

The normalized fs = 256, the total length of the signal is 512 and the signal to noise ratio (SNR) is 0 dB.

From Figure 6, VSST retains the time-varying information of the original signal and has a lossless inverse transformation, which brings great convenience for feature extraction and identification. Figure 7 is an enlarged partial view of VSST. It can be seen that most energy is concentrated near instantaneous frequency points. Compared to STFT and WVD, SST is a more accurate TF representation because of its perfect TF concentration.

5.2. Suppression on Multipath Effect by VSST

The multipath effect refers to a phenomenon that after electromagnetic wave propagates through different paths, each component arrives at the receiver in different time, which causes great interference due to the superposition of each phase [20]. The mulipath effect denotes as:

$$\begin{array}{c}\hfill f_m\left(t\right)=f\left(t\right)+E_m\ast f(t-\tau ),\end{array}$$

(22)

where $f_m$ is the signal with multipath effect, $E_m$ is the ratio of amplitude of multipath effect to that of signal. Figure 8 shows the SST of the simulation signal $sig\left(t\right)$ in Equation (21) with multipath effect under different $E_m$. Here, time delay has two conditions: $\tau =16$ and 128 sampling points. The SNR is 0 dB. From Figure 8a,d,g,j, it can be seen that multipath effect can cause a great adverse impact on SST feature. Figure 8c,f,i,l manifest that VSST can separate the multipath component from the signal and suppress it effectively.

5.3. Feature Extraction

It is clear from the various experimental results of SST that energy from the signal in TFD is mainly concentrated within a limited range, while some redundant information is distributed around the rest area. Based on this property, SST feature will be derived as follows: First, STFT together with VSST are performed on one signal to obtain the TF spectrum. Second, the energy sequence $E\left(n\right)$ and the phase trajectory sequence $f\left(n\right)$ are extracted from the SST domain. It is worth noting that in order to embody more detailed information, the compression bandwidth of the energy sequence can be appropriately extended. Finally, the two sequences are merged into one new feature matrix and then sent to the classifier for identification. The feature matrix is shown below:

$$\begin{array}{c}\hfill \mathbf{M}=\left[\begin{array}{cccc}f\left(1\right)& f\left(2\right)& \dots & f\left(n\right)\\ E\left(1\right)& E\left(2\right)& \dots & E\left(n\right)\end{array}\right]\end{array}$$

(23)

Take the SST feature of the real measured signal from DatasetII as an example (The description of DatasetII is given in the next section), see in Figure 9, the length of the signal is $N=500$, and the size of original spectrum matrix is $250\times 500$. In comparison, the size of the feature matrix extracted by SST is reduced to $2\times 500$.

The process of the proposed algorithm (including feature extraction and LRP interpretation) is given as follows and Figure 10:

Step1: Signal Preprocessing.

Preprocess the signals to obtain invalid and neat signals including demean, aligning and normalization. (Denoise is not considered in this paper since it may remove important details by accident, while SST can keep the main part of the signal within the narrow band and remove the rest of the signal with noise.) The obtained data usually contain abnormal value, noise, and interference which influence the identification in this algorithm. In addition, normalization is required in order to avoid the adverse impact of different signal values for identification.

Step2: Feature Extraction.

In this step, firstly it is required to perform STFT on the processed signals. Then SST is used to suppress the energy of STFT. In SST spectrum, instantaneous frequency and suppressed energy are extracted as a new TF feature matrix.

Step3: Identification.

Finally, the extracted SST feature in Step2 is delivered into DNN for training and identification. It should be noted that the structure of DNN for each type of data is identical. A better recognition rate means that the feature may include more useful information for training the network.

Step4: Interpretation.

Combined with the LRP algorithm, heatmaps of various signals can be obtained. Based on these heatmaps, we can manage to understand how the networks work, and even try to give a reasonable interpretation of its identification mechanism.

6. Experimental Results

In this section, the supervised classifier is applied to identify the four sets of measured radar data, and compare SST feature with other three features to verify the effectiveness of the proposed method. DatasetI contains 10 types of signals from 10 identical types of civil aviation meteorological radars. It includes 500 approximate single-frequency signal samples; the next three data sets are all collected from radar signal generators. DatasetII contains 10 types of 500 single frequency signal samples; DatasetIII contains 10 types of 500 LFM signal samples; DatasetIV contains 10 types of 500 LFM signal samples. The reason for choosing these four data sets is that they are very typical and representative real measured signals of civil radars. We randomly select a signal from each data set to perform STFT and SST in order to show the TF characteristics of each data set, shown in Figure 11.

Table 1. Processing time for each dataset.

Dataset	DatasetI	DatasetII	DatasetIII	DatasetIV
Processing time/ms	$64.60$	$65.42$	$65.86$	$76.12$

shows the processing time to obtain SST feature of each data set.

The extracted SST features of each data set are then sent to the classifier for training and identification. The training rate is set from $10\%$ to $70\%$, and deep neural network (DNN) with fully connected layers is used as a classifier. The recognition rates under different features are shown in

Table 2. Recognition rate for DatasetI.

Training Rate	10	20	30	40	50	60	70
PSE	$\mathbf{93}.\mathbf{10}$	$\mathbf{92}.\mathbf{18}$	$90.61$	$\mathbf{96}.\mathbf{43}$	$91.78$	$92.80$	$94.01$
AF-RS	$88.69$	$91.27$	$89.96$	$91.06$	$88.59$	$89.58$	$87.70$
CS-MASK	$92.01$	$90.75$	$92.87$	$90.57$	$90.03$	$93.17$	$90.16$
SST	$85.68$	$91.45$	$\mathbf{91}.\mathbf{56}$	$94.03$	$\mathbf{92}.\mathbf{71}$	$\mathbf{94}.\mathbf{45}$	$\mathbf{96}.\mathbf{02}$

,

Table 3. Recognition rate for DatasetII.

Training Rate	10	20	30	40	50	60	70
PSE	$56.60$	$73.87$	$86.67$	$\mathbf{93}.\mathbf{62}$	$93.01$	$86.69$	$\mathbf{95}.\mathbf{34}$
AF-RS	$62.00$	$63.25$	$69.73$	$71.33$	$84.67$	$80.02$	$86.00$
CS-MASK	$\mathbf{93}.\mathbf{54}$	$\mathbf{91}.\mathbf{26}$	$\mathbf{95}.\mathbf{98}$	$93.54$	$92.51$	$\mathbf{94}.\mathbf{19}$	$91.92$
SST	$67.18$	$62.40$	$78.74$	$89.45$	$81.33$	$84.59$	$84.42$

,

Table 4. Recognition rate for DatasetIII.

Training Rate	10	20	30	40	50	60	70
PSE	$41.56$	$62.50$	$82.54$	$90.67$	$91.60$	$93.50$	$92.67$
AF-RS	$52.39$	$64.75$	$74.86$	$84.33$	$81.60$	$83.00$	$87.33$
CS-MASK	$52.56$	$72.75$	$82.57$	$83.33$	$83.20$	$89.00$	$86.00$
SST	$\mathbf{88}.\mathbf{67}$	$\mathbf{93}.\mathbf{75}$	$\mathbf{95}.\mathbf{14}$	$\mathbf{98}.\mathbf{00}$	$\mathbf{99}.\mathbf{20}$	$\mathbf{98}.\mathbf{23}$	$\mathbf{98}.\mathbf{60}$

and

Table 5. Recognition rate for DatasetIV.

Training Rate	10	20	30	40	50	60	70
PSE	$80.55$	$82.90$	$89.15$	$92.95$	$96.45$	$93.05$	$96.15$
AF-RS	$70.65$	$71.25$	$72.95$	$73.55$	$76.05$	$86.95$	$84.05$
CS-MASK	$69.10$	$72.75$	$73.65$	$75.95$	$75.30$	$81.00$	$75.30$
SST	$\mathbf{89}.\mathbf{30}$	$\mathbf{97}.\mathbf{75}$	$\mathbf{98}.\mathbf{70}$	$\mathbf{98}.\mathbf{30}$	$\mathbf{98}.\mathbf{00}$	$\mathbf{98}.\mathbf{45}$	$\mathbf{99}.\mathbf{35}$

and Figure 12. It can be seen that for DatasetI, SST feature has only an inapparent advantage over the power spectrum (PSE), AF-RS and CS-MASK [7,8]. For the single-frequency signal in DatasetII, the effect of SST features is even slightly worse than PSE and CS-MASK. It is because the form of single-frequency signals is simple, and the energy is highly concentrated at the origin of AD domain, hence the information extracted by the mask is very complete. From the perspective of PSE, single-frequency signals are time-invariant signals, so one-dimensional energy information can represent signals enough, while SST feature may introduce redundant information which reduces the recognition rate.

In contrast, for DatasetIII and IV, the SST feature obtains the best recognition rate among these four features. This is because the modulation form of LFM signal is more complicated than single-frequency signal, and theoretically it needs more accurate information to identify. SST feature contains more detailed time-varying information of signal than other three features, so it can obtain a better representation of LFM signal.

To further verify the effectiveness of the SST feature, the LRP algorithm [10,11] is used to obtain the TF heatmaps in reverse. The parameters $\alpha$ and $\beta$ in Equation (11) are set as $\alpha =3$, and $\beta =-1.5$. The number of the hidden layers in auto-coder is 3, and $relu$ is selected as the activation function of each layer. Figure 13 shows heatmaps of deep network for LFM signals from DatasetIII and DatasetIV. The brightness of pixels represents the contribution to identification of deep network. It is apparent from Figure 13 that the brightest region is concentrated near the SST feature, meaning SST feature is akin to the part of spectrum which contributes most to the identification of deep network. In other words, the decisive information in deep network is not disarray of some pixels in spectrum, but some special parts of the spectrum make the most contribution. It can be observed that there exists a horizontal stripe and two vertical stripes in Figure 13a,b. The horizontal stripe may be caused by harmonic components in the signal, and the vertical stripes may be caused by phase hopping. In Figure 13c,d, a dim horizontal stripe appears within the first quarter of the entire time, which is probably due to the noise at the beginning from signal generators. The specific causes of this phenomenon need to be further studied. These experimental results demonstrate the validity of the SST feature, and also reflect that recognition of the deep network is possibly interpretable.

7. Conclusions

In this paper, we proposed a time-frequency feature extraction method based on SST. The energy sequence and phase trajectory sequence are merged into a feature matrix, which retains almost complete time-frequency information, and holds clear physical meaning with a small dimension. The proposed feature can obtain better recognition rates especially for complex modulation signals due to its high energy concentration and precise instantaneous frequency which can contain more beneficial information, meanwhile, suppress redundant or even negative information, such as noise and multipath effect, to classifier. The superiority of the method in comparison to PSE, AF-RS, and CS-MASK [7,8] is verified by the identification of the measured radar data. Furthermore, combined with the proposed feature and LRP algorithm [11], it is possible to understand and interpret the inner mechanism of the identification by deep networks.