Multiple-Antenna Emitters Identification Based on a Memoryless Power Amplifier Model

Abstract

Power amplifier (PA) nonlinearity is typically unique at the radio frequency (RF) front-end for particular emitters. It can play a crucial role in the application of specific emitter identification (SEI). In this paper, under the Multi-Input Multi-Output (MIMO) multipath communication scenario, two data-aided approaches are proposed to identify multi-antenna emitters using PA nonlinearity. Built upon a memoryless polynomial model, the first approach formulates a linear least square (LLS) problem and presents the closed-form solution of nonlinear coefficients in a MIMO system by means of singular value decomposition (SVD) operation. Another alternative approach estimates nonlinear coefficients of each individual PA through nonlinear least square (NLS) solved by the regularized Gauss–Newton iterative scheme. Moreover, there are some practical discussions of our proposed approaches about the mismatch of the order of PA model and the rank-deficient condition. Finally, the average misclassification rate is derived based on the minimum error probability (MEP) criterion, and the proposed approaches are validated to be effective through extensively numerical simulations.

1. Introduction

Specific emitter identification (SEI) is committed to distinguish individual radiation sources by using essential radio frequency fingerprint (RFF) features extracted from different emitters. It can be applied to military communication [1], physical layer authentication [2], and enhancement of the security in wireless network, such as very high frequency (VHF) radio networks, Wi-Fi networks, cognitive radios, cellular networks [3], and so on.

In general, based on different states of signals, the identifiable RFF features for SEI are usually extracted from the transient signal or the steady-state signal. The transient signal, actually the turn-on signal, carries unique and unintentional information that is advantageous to emitter identification, and the features underlying are mainly extracted from the instantaneous amplitude, phase, frequency, and energy envelope [4,5,6]. Nevertheless, it is difficult to capture the transient signal since the duration time is often too short to use. As for the steady-state signal, many researchers focus on extracting the distinguishable features that are generated by hardware imperfection of the components inside the radiation source through advanced signal processing techniques. On one hand, the statistical characteristics of the original RF signal, such as the high order spectrum, have been used as the features to identify different emitters [7,8]. On the other hand, the Time-Frequency Analysis (TFA) methods [9,10], Wavelet Transform (WT) [11,12], and Hilbert–Huang Transform (HHT) [13,14] are successively applied to extract the transform domain characteristics from the received RF signals. However, these methods have little knowledge of impairments inside the individual emitter, and the performance can be easily affected by the wireless channels. Hence, there are many additional works in literature to model the characteristic of hardware imperfection of the internal component such as digital to analog converter (DAC) [15], modulator [16,17], and power amplifier (PA) [18,19,20,21], etc., and to extract the unique RF front-end feature for a particular emitter. This paper concentrates on the extraction of PA nonlinear features to identify communication emitters with multiple antennas under Multiple-Input Multiple-Output (MIMO) multipath channels.

It is well known that the MIMO transmission scheme can improve the spectral efficiency by introducing additional spatial diversity. In practice, transceivers with multiple transmit and receive antennas have shown their powerful merits and served as the major mechanism for current and future wireless communication systems. Meanwhile, due to the increasing number of antennas for MIMO emitters, there will be more diverse PA nonlinearities available at the RF front end for SEI. A lot of research efforts have been afforded to analyze the nonlinearity incurred by PAs with the Saleh model [22], polynomial model [23], and Volterra model [24], etc. in MIMO systems, whereas they are mainly devoted to implementing predistortion or PA linearization. However, to the authors’ best knowledge, there are few open results of SEI in MIMO communication systems. In [21], Li uses a modified artificial bee colony (ABC) algorithm to estimate the coefficients of the Hammerstein model, a simplified version of the Volterra model, in a MIMO system. However, the impact of wireless propagation channels is not considered, and the ABC algorithm also appears complicated to obtain the optimal solution in the context of SEI. Recently, in [25,26], an estimation of signal parameters via rotational invariance technique (ESPRIT)-based approach, which takes advantage of the multiple antennas at the receiver to separate the RFF from wireless channel, is proposed for RFF estimation in orthogonal frequency division multiplexing (OFDM) systems, whereas it is only suitable for a Single-Input Multiple-Output (SIMO) system rather than the MIMO one. In [19,20], a memoryless polynomial model is used to characterize the nonlinearity of PA, and a data-aided iterative algorithm is proposed to estimate nonlinear coefficients of the PA model for SEI from the observations in both MIMO and single-input single-output (SISO) scenarios.

Treating the fact that all PAs of a multiple-antenna emitter are independent from each other and following a memoryless polynomial model, in this paper, we propose two data-aided solutions that are different from the idea of [19,20] in the MIMO multipath scenario. Given received signals, we extend the method in [27] to the MIMO multipath scenario; a closed-form solution of the nonlinear coefficients is thus obtained by combining the linear least square (LLS) and singular value decomposition (SVD) methods. An alternative but more effective approach is also presented through solving a nonlinear least square (NLS) problem with independent variables consisting of both channel coefficients and nonlinear coefficients of the PA model. Furthermore, we explicitly provide deep discussion on the parameter estimate bias for the general case of unknown order of the PA model. In particular, it is proved in this paper that the rank-deficient property of both the NLS and LLS problems are in association with the amplitude level of the training sequence. The average misclassification rate based on the minimum error probability (MEP) criterion is theoretically derived, and we finally verify the proposed approaches via extensive numerical simulations.

The rest of the paper is organized as follows. In Section 2, we introduce the memoryless nonlinear model of PA in the MIMO multipath system. In Section 3, we present the linear and nonlinear frameworks for SEI, respectively. Then, the practical discussions for the proposed approaches and the misclassification rate are separately given in Section 4 and Section 5. In Section 6, numerical results are presented to demonstrate the effectiveness. Section 7 summarizes the paper.

2. Preliminaries and Problem Formulation

This paper mainly considers the scenario that K communication emitters equipped with multiple antennas are identified through a multiple-antenna receiver in the MIMO multipath environment. Since constant modulus modulation schemes such as phase-shift keying (PSK) may generally introduce less nonlinear distortion, we thus consider how to extract the underlying PA nonlinearity of each emitter from the received signal, with the help of a quadrature amplitude modulation (QAM) training sequence. In other words, PA nonlinearity is treated as a unique RFF of the emitter to fulfill the identification task, and we assume QAM is used by emitters when communicating with the receiver.

2.1. Memoryless Nonlinear PA Modeling

Generally, the memoryless polynomial model is simple in expression and can describe the intermodulation distortion of an RF PA well, such as PA9440 amplifier [28], in the narrowband communication system. More specifically, the polynomial coefficients are corresponding to the intermodulation coefficient such as third-order interception point ($I{P}_3$) and fifth-order interception point ($I{P}_5$). In this work, we choose the memoryless polynomial model to characterize the nonlinear behavior of all PA units in a multi-antenna emitter; then, the relationship between baseband equivalent input and output of the PA [20] at jth antenna branch can be written as

$$x_j\left(n\right)=\sum _{p=1}^{(P+1)/2}{\alpha }_{2p-1,j}·{\left(s_j\left(n\right)\right)}^p·{\left(s_j^{*}\left(n\right)\right)}^{p-1},$$

(1)

where P denotes the max order of the PA model and can be configured as $P=3,5,7,\dots$. At the RF front-end of emitters, the frequency components that resulted from the even terms in the model can be removed by the bandpass filter, thus the even terms are ignored in Equation (1). $s_j\left(n\right)$ is the baseband equivalent input of the nonlinear system, denoting the nth QAM symbol transmitted over the jth antenna. The superscript “*” is the conjugate operator. ${\alpha }_{2p-1,j}$ denotes the normalized $(2p-1)$th order PA coefficient of the jth antenna and without loss of generality, we set ${\alpha }_{1,j}=1$ hereafter. $x_j\left(n\right)$ stands for the response of the nonlinear system.

In this paper, we suppose that each emitter is equipped with J antennas. Therefore, the normalized discrete-time baseband equivalent form of the nonlinear distortion model for the multi-antenna system can be expressed as

$$\left[\begin{array}{c}{x}_1\left(n\right)\hfill \\ x_2\left(n\right)\hfill \\ ⋮\hfill \\ x_J\left(n\right)\hfill \end{array}\right]=\left[\begin{array}{c}{\overline{s}}_1\left(n\right)\hfill \\ {\overline{s}}_2\left(n\right)\hfill \\ ⋮\hfill \\ {\overline{s}}_J\left(n\right)\hfill \end{array}\right]+\sum _{p=2}^{(P+1)/2}\left[\begin{array}{cccc}{\alpha }_{2p-1,1}& 0& \cdots & 0\\ 0& {\alpha }_{2p-1,2}& \dots & ⋮\\ ⋮& ⋮& \ddots & 0\\ 0& \cdots & 0& {\alpha }_{2p-1,J}\end{array}\right]\left[\begin{array}{c}{\left({\overline{s}}_1\left(n\right)\right)}^p·{\left({\overline{s}}_1^{*}\left(n\right)\right)}^{p-1}\hfill \\ {\left({\overline{s}}_2\left(n\right)\right)}^p·{\left({\overline{s}}_2^{*}\left(n\right)\right)}^{p-1}\hfill \\ ⋮\hfill \\ {\left({\overline{s}}_J\left(n\right)\right)}^p·{\left({\overline{s}}_J^{*}\left(n\right)\right)}^{p-1}\hfill \end{array}\right],$$

(2)

where ${\overline{s}}_j\left(n\right)$ is the normalized version of $s_j\left(n\right)$.

2.2. MIMO Multipath Channel

The propagation channel considered in this paper is a linear discrete MIMO system with J transmit antennas and R receive antennas. Given the length of the training sequence N, the signal received at the antenna r can be represented by

$${\mathbf{y}}_r={\mathbf{h}}_{r1}\otimes {\mathbf{x}}_1+{\mathbf{h}}_{r2}\otimes {\mathbf{x}}_2+\cdots +{\mathbf{h}}_{rJ}\otimes {\mathbf{x}}_J+{\mathbf{v}}_r,$$

(3)

where ⊗ denotes the convolution operation. ${\mathbf{x}}_j$ is the nonlinear distortion signal transmitted over the jth antenna, ${\mathbf{h}}_{rj}$ represents the channel impulse response between the jth transmit antenna to the rth receive antenna and remains time-invariant during the data receiving process, and ${\mathbf{v}}_r$ is the zero-mean additive white Gaussian noise received at the antenna r. Thus, we can unfold Equation (3) naturally in the form of

$${\mathbf{y}}_r={\mathbf{X}}^{\left(1\right)}{\mathbf{h}}_{r1}+{\mathbf{X}}^{\left(2\right)}{\mathbf{h}}_{r2}+\dots +{\mathbf{X}}^{\left(J\right)}{\mathbf{h}}_{rJ}+{\mathbf{v}}_r,$$

(4)

where the order of the channel is L, ${\mathbf{X}}^{\left(j\right)}\in {\mathbb{C}}^{(N+L-1)\times L}$ is a Toeplitz matrix populated by ${\mathbf{x}}_j$ in Equation (3). Note that the nonlinear coefficients are cross-coupled with channel coefficients in Equation (4); therefore, the major target arising in this paper is to get the accurate estimations of the separate nonlinear coefficient $\left\{{\alpha }_{2p-1,j}\right\}$ via the received signal $\left\{{\mathbf{y}}_r\right\}$ with the aid of training sequences in MIMO multipath scenarios, before we can identify a specific emitter.

3. The Proposed Estimation Approaches

As mentioned above, it has been reported in [19,20] that the nonlinear coefficients of PAs can be estimated with two stages. The first stage establishes the initial estimation of channel coefficients and nonlinear coefficients through some well-designed training sequences sorted by amplitude. Then, an iterative method is applied at the second stage to eliminate the estimate bias. In [27], a PA parameter estimator combined the best linear unbiased estimation (BLUE) and singular value decomposition (SVD) is proposed for the SISO system. Since the method in [27] can get a closed-form solution of the PA nonlinear coefficients and has no constraint on the ordering of the amplitude of training symbols, which is more practical compared to the one in [19], we extend it to the MIMO multi-path scenarios and mark it as linear method in MIMO (LMM). Furthermore, we propose an alternative method to extract PA parameters in a nonlinear least square (NLS) manner. It should be noted that both LMM and NLS approaches can be degraded into SISO systems.

3.1. The LMM Approach

Note that, in [27], if the product terms of the channel coefficients and the nonlinear coefficients are substituted by some new integrated variables, we can obtain a set of linear equations with regard to the new variables, thus the product terms can be readily solved by LLS. Then, the only thing left to us is to extract ${\alpha }_{2p-1,j}$ from the solution.

According to Equation (4), we first vectorize signals received by all R antennas into vector ${\mathbf{y}}_{vec}$, i.e., ${\mathbf{y}}_{vec}={\left[{\mathbf{y}}_1^T,{\mathbf{y}}_2^T,\dots ,{\mathbf{y}}_R^T\right]}^T\in {\mathbb{C}}^{R(N+L-1)}$, where the superscript “T” is the transpose operator. In addition, the signal at the receiver side now is

$${\mathbf{y}}_{vec}={\mathbf{D}}_s{\mathbf{h}}_{\alpha }+{\mathbf{v}}_{vec},$$

(5)

where ${\mathbf{v}}_{vec}$ is the corresponding reshaped noise vector. ${\mathbf{h}}_{\alpha }\in {\mathbb{C}}^{RJL(P+1)/2}$ is the integrated vector composed of all independent variables and can be represented as:

$${\mathbf{h}}_{\alpha }={\left[{{\mathbf{u}}_{1,0}}^T,\dots ,{{\mathbf{u}}_{1,L-1}}^T\left|\dots \left|{{\mathbf{u}}_{R,0}}^T,\dots ,{{\mathbf{u}}_{R,L-1}}^T\right.\right.\right]}^T,$$

(6)

with

$${\mathbf{u}}_{r,l_h}={\mathbf{W}}_{r,l_h}{\mathit{\alpha }}_{vec},$$

(7)

$${\mathit{\alpha }}_j={\left[1,{\alpha }_{3,j},{\alpha }_{5,j},\dots ,{\alpha }_{P,j}\right]}^T,$$

(8)

$${\mathit{\alpha }}_{vec}={\left[{{\mathit{\alpha }}_1}^T,{{\mathit{\alpha }}_2}^T,\dots ,{{\mathit{\alpha }}_J}^T\right]}^T,$$

(9)

$${\mathbf{W}}_{r,l_h}=blkdiag\left(h_{r1}\left(l_h\right){\mathbf{I}}_{\alpha },h_{r2}\left(l_h\right){\mathbf{I}}_{\alpha },\dots ,h_{rJ}\left(l_h\right){\mathbf{I}}_{\alpha }\right),$$

(10)

in which $l_h=0,\dots ,L-1$, ${\mathbf{I}}_{\alpha }\in {\mathbb{C}}^{(P+1)/2\times (P+1)/2}$ is a unit matrix and $blkdiag(·)$ is the block diagonalization function. In addition, ${\mathbf{D}}_s$ is defined by

$${\mathbf{D}}_s=blkdiag\left(\underset{R}{\underbrace{{\mathbf{d}}_s,{\mathbf{d}}_s,\dots ,{\mathbf{d}}_s}}\right),$$

(11)

with ${\mathbf{d}}_s\in {\mathbb{C}}^{(N+L-1)\times JL(P+1)/2}$ being constructed by the known training sequence. To elaborate the process, we provide a numerical example in Appendix A.

Consequently, the least square (LS) estimation of ${\mathbf{h}}_{\alpha }$ in Equation (5) is

$${\widehat{\mathbf{h}}}_{\alpha }={\mathbf{D}}_s^{†}{\mathbf{y}}_{vec},$$

(12)

where the superscript “†” denotes the pseudo-inverse operation. It is worth noting that, compared to the BLUE method in [27], the least square solution of Equation (12) has no requirement on the estimation of noise power. Obviously, the condition of a unique solution to Equation (12) is natural that ${\mathbf{D}}_s$ is full column rank and $N\ge JL(P+1)/2-L+1$ is satisfied.

Afterwards, we define

$$\mathbf{Q}=\left[{\widehat{\mathbf{u}}}_{1,0},\dots ,{\widehat{\mathbf{u}}}_{1,L-1}\left|\dots \left|{\widehat{\mathbf{u}}}_{R,0},\dots ,{\widehat{\mathbf{u}}}_{R,L-1}\right.\right.\right]$$

(13)

and further perform SVD on the matrix $\mathbf{Q}\in {\mathbb{C}}^{J(P+1)/2\times RL}$ to get a closed-form estimation of the PA parameters in the MIMO system. Since the normalized first-order nonlinear coefficients are assumed equal to 1, the execution steps of PA nonlinear coefficients estimator can be summarized as:

(1)

Reshape the observations according to Equation (5), and estimate ${\mathbf{h}}_{\alpha }$ according to Equation (12).

(2)

Reshape the ${\widehat{\mathbf{h}}}_{\alpha }$ into the matrix $\mathbf{Q}$, then perform SVD on ${\mathbf{Q}}_j={\mathbf{U}}_j{\sum }_j{\mathbf{V}}_j^H$, where ${\mathbf{Q}}_j$ is a submatrix consisting of the $\left((j-1)(P+1)/2+1\right)$th to $\left(j(P+1)/2\right)$th rows of the matrix $\mathbf{Q}$.

(3)

Estimate the nonlinear coefficients of jth transmitting antenna as follows:

$${\widehat{\mathit{\alpha }}}_j=\frac{1}{{\mathbf{U}}_j^{(1,1)}}{\mathbf{U}}_j^{(:,1)},$$

(14)

where ${\mathbf{U}}_j^{(:,1)}$ and ${\mathbf{U}}_j^{(1,1)}$ are the first column and first element of ${\mathbf{U}}_j$, respectively.

3.2. The NLS Approach

Bearing in mind the received signal in Equation (4) for MIMO multipath transmission system, we can now further expand the expression with a series of nonlinear equations due to the existence of product term of the nonlinear coefficient ${\alpha }_{2p-1,j}$ and the channel coefficient $h_{rj}\left(l_h\right)$. As a consequence, the problem of nonlinear coefficients estimation can thus be alternatively transformed into a NLS optimization one when introducing a training sequence with length N.

In order to get a more robust solution, we choose to solve a constrained NLS optimization problem with $q=RJL+J\left(P-1\right)/2$ independent variables, and the cost function can be given by:

$$\underset{\mathbf{z}}{min}{G}_{\gamma }\left(\mathbf{z}\right)=\underset{\mathbf{z}}{min}\frac{1}{2}{∥\mathbf{g}\left(\mathbf{z}\right)∥}_2^2+\frac{1}{2}{∥\gamma \mathbf{z}∥}_2^2,$$

(15)

where ${∥·∥}_2$ denotes 2-norm. $\mathbf{z}\in {\mathbb{C}}^q$ is the vector consisting of the independent variables and can be described as $\mathbf{z}={\left[{\mathbf{z}}_h^T,{\mathbf{z}}_{\alpha }^T\right]}^T$, in which ${\mathbf{z}}_{\alpha }$ and ${\mathbf{z}}_h$ are respectively as

$${\mathbf{z}}_{\alpha }={\left[{\alpha }_{3,1},\dots ,{\alpha }_{P,1},\dots ,{\alpha }_{3,J},\dots ,{\alpha }_{P,J}\right]}^T$$

(16)

$${\mathbf{z}}_h={\left[{{\mathbf{H}}_1}^T,{{\mathbf{H}}_2}^T,\dots ,{{\mathbf{H}}_R}^T\right]}^T,$$

(17)

with

$${\mathbf{H}}_r={\left[h_{r1}\left(0\right),\dots ,h_{rJ}\left(0\right),\dots ,h_{r1}(L-1),\dots ,h_{rJ}(L-1)\right]}^T.$$

(18)

${∥\gamma \mathbf{z}∥}_2^2/\phantom{{∥\gamma z∥}_2^{2}2}2$ is the regularization. $\mathbf{g}\left(\mathbf{z}\right)\in {\mathbb{C}}^{R(N+L-1)}$ denotes the residual function: ${\mathbb{C}}^q\to {\mathbb{C}}^{R(N+L-1)}$ with $R(N+L-1)\ge q$, that is,

$$\mathbf{g}\left(\mathbf{z}\right)={\left[\Delta {{\mathbf{y}}_1}^T,\Delta {{\mathbf{y}}_2}^T,\dots ,\Delta {{\mathbf{y}}_R}^T\right]}^T,$$

(19)

with

$$\Delta {\mathbf{y}}_r={\mathbf{y}}_r-\left({\mathbf{X}}^{\left(1\right)}{\mathbf{h}}_{r1}+{\mathbf{X}}^{\left(2\right)}{\mathbf{h}}_{r2}+\dots +{\mathbf{X}}^{\left(J\right)}{\mathbf{h}}_{rJ}\right).$$

(20)

According to [29], the regularized Gauss–Newton iterative method is introduced here to figure out $\mathbf{z}$, i.e.,

$${\mathbf{z}}_{i+1}={\mathbf{z}}_i+\Delta {\mathbf{z}}_i$$

(21)

with

$$\Delta {\mathbf{z}}_i=-{\left({\left(\mathbf{J}\left({\mathbf{z}}_i\right)\right)}^H\mathbf{J}\left({\mathbf{z}}_i\right)+{\gamma }^2{\mathbf{I}}_q\right)}^{-1}{\left[\begin{array}{c}\mathbf{J}\left({\mathbf{z}}_i\right)\\ \gamma {\mathbf{I}}_q\end{array}\right]}^H\left[\begin{array}{c}\mathbf{g}\left({\mathbf{z}}_i\right)\\ \gamma {\mathbf{z}}_i\end{array}\right].$$

(22)

The superscript “H” denotes conjugate transposition, and ${\mathbf{I}}_q$ is a $q\times q$ unit matrix. $\mathbf{J}\left(\mathbf{z}\right)$ is the Jacobian matrix of $\mathbf{g}\left(\mathbf{z}\right)$, that is, the first-order derivative of $\mathbf{g}\left(\mathbf{z}\right)$ with respect to $\mathbf{z}$.

In general, given the Jacobian matrix $\mathbf{J}\left(\mathbf{z}\right)$ with full column rank, one can apply the regularized Gauss–Newton method in Equation (21) to optimize problem Equation (15) and can eventually obtain the optimal solution ${\mathbf{z}}_{opt}$. In addition, the appropriate regularization factor $\gamma$ can improve the condition number of the inverse matrix in Equation (22), which guarantees the robustness of the NLS optimization.

4. Practical Discussions on the Proposed Approaches

As mentioned in Section 3, we can get a closed-form solution for the problem of nonlinear coefficient estimation through the LMM approach. Alternatively, an iterative NLS approach can also be applicable to extract PA parameters from observations. However, in practice, there are two main factors that affect the accuracy of our proposed algorithms. The first factor is the case where the order of PA model is mismatched between the transmitting and receiving ends. The second one is the case where the matrix ${\mathbf{D}}_s$ and the Jacobian matrix $\mathbf{J}\left(\mathbf{z}\right)$ are rank-deficient. In the next Section 4.1 and Section 4.2, we first theoretically analyze the impact of the mismatched model. After that, we present the rank-deficient condition of the ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ matrix in Section 4.3.

4.1. Overdetermined Order of the PA Model

In this subsection, we consider the impact of the overdetermined order of the PA model on the estimation accuracy of the nonlinear coefficients. Assume that all the PAs of an emitter actually have the same model order as P. If we obtain an overdetermined order of the PA model beforehand, e.g., $P_0$ with $(P_0>P)$, the expected observation of ${\mathbf{y}}_{vec}={\mathbf{D}}_s^{(1,\dots ,P)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P)}+{\mathbf{v}}_{vec}$ may be formulated as ${\mathbf{y}}_{vec}={\mathbf{D}}_s^{(1,\dots ,P_0)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}+{\mathbf{v}}_{vec}$ instead. Since ${\mathbf{D}}_s^{(1,\dots ,P)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P)}$ equals $\left({\mathbf{D}}_s^{(1,\dots ,P_0)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}-{\mathbf{D}}_s^{(P+2,\dots ,P_0)}{\mathbf{h}}_{\alpha }^{(P+2,\dots ,P_0)}\right)$, the estimation of ${\mathbf{h}}_{\alpha }$ can thus be expressed as follows:

$$\begin{array}{cc}{\widehat{\mathbf{h}}}_{\alpha }\hfill & ={\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}{\mathbf{y}}_{vec}\hfill \\ & ={\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}\left({\mathbf{D}}_s^{(1,\dots ,P)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P)}+{\mathbf{v}}_{vec}\right)\hfill \\ & ={\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}-{\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}\left({\mathbf{D}}_s^{(P+2,\dots ,P_0)}{\mathbf{h}}_{\alpha }^{(P+2,\dots ,P_0)}\right)+{\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}{\mathbf{v}}_{vec},\hfill \end{array}$$

(23)

where ${\mathbf{D}}_s^{(a,\dots ,b)}$ denotes a matrix ${\mathbf{D}}_s$ that is populated by the ath to the bth order versions of training sequence, and ${\mathbf{h}}_{\alpha }^{(a,\dots ,b)}$ denotes a vector ${\mathbf{h}}_{\alpha }$ that is populated by the corresponding ath to the bth order nonlinear coefficients. Some numerical examples for the construction process are also provided in Appendix A.

Note that the $(P+2)$th to the $P_0$th order the nonlinear coefficients and ${\mathbf{h}}_{\alpha }^{(P+2,\dots ,P_0)}$ should be all zeros in this case, then, we have ${\widehat{\mathbf{h}}}_{\alpha }={\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}+{\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}{\mathbf{v}}_{vec}$, which is equivalent to the unbiased estimation of the ${\mathbf{h}}_{\alpha }^{(1,\dots ,P)}$, and the additive white noise has no effect on the unbiasedness of the estimation results. Therefore, we can conclude that, if the order of the PA model is overdetermined, the estimation of the nonlinear coefficients obtained by the LMM approach is still unbiased.

4.2. Underdetermined Order of the PA Model

In the sequel, we attempt to analyze the case that the order of the PA model is underdetermined, i.e., $P_0$ is lower than the actual P. It is easy to understand that ${\mathbf{D}}_s^{(1,\dots ,P)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P)}$ now equals $\left({\mathbf{D}}_s^{(1,\dots ,P_0)}{\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}+{\mathbf{D}}_s^{(P_0+2,\dots ,P)}{\mathbf{h}}_{\alpha }^{(P_0+2,\dots ,P)}\right)$, so that the estimation of ${\mathbf{h}}_{\alpha }$ can be written as:

$$\begin{array}{c}{\widehat{\mathbf{h}}}_{\alpha }={\mathbf{h}}_{\alpha }^{(1,\dots ,P_0)}+{\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}\left({\mathbf{D}}_s^{(P_0+2,\dots ,P)}{\mathbf{h}}_{\alpha }^{(P_0+2,\dots ,P)}\right)+{\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}{\mathbf{v}}_{vec}.\hfill \end{array}$$

(24)

Attention must be paid to non-zero ${\mathbf{h}}_{\alpha }^{(P_0+2,\dots ,P)}$ in Equation (24) due to the non-zero nonlinear PA coefficients through the $(P_0+2)$th to the Pth order. The resulting ${\widehat{\mathbf{h}}}_{\alpha }$ is a biased estimation w.r.t. ${\mathbf{h}}_{\alpha }^{(1,\dots ,P)}$ and the estimate bias equals ${\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}\left({\mathbf{D}}_s^{(P_0+2,\dots ,P)}{\mathbf{h}}_{\alpha }^{(P_0+2,\dots ,P)}\right)$.

Furthermore, if we define a guess as:

Guess 1.

If the training squences are discrete amplitude communication symbols such as PSK and QAM, then

$${\left({\mathbf{D}}_s^{(1,\dots ,P_0)}\right)}^{†}{\mathbf{D}}_s^{(P_0+2,\dots ,P)}=blkdiag\left(\underset{RJL}{\underbrace{\mathbf{c},\mathbf{c},\dots ,\mathbf{c}}}\right),$$

(25)

where$\mathbf{c}\in {\mathbb{C}}^{m\times n}$is a constant matrix with$m=\frac{P_0+1}{2}$and$n=\frac{P-P_0}{2}$.

Then, the estimate bias of the nonlinear PA coefficients can be given as the following Proposition 1.

Proposition 1.

If Guess 1 holds true, the estimate bias of the nonlinear PA coefficients obtained by SVD according to Equation (14) is shown as

$$\left\{\begin{array}{}\mathrm{(26a)}& \Delta {\alpha }_{3,j}=\frac{(c_{21}{\alpha }_{P_0+2,j}+\dots +c_{2n}{\alpha }_{P,j})-(c_{11}{\alpha }_{P_0+2,j}+\dots +c_{1n}{\alpha }_{P,j}){\alpha }_{3,j}}{1+(c_{11}{\alpha }_{P_0+2,j}+\dots +c_{1n}{\alpha }_{P,j})}\hfill \mathrm{(26b)}& ⋮\hfill \mathrm{(26c)}& \Delta {\alpha }_{P_0,j}=\frac{(c_{m1}{\alpha }_{P_0+2,j}+\dots +c_{mn}{\alpha }_{P,j})-(c_{m1}{\alpha }_{P_0+2,j}+\dots +c_{mn}{\alpha }_{P,j}){\alpha }_{P_0,j}}{1+(c_{11}{\alpha }_{P_0+2,j}+\dots +c_{1n}{\alpha }_{P,j})}.\hfill \end{array}\right.$$

Proof of Proposition 1.

See Appendix B. □

Therefore, we can conclude that, when the order of the PA model is underdetermined, the bias term of each nonlinear coefficient is related to training sequence and the higher-order nonlinear coefficients. In addition, the high-order nonlinear coefficients (i.e., ${\alpha }_{P_0+2,j}$ to ${\alpha }_{P,j}$) are still unknown.

4.3. Rank Deficiency Condition of the Proposed Approaches

As mentioned earlier, it is impossible to determine a unique solution for the LLS problem in Equation (12) when the matrix ${\mathbf{D}}_s$ is rank-deficient. In addition, for the NLS problem in Equation (15), it is also hard for the regularized Gauss–Newton method in Equation (21) to find the global optimal solution when the Jacobian matrix $\mathbf{J}\left(\mathbf{z}\right)$ is rank-deficient at every $\mathbf{z}$. However, since QAM is employed by all emitters, the training sequence has only a limited level of amplitude. We find that the rank attributes of ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ are both associated with the maximum order of the PA model. The specific relationship can be revealed by Proposition 2 as follows.

Proposition 2.

Given QAM symbols with M modulus values as the training sequence and the maximum P order of the PA model, the matrix ${\mathbf{D}}_s$ is rank-deficient and the Jacobian matrix $\mathbf{J}\left(\mathbf{z}\right)$ is also rank-deficient at every $\mathbf{z}$, if $M<(P+1)/2$.

Proof of Proposition 2.

See Appendix C. □

According to the Proposition 2, the ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ are full rank as long as the amplitude type of signal is sufficient. Therefore, our proposed approach can be applied in some single-carrier communication systems with higher order QAM modulation such as 16-QAM, 64-QAM and 256-QAM, which have 3, 9, and 32 different modulus values, respectively. Actually, it appears that our proposed approaches are readily suitable for other popular wireless communication systems, such as MIMO-OFDM systems, where there will be no rank-deficient problem of ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ in nature since the amplitude of the transmitted signal is generally continuous. Moreover, it is not necessary to estimate the channel order in a MIMO-OFDM system due to its ability to resist multipath effects, which will make the proposed approach more practical.

5. Error Rate Analysis For Classification

In this paper, we apply a minimum error probability (MEP) criterion [30] based on Bayesian theory to classify different emitters, and the RFF feature of each emitter is composed of the estimated nonlinear coefficients and can be expressed as follows:

$$\mathbf{a}={\left[{{\mathbf{a}}_3}^T,{{\mathbf{a}}_5}^T,\dots ,{{\mathbf{a}}_P}^T\right]}^T,$$

(27)

with

$${\mathbf{a}}_p={[Re\left({\widehat{\alpha }}_{p,1}\right),\dots ,Re\left({\widehat{\alpha }}_{p,J}\right),Im\left({\widehat{\alpha }}_{p,1}\right),\dots ,Im\left({\widehat{\alpha }}_{p,J}\right)]}^T,$$

(28)

where $p=3,5,\dots ,P$.

Here, we take the case where there are two emitters as an example to give the derivation of the decision criteria, and the binary hypothesis test model can be considered as

$$\left\{\begin{array}{}\mathrm{(29a)}& C_1:\mathbf{a}={\mathbf{m}}_1+{\mathbf{v}}_{\alpha },\hfill \mathrm{(29b)}& C_2:\mathbf{a}={\mathbf{m}}_2+{\mathbf{v}}_{\alpha },\hfill \end{array}\right.$$

where $C_i$ denotes the category i, $i=1,2$; ${\mathbf{v}}_{\alpha }$ is composed of the residual additive Gaussian noise in the estimations of the nonlinear coefficients and we assume that each element of ${\mathbf{v}}_{\alpha }$ obeys a Gaussian distribution with a zero mean and a variance of ${\delta }^2$; ${\mathbf{m}}_1\in {\mathbb{C}}^{J(P-1)}$ and ${\mathbf{m}}_2\in {\mathbb{C}}^{J(P-1)}$ are respectively the mean vectors of the estimated feature vector $\mathbf{a}$ for the two emitters, which can be obtained from the samples collected offline. Thus, the decision rule can be derived based on MEP criterion as follows:

$$\left\{\begin{array}{}\mathrm{(30a)}& {({\mathbf{m}}_2-{\mathbf{m}}_1)}^T\mathbf{a}>\frac{{{\mathbf{m}}_2}^T{\mathbf{m}}_2-{{\mathbf{m}}_1}^T{\mathbf{m}}_1}{2},\mathbf{a}\in C_2,\hfill \mathrm{(30b)}& {({\mathbf{m}}_2-{\mathbf{m}}_1)}^T\mathbf{a}<\frac{{{\mathbf{m}}_2}^T{\mathbf{m}}_2-{{\mathbf{m}}_1}^T{\mathbf{m}}_1}{2},\mathbf{a}\in C_1,\hfill \end{array}\right.$$

where the test statistic is $a_{ts}={({\mathbf{m}}_2-{\mathbf{m}}_1)}^T\mathbf{a}$ and the decision threshold is $thr=({{\mathbf{m}}_2}^T{\mathbf{m}}_2-{{\mathbf{m}}_1}^T{\mathbf{m}}_1)/\phantom{({{\mathbf{m}}_2}^T{\mathbf{m}}_2-{{\mathbf{m}}_1}^T{\mathbf{m}}_1)2}2$.

For simplicity, if we assume that the variables in $\mathbf{a}$ are independent of each other, then the test statistic $a_{ts}\sim \mathcal{N}\left(m_{ts},{\delta }_{ts}^2\right)$, for the i-th emitter, the mean is as $m_{ts}^{\left(i\right)}={({\mathbf{m}}_2-{\mathbf{m}}_1)}^T{\mathbf{m}}_i,i=1,2$, and the variance is as ${\delta }_{ts}^2={\displaystyle \sum _{k=1}^{J(P-1)}}{\delta }^2{\beta }_k^2$ with ${\beta }_k$ being the k-th element of $({\mathbf{m}}_2-{\mathbf{m}}_1)$. As a result, with the assumption of equally probable hypotheses, the average misclassification rate based on MEP criterion can be derived as

$$\begin{array}{cc}{\overline{P}}_e\hfill & =\frac{1}{2}\left[P_r\left(\left.a_{ts}>thr\right|C_1\right)+P_r\left(\left.a_{ts}<thr\right|C_2\right)\right]\hfill \\ & =\Phi \left(-\frac{{({\mathbf{m}}_2-{\mathbf{m}}_1)}^T({\mathbf{m}}_2-{\mathbf{m}}_1)}{2{\delta }_{ts}}\right)\hfill \\ & =\Phi \left(-\frac{\sqrt{{({\mathbf{m}}_2-{\mathbf{m}}_1)}^T({\mathbf{m}}_2-{\mathbf{m}}_1)}}{2\delta }\right),\hfill \end{array}$$

(31)

where $\Phi \left(x\right)={\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}dx$ and $\Phi (-x)=1-\Phi \left(x\right)$. When x increases, $\Phi (-x)$ decreases.

As for the case where the estimations of the nonlinear coefficients are unbiased, the mean values of the RFF feature for the two emitters are readily as

$${\mathbf{m}}_i={[{{{\mathbf{m}}_i}^{\left(3\right)}}^T,{{{\mathbf{m}}_i}^{\left(5\right)}}^T,\dots ,{{{\mathbf{m}}_i}^{\left(P\right)}}^T]}^T$$

(32)

with

$${{\mathbf{m}}_i}^{\left(p\right)}={[Re\left({\alpha }_{p,1}^{\left(i\right)}\right),\dots ,Re\left({\alpha }_{p,J}^{\left(i\right)}\right),Im\left({\alpha }_{p,1}^{\left(i\right)}\right),\dots ,Im\left({\alpha }_{p,J}^{\left(i\right)}\right)]}^T,$$

(33)

where ${\alpha }_{p,j}^{\left(i\right)}$ is the nonlinear coefficient of the i-th emitter. However, for the case of biased estimations, ${\tilde{\mathbf{m}}}_i={\mathbf{m}}_i+\Delta {\mathbf{m}}_i$, where

$$\begin{array}{c}\Delta {\mathbf{m}}_i={[{\Delta {{\mathbf{m}}_i}^{\left(3\right)}}^T,{\Delta {{\mathbf{m}}_i}^{\left(5\right)}}^T,\dots ,{\Delta {{\mathbf{m}}_i}^{\left(P\right)}}^T]}^T\hfill \end{array}$$

(34)

with

$$\begin{array}{c}\Delta {{\mathbf{m}}_i}^{\left(p\right)}={[Re\left([\Delta {\alpha }_{p,1}^{\left(i\right)},\dots ,\Delta {\alpha }_{p,J}^{\left(i\right)}]\right),Im\left([\Delta {\alpha }_{p,1}^{\left(i\right)},\dots ,\Delta {\alpha }_{p,J}^{\left(i\right)}]\right)]}^T\hfill \end{array}$$

(34)

being the mean of the biases with $i=1,2$ in the estimation of RF fingerprint features. Therefore, according to Equation (31), the average misclassification rate for this case can be obtained by:

$$\begin{array}{cc}{\overline{P}}_e\hfill & =\Phi \left(-\frac{{({\tilde{\mathbf{m}}}_2-{\tilde{\mathbf{m}}}_1)}^T({\tilde{\mathbf{m}}}_2-{\tilde{\mathbf{m}}}_1)}{2\sqrt{{\displaystyle \sum _{k=1}^{J(P-1)}}{\delta }^2{\tilde{\beta }}_k^2}}\right)\hfill \\ & =\Phi \left(-\frac{\sqrt{{\mathit{\theta }}_{\mathbf{1}}+{\mathit{\theta }}_{\mathbf{2}}+{\mathit{\theta }}_{\mathbf{3}}}}{2\delta }\right),\hfill \end{array}$$

(36)

where $\Delta {\beta }_k$ is the k-th element of $(\Delta {\mathbf{m}}_2-\Delta {\mathbf{m}}_1)$ and

$$\{\begin{array}{}\mathrm{(37a)}& {\mathbf{\theta }}_{\mathbf{1}}={({\mathbf{m}}_2-{\mathbf{m}}_1)}^T({\mathbf{m}}_2-{\mathbf{m}}_1)\hfill \mathrm{(37b)}& {\mathbf{\theta }}_{\mathbf{2}}=2{({\mathbf{m}}_2-{\mathbf{m}}_1)}^T(\Delta {\mathbf{m}}_2-\Delta {\mathbf{m}}_1)\hfill \mathrm{(37c)}& {\mathbf{\theta }}_{\mathbf{3}}={(\Delta {\mathbf{m}}_2-\Delta {\mathbf{m}}_1)}^T(\Delta {\mathbf{m}}_2-\Delta {\mathbf{m}}_1).\hfill \end{array}$$

Therefore, we can conclude as follows:

(1)

The ${\overline{P}}_e$ decreases as the variance of the additive noise decrease.

(2)

The ${\overline{P}}_e$ decreases as the difference between the mean values of the RF fingerprint feature for the two emitters increases, which indicates that the PA parameters of each emitter should be designed to be as different as possible in a bid to achieve better performance.

(3)

According to Equation (31), as for a fixed P, more nonlinear coefficients are used as features will make the $\sqrt{{({\mathbf{m}}_2-{\mathbf{m}}_1)}^T({\mathbf{m}}_2-{\mathbf{m}}_1)}$ larger, which obviously leads to better classification performance.

(4)

According to Equation (36), when there are biases in the estimations, the additional terms ${\mathit{\theta }}_{\mathbf{2}}$ and ${\mathit{\theta }}_{\mathbf{3}}$ may cause ${\overline{P}}_e$ to decrease compared to the case where there is no bias.

More generally, as for K ($K>2$) emitters, with the assumption of equally probable hypotheses, the average misclassification rate based on MEP criterion can be represented as

$${\overline{P}}_e=1-{\overline{P}}_c=1-\frac{1}{K}\left[P_r\left(D_1\left|C_1\right.\right)+P_r\left(D_2\left|C_2\right.\right)+\cdots +P_r\left(D_K\left|C_K\right.\right)\right],$$

(38)

where ${\overline{P}}_c$ is the average correct classification rate; $D_k(k=1,2,\dots ,K)$ indicates the discriminant domain of the kth class and it shows as Equation (39)

$$\left\{\begin{array}{c}{\left({\mathbf{m}}_k-{\mathbf{m}}_1\right)}^T\mathbf{a}>\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_1^T{\mathbf{m}}_1\right)/\phantom{\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_1^T{\mathbf{m}}_1\right)2}2\hfill \\ ⋮\hfill \\ {\left({\mathbf{m}}_k-{\mathbf{m}}_{k-1}\right)}^T\mathbf{a}>\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_{k-1}^T{\mathbf{m}}_{k-1}\right)/\phantom{\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_{k-1}^T{\mathbf{m}}_{k-1}\right)2}2\hfill \\ {\left({\mathbf{m}}_k-{\mathbf{m}}_{k+1}\right)}^T\mathbf{a}>\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_{k+1}^T{\mathbf{m}}_{k+1}\right)/\phantom{\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_{k+1}^T{\mathbf{m}}_{k+1}\right)2}2\hfill \\ ⋮\hfill \\ {\left({\mathbf{m}}_k-{\mathbf{m}}_K\right)}^T\mathbf{a}>\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_K^T{\mathbf{m}}_K\right)/\phantom{\left({\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_K^T{\mathbf{m}}_K\right)2}2.\hfill \end{array}\right.$$

(39)

Therefore, the $P_r\left(D_k\left|C_k\right.\right)$ can be obtained by solving the multiple integrals in Equation (40)

$$\begin{array}{c}{P}_r\left(D_k\left|C_k\right.\right)\hfill \\ ={\int }_{\frac{{\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_K^T{\mathbf{m}}_K}{2}}^{+\infty }\cdots {\int }_{\frac{{\mathbf{m}}_k^T{\mathbf{m}}_k-{\mathbf{m}}_1^T{\mathbf{m}}_1}{2}}^{+\infty }f\left(w_1\right)\dots f\left(w_{k-1}\right)f\left(w_{k+1}\right)\dots f\left(w_K\right){\mathbf{d}}_{w_1}\dots {\mathbf{d}}_{w_{k-1}}{\mathbf{d}}_{w_{k+1}}\dots {\mathbf{d}}_{w_K},\hfill \end{array}$$

(40)

where $f\left(w_i\right)=\frac{1}{\sqrt{2\pi }{\delta }_{w_i}}exp(-\frac{{(w_i-{({\mathbf{m}}_k-{\mathbf{m}}_i)}^T{\mathbf{m}}_k)}^2}{2{\delta }_{w_i}^2})$ with $w_i={({\mathbf{m}}_k-{\mathbf{m}}_i)}^T\mathbf{a}$ ($i=1,2,\dots ,K$ and $i\ne k$), and ${\delta }_{w_i}$ is a standard deviation similar to the ${\delta }_{ts}$ in Equation (31).

6. Numerical Results

6.1. Simulation Setting

In the following simulations, if not specified, we use a $2\times 2$ MIMO channel with Rayleigh multipath fading model, and the number of paths L is set to 2. The training sequences are QAM symbols with $N=35$ for proper complexity and performance trade-off. The regularization factor $\gamma$ is set to 0.05 for the trade-off between the variance and bias in the estimation of NLS approach. The complex nonlinear coefficients lead to AM–AM and AM–PM distortion, where the AM–AM distortion is mainly caused by the $I{P}_3$ and $I{P}_5$, and higher order intermodulation distortion is usually ignored [28,31]. Generally, the actual communication system has specific requirements for the out-of-band spectral emission level of RF signals. Therefore, for the rationality of PA parameters in reality, we select the PA parameters that AM–AM characteristic obeys the method in [28,31] and the AM–PM is obtained by slightly adjusting the phase of the parameter in [20]. The nonlinear PA coefficients are displayed in

Table 1. Nonlinear coefficients of the separate PA for each RF chain.

	${\mathit{\alpha }}_{3,1}$	${\mathit{\alpha }}_{3,2}$	${\mathit{\alpha }}_{5,1}$	${\mathit{\alpha }}_{5,2}$
Emitter1	$0.0210-0.0365i$	$0.0290-0.0503i$	$-0.0099-0.0011i$	$-0.0158-0.0018i$
Emitter2	$0.0182+0.0028i$	$0.0131+0.0020i$	$0.0014-0.0008i$	$0.0009-0.0005i$
Emitter3	$0.0780-0.0488i$	$0.0566-0.0354i$	$-0.0111-0.0384i$	$-0.0070-0.0242i$
Emitter4	$0.0270-0.0247i$	$0.0196-0.0179i$	$0.0001-0.0063i$	$0.0000-0.0040i$
Emitter5	$0.0241+0.0393i$	$0.0174+0.0285i$	$0.0023-0.0097i$	$0.0014-0.0061i$
Emitter6	$0.0514-0.0520i$	$0.0372-0.0377i$	$0.0112-0.0226i$	$0.0070-0.0143i$
Emitter7	$0.0129-0.0167i$	$0.0102-0.0132i$	$-0.0006-0.0024i$	$-0.0004-0.0015i$
Emitter8	$0.0246+0.0226i$	$0.0195+0.0179i$	$0.0011-0.0062i$	$0.0007-0.0039i$
Emitter9	$0.0335+0.0410i$	$0.0266+0.0326i$	$0.0081-0.0137i$	$0.0051-0.0086i$
Emitter10	$0.0256-0.0069i$	$0.0204-0.0055i$	$0.0022-0.0033i$	$0.0014-0.0021i$

, and it has been verified that the out-of-band emission level of the amplified signals is about 40 dBr, which meets the requirement of general agreement.

In all experiments, we use the Normalized Mean Squared Error (NMSE) to evaluate the estimation accuracy of the nonlinear coefficients, i.e.,

$$er{r}_{\alpha }=10log10\left(E\left\{\frac{{∥{\mathit{\alpha }}_{est}-{\mathit{\alpha }}_{true}∥}_2^2}{{∥{\mathit{\alpha }}_{true}∥}_2^2}\right\}\right),$$

(41)

where ${\mathit{\alpha }}_{true}={\left[{\alpha }_{3,1},{\alpha }_{3,2},{\alpha }_{5,1},{\alpha }_{5,2},\dots ,{\alpha }_{P,1},{\alpha }_{P,2}\right]}^T$, ${\mathit{\alpha }}_{est}$ is the estimation of the ${\mathit{\alpha }}_{true}$, and $E\{·\}$ denotes the Mathematical Expectation.

6.2. Simulation Results

6.2.1. The Impact of SNR

Note that Liu’s algorithm in [19] assumes all PAs in one emitter are the same. To facilitate the comparison of the proposed approaches with Liu’s algorithm, we extend Liu’s algorithm to be suitable for our MIMO scenario. For brevity, we use the name Modified Liu Algorithm (MLA) to represent the modified version in the following. Here, we set P to 5, and the training sequence is 16-QAM symbols. At first, we give comparisons among the MLA, LMM, and NLS in Figure 1 and Figure 2, where the third-order and fifth-order coefficients are combined as a classification feature, from which one can see that the performance of the LMM and NLS are apparently better than that of MLA, especially at low signal to noise ratio (SNR) regime. One possible reason lies in that when P equals 5 and the training sequence is composed of 16-QAM symbols, the ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ are both full column rank, whereas the matrix populated by the training sequence for the initial estimation of MLA method is rank-deficient regardless of how long the training sequence is. In addition, we notice that the NLS approach also performs better than the LMM one, which is because the introduction of regularization during the process of optimization. Furthermore, in order to explore the impact of iterations on performance of NLS approach, when SNR is 25 dB, we give some results in Figure 3 and Figure 4 to demonstrate the convergence speed of NLS approach, where 19 iterations are needed to obtain the optimal solution.

In Figure 5, we fit a PA with $P=5$, when third-order and fifth-order coefficients are combined as a classification feature, the classification performance is better than when only the third-order coefficient is used. This confirms the analysis that, for a fixed P, features with more nonlinear coefficients can effectively improve classification performance in Section 5. In the next simulations, if not specified, we use the estimated third-order and fifth-order coefficients together as classification features.

Moreover, in Figure 6 and Figure 7, we present results of the proposed approaches under a $4\times 4$ MIMO scenario; as you can see, more RF channels bring more nonlinear coefficients, which is beneficial to improve classification performance. Here, for convenience, we take the nonlinear coefficients of Emitter2 and Emitter4 in Table 1 as those of the first 4-antenna emitter, and take the nonlinear coefficients of Emitter3 and Emitter6 in Table 1 as those of the second 4-antenna emitter.

6.2.2. Identification for Multiple Emitters

In this simulation, we test the performance of proposed methods under the multi-user case. Here, we set the maximum order of PAs in individual emitters to 5, and 16-QAM symbols are used as the training sequence. The detailed nonlinear coefficients of each emitter are displayed in Table 1. The results are shown in Figure 8 and Figure 9, where the SNR is set to 25 dB, and each method is performed by 1000 Monte Carlo simulations.

As we can see from the Figure 8 and Figure 9, our proposed approaches are also robust and applicable to the multi-user case, and average misclassification rates of them are both increased as the number of emitters increased. Obviously, the greater the difference in nonlinear characteristics between emitters, the higher the resolution of the proposed approaches. In this simulation, as far as the nonlinear coefficients in Table 1 are concerned, the proposed approaches can resolve up to seven emitters at most.

6.2.3. The Confirmation of Practical Discussions

Hereafter, we set the number of emitters to 2 (i.e., emitter1 and emitter2), the maximum order of the PA model P to 5. The first to fifth order parameters are selected from Table 1. When we use a PA model with $P_0=3$ to fit the nonlinearity of the emitters, the results of Figure 10 confirm the conclusion in the Section 4.2 that the underdetermined order of the PA model indeed leads to the biases in the estimations of PA parameters. However, the results in Figure 11 show that the classification performance can be improved by using the biased estimations as features, which reveals that the biases in Equation (36) can make the difference between the classification features of emitters larger. Where the training sequence is 16-QAM symbols, the estimated third-order are used as classification features, and the legends “Well-estimated LMM(NLS)” in Figure 10 and Figure 11 denote that results obtained from a well-determined order of PA model with $P_0=P=3$.

Then, according to the Proposition 2, when we use a PA model with $P_0=P=5$ to fit the nonlinearity of emitters, if 4-QAM, which is a constant modulus modulation, is adopted by each emitter, then it is clear that the ${\mathbf{D}}_s$ and $\mathbf{J}\left(\mathbf{z}\right)$ are both rank-deficient, and there is no unique solution in this case. Here, we compare the estimation accuracy of when taking 4-QAM and 16-QAM symbols as training sequence; the results are shown in Figure 12, which obviously corroborates the correctness of Proposition 2.

6.3. Experimental Results

In this subsection, we design a preliminary verification experiment to explore the effectiveness of our proposed approaches in reality. Due to the limitation of single-channel acquisition, we limit the proposed approaches in the SISO-OFDM system and validate them in a 802.11.g-based wireless local area network (WLAN) that is very common in real life. Therefore, we build an experiment platform, which is shown in Figure 13, to collect the measured router data. In this platform, we first use a LeCroy WaveMaster 813Zi-A Oscilloscope (Chestnut Ridge, NY, USA) equipped with a single antenna to acquire the RF signals from three TL-WR740 routers communicating with a smart phone on the 2.412 GHz channel, respectively. In this experiment, we use the extended NLS approach as an example to estimate nonlinear coefficients of each router, here we use the “wlan toolbox” in MATLAB R2019a (MathWorks, Natick, Massachusetts, USA) to perform pre-processing such as timing, synchronization, and de-frequency offset on the acquired RF signals. Since the nonlinearity of PA mainly caused by $I{P}_3$, we set $P=3$. Finally, a MEP-based classifier is used to identify the individual router. Note that the memoryless polynomial model may not be able to describe the nonlinearity of the PA in a broadband WLAN system, whereas

Table 2. The mean and variance of the estimated nonlinear coefficients in Figure 14 for the three routers.

	Router1	Router2	Router3
Mean	$0.0522-0.0013i$	$0.0579+0.0008i$	$0.0441-0.0000i$
Variance	$7.1501e^{-6}$	$2.6937e^{-6}$	$3.9458e^{-6}$

indicates that, according to the estimated PA coefficients in Figure 14, the mean is far greater than the variance for each router, thus the three routers are identifiable based on the extended NLS approach. Moreover, the average misclassification rates of the three routers are all 0. Therefore, we can conclude that the mismath of PA model does not affect classification performance. In addition, we also compare the power spectral density among measured data and simulated data for three routers in Figure 15, Figure 16 and Figure 17, respectively, where the legend “Measured baseband OFDM symbols” denotes downconverted acquired RF signals, the legend “Amplified simulated baseband OFDM symbols” denotes simulated baseband OFDM signals amplified by the PA with measured nonlinear coefficients, and the legend “Raw simulated baseband OFDM symbols” denotes raw simulated baseband OFDM signals, all of their power being normalized. In order to explain the results in Figure 15, Figure 16 and Figure 17, we calculate the NMSE between the PSDs of “Amplified simulated baseband OFDM symbols” and that of “Measured baseband OFDM symbols” for each router, and the NMSEs of three routers are, respectively, $-7.6472$, $-7.5301$ and $-7.3589$ dB, which reveal that the PSD of the signal reconstructed by using the estimated PA coefficient can well fit that of the measured signal.

7. Conclusions

This paper investigates the SEI scheme for multiple-antenna communication emitters, using PA nonlinearity as RFF features with the assumption that all PAs of a multiple-antenna emitter are independent from each other. Both the LMM and the NLS approaches are proposed to estimate the nonlinear coefficients in association with the memoryless polynomial PA model, where a closed-form solution can be obtained by the LMM approach, and the alternative NLS approach achieves better performance by adopting a regularized Newton–Gauss scheme. Practical discussion on the PA model mismatch is presented, and some theoretical results about the estimate bias and rank-deficient condition are provided to guide the design and implementation of the SEI over MIMO channels. In addition, an error rate analysis is also introduced for the MEP classifier. Simulation results demonstrate that the proposed approaches outperform the other existing schemes, especially in the rank-deficient case, and are effective to deal with SEI in MIMO communication systems. Moreover, the proposed approaches are verified to be effective on a 802.11.g-based experiment platform.