Blind Modulation Identification of Underwater Acoustic MPSK Using Sparse Bayesian Learning and Expectation Maximization

Abstract

This paper presents a likelihood-based algorithm for identifying different phase shift keying (PSK) modulations, i.e., BPSK, QPSK, and 8PSK. This algorithm selects the modulation type that maximizes a loglikelihood function that is based on the known original constellation associated with the constellation of the received signals for the candidate modulation types. However, there are two problems in non-cooperative underwater acoustic Multiple Phase Shift Keying (MPSK) modulation identification based on the likelihood method. One is the original constellation, which as prior information is unknown. The other is the underwater acoustic multipath channel makes the constellation distort seriously. In this paper, we solved these problems by combining sparse bayesian learning (SBL) with expectation maximization (EM). The specific steps are as follows. Firstly, blind channel equalization can be achieved by channel impulse response (CIR), which is estimated by sparse bayesian learning in single input multi output (SIMO) underwater acoustic channel. Subsequently, we used expectation maximization to compensate amplitude attenuation and phase offset, as the original constellation of MPSK is unknown. Finally, modulation can be successfully identified by the Quasi Hybrid Likelihood Ratio Test (QHLRT). The simulation results show that the channel estimation method based on SBL can eliminate the influence of channel effectively, and the EM algorithm can make the received constellation converge to the preset constellation in the case of unknown original transmit constellation, which effectively solves these two problems. We use the proposed SBL-EM-QHLRT method to achieve an identification rate of more than 95% in underwater acoustic multipath channels with Signal to Noise Ratio (SNR) higher than 15 dB, which provides a new idea for modulation identification of non-cooperative underwater acoustic MPSK.

1. Introduction

This paper deals with the blind modulation identification of signals transmitting through underwater acoustic multipath channel. We assume that the received signal is Multiple Phase Shift Keying (MPSK), i.e., BPSK, QPSK, and 8PSK, but do not assume any prior knowledge of the modulation order.

Underwater acoustic MPSK is widely used in information transmission; it is mainly used in high-speed communication scenarios [1,2]. MPSK can also be used as a part of Direct Sequence Spread Spectrum (DSSS) to achieve long range data transmission [3,4], and has been abundantly used for underwater Orthogonal Frequency Division Multiplexing (OFDM) subcarrier modulation [5,6,7]. With the development of underwater acoustic communication technology and the increasing demand in the field of underwater acoustic signal detection, modulation identification of non-cooperative MPSK has become a research hotspot. MPSK modulation identification is an important basis for parameter estimation, scrambling, and blind demodulation [8,9,10], particularly in military applications. Various authors have discussed the MPSK modulation identification in the literature and have proposed many schemes; however, most of the application scenarios are focused on the radio. These methods can be divided into two different groups: likelihood function-based and cyclic cumulant-based methods.

Regarding the MPSK modulation identification, various methods have been proposed based on higher-order cumulants and cyclic cumulants [11,12,13,14,15,16]. These methods assumed that the signals were in baseband.But in practical application, it is necessary to estimate the parameters of communication system in order to achieve baseband processing. Because the estimated parameters are susceptible to the channel, the higher-order cumulants and cyclic cumulants based method is not effective in practical application.

A different class of approaches to PSK modulation identification that employed the likelihood-based method. In [17,18], the method indicated that the system needed high SNR in order to obtain accurate identification results and this method is only suitable for noise channel. A similar likelihood-based method for identifying MPSK signals was described in [19]. However, this method needs to know the parameters of MPSK modulation. In [20], the authors introduced a likelihoodbased method for BPSK, QPSK, and 8PSK in multipath channels. This method does not take into account that the original constellation is unknown and the frequency-selective channel was slow-varying. In [21], a new MPSK modulation identification method based on likelihood function is proposed. This method identify modulation without knowledge of the symbol rate, carrier frequency and noise power in the received signals. However, this method is suitable for noise channel, its performance is poor in multipath channels. In non-cooperative environment, it is difficult to estimate the channel impulse response (CIR) of multipath channel. Therefore, the likelihood-based modulation identification method is not suitable for multipath channels generally. In [19], the authors concluded that the likelihood-based method performs well in modulation identification when the CIR and noise power are known. Therefore, the likelihood based method draws researchers’ attention when the CIR and noise power are unknown in the multipath fading channel. Recently, lots of scholars used the method of expectation maximization (EM) to estimate the CIR and noise power in flat fading channel in view of MPSK and MQAM signals and then they used QHLRT to identify modulation [22,23,24].Through these analysis and studies, we found EM-QHLRT can increase the identification rate effectively and make us come up with some new ideas about modulation identification in fading channel. However, EM-QHLRT is only suitable for flat fading channels, but not for underwater acoustic multipath channels. Consequently, the likelihood based methods need concern how to reduce the impact of multipath fading channels.

In [25,26], researchers used the cross-relations in the multiple reception in order to complete channel estimation in single input multi output (SIMO) channel. In [25], a least square blind channel estimation method is proposed and the necessary condition for channel estimation based on cross-relations is analyzed in SIMO channel. In [26], an adaptive blind channel estimation algorithm is proposed. However, these blind channel estimation methods have huge amount of computation to do and do not consider the sparsity of underwater acoustic channel. In recent years, scholars have done lots of research on channel estimation methods based on sparse bayesian learning (SBL) and have achieved good result [27,28,29]. However, it is difficult to apply the SBL-based channel estimation method in the non-cooperative condition due to the lack of prior information. Moreover, the identification will fail because the original constellation cannot be achieved when the channel estimation result is applied to the likelihood method directly.

In this paper, we consider likelihood-based identification of modulation types including BPSK, QPSK, and 8PSK. Firstly, we use the cross-relations and SBL to estimate the CIR. Then, EM algorithm is used to further process the received signal after channel equalization. Finally, modulation can be identified successfully by likelihood-based method. The proposed method can effectively improve the identification rate in the underwater acoustic multipath channel. The proposed SBL-EM-QHLRT method provide the following contributions in this regard:

• Aiming at the problem that the traditional likelihood-based modulation identification method cannot achieve the identification in underwater acoustic multipath channel, we propose a blind channel estimation method based on SBL, which can eliminate the impact of multipath on the identification effectively.
• Because of the error of SBL blind channel estimation and the unknown original constellation mapping of MPSK, we model the signal processed by SBL as Gaussian mixture model (GMM), and use EM in orr to correct the constellation and compensate the error of channel estimation.

The rest of the paper is organized, as follows. Section 2 shows a SIMO underwater acoustic system model, and the cross-relation in SIMO channel is derived. A Novel method named SBL-EM-QHLRT is developed in Section 3. Section 4 discusses the performance of SBL-EM-QHLRT method through numerical simulation. Conclusions are drawn in Section 5.

2. System Model

The SIMO channel model which is showed in Figure 1 can be described as follows:

$$\left\{\begin{array}{c}{\mathrm{x}}_1\left(n\right)=\mathrm{s}\left(n\right)\ast {\mathrm{h}}_1\left(n\right)+{\mathrm{w}}_1\left(n\right)\hfill \\ ⋮\hfill \\ {\mathrm{x}}_K\left(n\right)=\mathrm{s}\left(n\right)\ast {\mathrm{h}}_K\left(n\right)+{\mathrm{w}}_K\left(n\right)\hfill \end{array}\right.$$

(1)

where * represents convolution, K is the number of outputs, $w\left(n\right)={\left[{\mathrm{w}}_1\left(n\right),{\mathrm{w}}_2\left(n\right),\dots ,{\mathrm{w}}_K\left(n\right)\right]}^T$ represents the noise which obey the Gaussian distribution with mean value of zero and variance of $\sigma$. $\mathrm{s}\left(n\right)$ is the transmission signal of length N, $x\left(n\right)={\left[{\mathrm{x}}_1\left(n\right),{\mathrm{x}}_2\left(n\right),\dots ,{\mathrm{x}}_K\left(n\right)\right]}^T$ represents the received signal, $h\left(n\right)={\left[{\mathrm{h}}_1\left(n\right),{\mathrm{h}}_2\left(n\right),\dots ,{\mathrm{h}}_K\left(n\right)\right]}^T$ represents the CIR, We assume that the length of CIR is L. Without considering the influence of noise, the received signals of every output have the following cross-relationship:

$${\mathrm{x}}_i\left(n\right)\ast {\mathrm{h}}_j\left(n\right)={\mathrm{x}}_j\left(n\right)\ast {\mathrm{h}}_i\left(n\right),1\le i\ne j\le K$$

(2)

Equation (2) can be expressed as a matrix form, as shown as:

$$X_i\left(n\right){\mathrm{h}}_j\left(n\right)=X_j\left(n\right){\mathrm{h}}_i\left(n\right)$$

(3)

where $X_i\left(n\right)=\left[{\mathrm{x}}_i^1\left(n\right),{\mathrm{x}}_i^2\left(n\right),\dots ,{\mathrm{x}}_i^l\left(n\right),\dots ,{\mathrm{x}}_i^L\left(n\right)\right]$, ${\mathrm{x}}_i^l\left(n\right)={\mathrm{x}}_i(n-l)$ is the cyclic convolution matrix of size $N\times L$.

According to Equation (3):

$${\mathrm{X}}^k\left(n\right)=\left[\begin{array}{c}0\cdots 0X_{k+1}\left(n\right)-X_k\left(n\right)\cdots 0\hfill \\ ⋮\ddots ⋮⋮0\ddots ⋮\hfill \\ \underset{k-1blocks}{\underbrace{0\cdots 0}}\underset{K-k+1blocks}{\underbrace{X_K\left(n\right)0\cdots -X_k\left(n\right)}}\hfill \end{array}\right]$$

(4)

Every block of Equation (4) is a matrix of size $N\times L$. We can obtain the following, as shown as:

$$\mathrm{X}\left(n\right)\mathrm{h}\left(n\right)=0,\mathrm{X}\left(n\right)={\left[{\mathrm{X}}^1\left(n\right),\dots ,{\mathrm{X}}^{K-1}\left(n\right)\right]}^T$$

(5)

Adding constraints by setting ${\mathrm{h}}_1\left(⌊L/3⌋\right)=1$, where $⌊·⌋$ represents the operation of rounding down. The purpose of the constraints is to avoid zero solution, which will introduce the delay and gain to the final channel estimation results inevitably. However, in practice, the delay and gain will not affect the performance of the system. Equation (5) can be expressed as:

$$-\mathrm{b}=\Psi \tilde{h}$$

(6)

where $\mathrm{b}$ is the column $⌊L/3⌋$ of the matrix $\mathrm{X}\left(n\right)$, $\Psi$ is the matrix of all the columns in $\mathrm{X}\left(n\right)$ except $⌊L/3⌋$, $\tilde{h}$ is all elements except ${\mathrm{h}}_1\left(⌊L/3⌋\right)=1$ in h. The CIR obtained by Equation (6) needs to meet two identifiable condition, which is described in detail in [25]. One condition is that the transmission function of each subchannel in $h\left(n\right)={\left[{\mathrm{h}}_1\left(n\right),{\mathrm{h}}_2\left(n\right),\dots ,{\mathrm{h}}_K\left(n\right)\right]}^T$ has no common zero, and the other condition is that the autocorrelation matrix of the transmission signal $\mathrm{s}\left(n\right)$ is full rank.

3. Method of SBL-EM-QHLRT Modulation Identification

3.1. Blind Channel Estimation Based on SBL

We can achieve the process of calculating CIR by using the method of SBL [30], which is based on Equation (6), due to the sparsity of underwater acoustic channel. The specific process is as follows. We replace $-\mathrm{b}$ with t and assume that the channel has noise that obeys the Gaussian distribution with mean value of 0 and variance of $\sigma$, Equation (6) can be expressed as:

$$t=\Psi \tilde{h}+\epsilon$$

(7)

The likelihood function of t can be expressed as:

$$p(t|\tilde{h},\sigma )={(\pi \sigma )}^{-N\ast K(K-1)/2}exp(-{∥t-\Psi \tilde{h}∥}^2/\sigma )$$

(8)

The maximum a posteriori (MAP) of $\tilde{h}$ can be expressed as:

$${\tilde{h}}_{MAP}=\underset{\tilde{h}}{argmax}\left[-{∥t-\Psi \tilde{h}∥}^2/\sigma +lnp(\tilde{h})\right]$$

(9)

The first term on the right of Equation (9) comes from the likelihood function and the second term comes from the prior distribution. SBL assumes the prior distribution of parameters firstly and then calculates the maximum posterior estimate. Assuming the prior distribution of CIR, as:

$$p(\tilde{h}|\alpha )={\displaystyle \sum _{i=0}^{K\ast L-1}}CN(0,{\alpha }_i^{-1})$$

(10)

where $\alpha =[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_i,\dots ,{\alpha }_{KL-1}]$ is called hyper-parameter and corresponds to the elements in $\tilde{h}$ one by one. The value of element in $\tilde{h}$ corresponding to $\alpha$ is zero When the value of element in $\alpha$ tends to infinity. In order to calculate MAP, we need to calculate $p(\tilde{h},\alpha ,\sigma |t)$:

$$p(\tilde{h},\alpha ,\sigma |t)=p(\tilde{h}|\alpha ,\sigma ,t)p(\alpha ,\sigma |t)$$

(11)

The first term on the right of Equation (11) is posterior distribution of $\tilde{h}$. According to the Bayesian rule:

$$\begin{array}{c}p(\tilde{h}|\alpha ,\sigma ,t)=p(t|\tilde{h},\sigma )p(\tilde{h}|\alpha )/p(t|\alpha ,\sigma )=(2\pi {)}^{-ML-1}{\left|\sum \right|}^{-1/2}exp\left\{-\frac{1}{2}{\left(\tilde{h}-\mu \right)}^T\sum ^{-1}\left(\tilde{h}-\mu \right)\right\}\hfill \end{array}$$

(12)

where $\sum ={\left(A+{\sigma }^{-1}{\mathsf{\Psi }}^T\Psi \right)}^{-1}$, $A=diag({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_i,\dots ,{\alpha }_{KL-1})$ is post-covariance matrix, $\mu ={\sigma }^{-1}\sum {\Psi }^{T}t$ is posterior mean. If we get the value of $\alpha$ and $\sigma$, we can bring them into $\mu ={\sigma }^{-1}\sum {\Psi }^{T}t$ to get the MAP of $\tilde{h}$ [31]. The second term on the right of Equation (15) can be expressed as:

$$p(\alpha ,\sigma |2t)=p(t|\alpha ,\sigma )p\left(\alpha \right)p\left(\sigma \right)/p\left(t\right)$$

(13)

Assuming $\alpha$ and $\sigma$ obey the uniform distribution and we can achieve the estimation of $\alpha$ and $\sigma$ by maximizing $p\left(t\right|\alpha ,\sigma )$:

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(t\right|\alpha ,\sigma )& =\int p\left(t\right|\tilde{h},\sigma )p\left(\tilde{h}\right|\sigma )d\tilde{h}\hfill \\ \hfill & ={\left(2\pi \right)}^{-NK(K-1)/2}{\left|\sigma I+\Psi A^{-1}{\Psi }^T\right|}^{-1/2}\hfill \\ \hfill & \times \left\{-\frac{1}{2}{t}^{-T}{\left(\sigma I+\Psi A^{-1}{\Psi }^T\right)}^{-1}t\right\}\hfill \end{array}\hfill \end{array}$$

(14)

The derivatives of $p\left(t\right|\alpha ,\sigma )$ with respect to $\alpha$ and $\sigma$, while setting the derivatives to zero, we can get the new estimates for $\alpha$ and $\sigma$ as:

$${\alpha }_i^{new}=\left(1-{\alpha }_i\sum _{ii}\right)/{\mu }_i^2$$

(15)

$${\sigma }^{new}={∥t-\Psi \mu ∥}^2/\left(NK(K-1)/2-\sum _i(1-{\alpha }_i\sum _{ii})\right)$$

(16)

where ${\sum }_{ii}$ is the $ith$ diagonal element of ∑, ${\mu }_i$ is the $ith$ element of $\mu$. According to the above analysis, SBL sets the initial value of $\alpha$ and $\sigma$, threshold and iteration times. The iterative process can be showed as follow. Firstly, we calculate ∑ and $\mu$, and achieve new $\alpha$ and $\sigma$. Secondly, we set the corresponding value in $\tilde{h}$ to zero when the value in $\alpha$ is greater than the threshold. Then we repeat the iterative process for set times. Finally, $\tilde{h}$ is only a small number of non-zero values and the channel equalization is realized.

3.2. Modulation Identification Method of EM-QHLRT

The performance of likelihood based identification method will be greatly affected because of the error of channel estimation and the unknown original constellation. Therefore, EM is used to compensate the error in channel estimation and the phase. Signal after channel equalization is modelled as GMM. We assume the number of Gaussian models is M and the probability density function of GMM can be expressed as [32]:

$$f\left(r\right)={\displaystyle \sum _{m=1}^M}\frac{1}{M}CN\left(r\right|s_m,h,\sigma )$$

(17)

where $h=a{e}^{j\phi }$ is channel coefficients, $\alpha$ is amplitude attenuation, $\phi$ is frequency offset, $CN\left(r\right|s_m,h,\sigma )$ is the probability density function of the $mth$ Gaussian model, $s_m$ is the symbol of the original constellation. Step Expectation aims to estimate the likelihood value of symbols regarding to every Gaussian model. Setting the hidden variable $z_{n,m}$ to represent that symbol $r_n$ comes from the $mth$ Gaussian model.

$$z_{n,m^{\prime }}=\left\{\begin{array}{c}1\mathrm{if}{m}^{\prime }=\underset{m}{\mathrm{argmin}}{\left|r_n-h{s}_m\right|}^2\hfill \\ 0other\hfill \end{array}\right.$$

(18)

Equation (17) can be expressed as:

$$\begin{array}{c}\begin{array}{cc}\hfill f\left(r,z\left|s,h,\sigma \right.\right)& ={\displaystyle \prod _{n=1}^N}{\displaystyle \sum _{m=1}^M}\frac{1}{M}N\left(r_n\left|s_m,h,\sigma \right.\right)z_{n,m}\hfill \\ \hfill & =\frac{1}{M^N}{\displaystyle \prod _{n=1}^N}{{\displaystyle \prod _{m=1}^M}N\left(r_n\left|s_m,h,\sigma \right.\right)}^{z_{n,m}}\hfill \end{array}\hfill \end{array}$$

(19)

Taking logarithm for Equation (19):

$$lnf\left(r,z\left|s,h,\sigma \right.\right)=-NlnM+{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{z}_{n,m}lnN\left(r_n\left|s_m,h,\sigma \right.\right)$$

(20)

We define Q function as:

$$\begin{array}{c}\begin{array}{cc}\hfill Q(h,\sigma ,h^i,{\sigma }^i)& ={\mathrm{E}}_{\mathrm{z}}\left[lnf\left(r,z\left|s,h,\sigma \right.\right)\left|\mathrm{r},h^i,{\sigma }^i\right.\right]\hfill \\ \hfill & =-NlnM+{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}\mathrm{E}\left(z_{n,m}\left|\mathrm{r},h^i,{\sigma }^i\right.\right)\hfill \\ \hfill & \times \left[ln\left(\pi \sigma \right)+{\left|r_n-h{s}_m\right|}^2/\sigma \right]\hfill \end{array}\hfill \end{array}$$

(21)

where

$$\begin{array}{c}\mathrm{E}(z_{n,m}\left|r_n,h^i,{\sigma }^i\right.)=f(z_{n,m}=1\left|r_n,h^i,{\sigma }^i\right.)\hfill \\ =exp(-\frac{{\left|r_n-h^i{s}_m\right|}^2}{{\sigma }^i})/{\displaystyle \sum _{k=1}^M}exp(-\frac{{\left|r_n-h^i{s}_k\right|}^2}{{\sigma }^i})\hfill \end{array}$$

(22)

where $h^i$ and ${\sigma }^i$ represents the channel coefficients and noise power of the previous iteration, assuming the value of Equation (20) is ${\gamma }_{n,m}$. We replace h with $\mathrm{Re}\left\{h\right\}+j\mathrm{Im}\left\{h\right\}$ and Q function can be expressed as:

$$\begin{array}{c}\begin{array}{cc}\hfill Q(h,\sigma ,h^i,{\sigma }^i)& =-NlnM+{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}\hfill \\ \hfill & \times \left[ln\left(\pi \sigma \right)+\frac{{\left|r_n-(\mathrm{Re}\left\{h\right\}+j\mathrm{Im}\left\{h\right\})s_m\right|}^2}{\sigma }\right]\hfill \end{array}\hfill \end{array}$$

(23)

Step Maximization aims to maximize the Q function. The derivatives of $Q(h,\sigma ,h^i,{\sigma }^i)$ with respect to the $\mathrm{Re}\left\{h\right\}$, $\mathrm{Im}\left\{h\right\}$ and $\sigma$, while setting the derivatives to zero, we get the new estimates, as shown as:

$$\mathrm{Re}\{\widehat{h}\}=\frac{{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}\left[\mathrm{Re}\left\{r_{n,m}\right\}\mathrm{Re}\left\{s_m\right\}+\mathrm{Im}\left\{r_{n,m}\right\}\mathrm{Im}\left\{s_m\right\}\right]}{{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}{\left|s_m\right|}^2}$$

(24)

$$\mathrm{Im}\{\widehat{h}\}=\frac{{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}\left[\mathrm{Im}\left\{r_{n,m}\right\}\mathrm{Re}\left\{s_m\right\}+\mathrm{Re}\left\{r_{n,m}\right\}\mathrm{Im}\left\{s_m\right\}\right]}{{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}{\left|s_m\right|}^2}$$

(25)

$$\widehat{\sigma }=\left({\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}{\left|r_n-h{s}_m\right|}^2\right)/\left({\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{\gamma }_{n,m}\right)$$

(26)

Repeat Equations (22)–(26) until the end of the iteration. After the channel coefficient and noise power are estimated by EM, the modulation identification based on QHLRT can be expressed as:

$$f_{QHLRT}^{\left(i\right)}={\displaystyle \prod _{n=1}^N}{\displaystyle \sum _{m=1}^{M_i}}\frac{1}{M_i}\frac{1}{\pi \widehat{\sigma }}exp\left\{-{\left|r_n-\widehat{h}{s}_m^{\left(i\right)}\right|}^2/\widehat{\sigma }\right\}$$

(27)

The result can be achieved by comparing the values of (27) under BPSK, QPSK, and 8PSK:

$$\widehat{i}=\underset{i}{argmax}ln\left\{f_{QHLRT}^{\left(i\right)}\right\}$$

(28)

where $n\left\{·\right\}$ is the natural logarithm.

4. Numerical Results

4.1. Parameter Setting

The simulation uses the channel model of single input and three outputs, and then we get three different CIR of underwater acoustic channel [33]. The generation process of channel is as follows. We assume the number of paths in three channels is 6, the average multipath delay is 14 ms, and the time delay of adjacent paths obeys the exponential distribution. We can achieve six random variables which are negative exponential distributed, then we calculate the sampling point positions $\mathrm{delay}\_\mathrm{p}$ corresponding to six paths. We also assume the amplitudes of paths are Rayleigh distributed with the average power decreasing exponentially with delay, where the difference between $T\_guard=20$ ms is $p=20$ dB. Therefore, we can conclude that:

$$delta\_k=\alpha e^{(-delay\_p/(fs\ast t_{-}\left\{rms\right\}))}$$

(29)

where $\alpha =1$, $t_{-}\left\{rms\right\}=T\_guard/log({10}^{p/10})$. Now, we can achieve random variables that are Rayleigh distributed in each path, and then we assume that the value of random variable is $(a+bj)$. The amplitude of each path can be expressed as:

$$delta=\left|\sqrt{0.5\ast delta\_k}\ast (a+bj)\right|$$

(30)

where $\Vert$ is the modulus of complex number. Finally, the sparse multipath channel is generated by normalizing the amplitudes of all paths [34,35]. The noise included in simulation is the gaussian white noise in band, and a fixed frequency offset of 300 Hz is added in simulation in order to verify the phase compensation capability of the algorithm, while the constellation may rotate by down sampling, which make the receiving signal become the constellation. The sampling frequency of the communication system is 48 kHz, the carrier frequency is 8 kHz, and the bandwidth is 6 kHz. We select 100 groups of signals from BPSK, QPSK, and 8PSK, respectively, to calculate the identification rate. The simulation parameter settings are shown in

Table 1. Parameters of Underwater Acoustic Multiple Phase Shift Keying (MPSK).

sampling frequency $f_s/\mathrm{kHz}$	48
Bandwidth$/\mathrm{kHz}$	6
Carrier frequency $f_c/\mathrm{Hz}$	8
frequency offset $f_d/\mathrm{Hz}$	300

.

4.2. Performance of SBL Blind Channel Estimation

The performance of blind channel estimation of SBL is verified by simulation, the multipath length of three outputs is assumed to be 900 sampling points. We use normalized mean square error (NMSE) to evaluate the performance of SBL channel estimation.

Figure 2 shows the estimated CIR of the three channels. Figure 3 shows the difference between the original CIR and the CIR estimated by SBL. The first 900 sampling points correspond to channel 1, the 901–1800 sampling points correspond to channel 2, and the 1801–2700 sampling points correspond to channel 3. Figure 4 shows the NMSE of three channels. Refer to Figure 2, Figure 3 and Figure 4, by using SBL, we can achieve the estimation of underwater acoustic sparse channel accurately with the increase of SNR. In addition, we also simulate the performance of channel estimation under time-varying channel. We set the simulation SNR to 20 dB. After each time interval, we can achieve 100 groups of random sparse channels [21], and then we use the proposed method to estimate these channels. Finally, we take the average to be the result as it is shown in Figure 5. From Figure 5, we can find that channel estimation accuracy is low as the channels’ NMSE is higher than −10 dB when the time interval is 0.05 s or lower than 0.05 s. By contrast, the channels’ NMSE is lower than −20 dB when the time interval is higher than 0.2 s, which shows that the method proposed in this paper is still effective when the time interval is about 0.2 s. Refer to Figure 4 and Figure 5, due to the constraints that we set ${\mathrm{h}}_1\left(⌊L/3⌋\right)=1$, the gain of channel 1 achieved by channel estimation cannot be further eliminated. Therefore, we can know the estimation of channel 2 and channel 3 are more accurate than channel 1.

In this section, we verify the performance of the SBL-based blind channel estimation method. Refer to Figure 2, Figure 3, Figure 4 and Figure 5, we found that the SBL-based channel estimation method can show good performance. After we obtained the estimated CIR, the channel equalization can be used to eliminate channel interference. However, even the channel interference is eliminated, the likelihood-based method cannot effectively identify MPSK signals in the non-cooperative condition.

4.3. Performance of SBL-EM-QHLRT

Two problems need to be considered in the likelihood method. The one is that the original constellation is unknown, and the other one is that the constellation will rotate when obvious frequency offset exist in the underwater acoustic and down sampling will also make the constellation rotate even the original constellation can be known. Obviously, the likelihood method cannot be used to identify modulation directly. Therefore, we use EM to solve these problems. EM sets the original constellation firstly, and then it estimates the amplitude attenuation and phase offset of the channel under all possible modulation modes.

The simulation process is as follows. QPSK modulation is adopted in simulation and the constellation is [0.707 + 0.707j, 0.707 − 0.707j, −0.707 − 0.707j, −0.707 + 0.707j], as it is shown in Figure 6a–c show the constellation of the received signal and the constellation after SBL channel estimation and equalization, respectively. By comparison original constellation and the constellation shown in Figure 6c, we can find that constellation will rotate after SBL channel estimation and equalization. In EM, the constellation of MPSK is arbitrarily determined, we set [1, 1j, −1, −1j] as original constellation in QPSK. We also need to set original constellation in BPSK and 8PSK. Figure 6d–f show the result of constellation compensated by EM under BPSK, QPSK and 8PSK. It can be found that constellation in Figure 6c has converged to constellation in Figure 6e that the influence of phase shift has been eliminated. Furthermore, the influence of the channel is further eliminated, which makes the constellation morepolymerized. Now, these two problems have been solved and the likelihood-based method can be used to identify MPSK modulation successfully. Next, we illustrate the influence of parameters of EM algorithm on its performance.

The influence of different iteration times on the performance of EM method is analysed by this simulation. We set the number of symbols of EM to 300 and the iteration times of EM to 4, 6, 8, and 10 respectively. Figure 7 is the identification rate curve of iteration times under different SNR, and we can find the identification rate is relatively better When iteration times is 8 and 10. Besides, taking the amount of computation into consideration. Therefore, we set iteration times to 8. Based on the above analysis, we can simulate the influence of different number of symbols on EM when the iteration times is 8. We set the number of symbols to 150, 300, and 450, respectively. The Figure 8 showed that the more the symbol number is, the better the identification rate is. When considering the amount of computation, we set the number of symbols to 300 finally.

In the end, we compared the identification rate by using SBL-QHLRT with frequency offset and without frequency offset, also by using EM-QHLRT without SBL and SBL-EM-QHLRT. As shown in the Figure 9 (here we have set EM iteration times to 8 and the number of symbols to 300). SBL-QHLRT with frequency offset cannot identify MPSK modulation successfully because the constellation rotation will make the QHLRT failed. EM-QHLRT also failed to identify MPSK modulation because EM is only suit for flat fading channel. SBL-QHLRT method with known original constellation and have no frequency offset can be useful, but it will be affected by channel estimation error. We can apparently find that SBL-EM-QHLRT can overcome the problem of constellation rotation and reduce the influence of underwater acoustic fading channel. Figure 10 shows the identification rate of BPSK, QPSK, and 8PSK under SBL-EM-QHLRT, of which BPSK has the highest identification rate and 8PSK has the lowest.

The simulation results show that EM can effectively improve the identification rate of underwater acoustic MPSK. When the number of iterations and the number of symbols of EM meet the requirements, EM can further eliminate the influence of channel and solve the problem of unknown original constellation. Therefore, the likelihood-based identification method can successfully identify MPSK modulation without any prior information.

5. Conclusions

This paper developed a novel likelihood-based modulation identification algorithm for MPSK signals in underwater acoustic multipath channels. The simulation results showed that the proposed method can identify the modulation mode of MPSK well in underwater acoustic multipath channel. In the aspect of blind channel estimation, the NMSE of SBL based blind channel estimation algorithm is below −20 dB when the SNR is higher than 15 dB. In addition, EM algorithm makes the received constellation more convergent and solves the problem of unknown original constellation. In view of the simulation results, some discussions for the SBL-EM-QHLRT method are given, as follows:

In the SBL channel estimation method, the channel 2 and 3 are more accurately estimated than channel 1 due to the setting constraint conditions on channel 1. Adding constraints will lead to gain and delay in the estimated CIR. Channel 2 and channel 3 can be normalized to eliminate the gain, but channel 1 cannot eliminate the gain due to constraints, which results in poor performance of channel 1 as compared with channel 2 and 3. The estimated result of channel 1 does not affect the identification rate, because we only need to take the received signal of one of the paths for processing.

The SBL blind channel estimation method proposed in this paper can also be applied to other MPSK modulation identification methods. After eliminating the channel interference, the features that are used for identification in these methods become more apparent.

From the simulation results, when the iteration times and symbol length of EM increases, the identification rate will be better, but the computational efficiency should be considered in practical use.

Further improvements of the blind likelihood-based modulation identification method for signals in different channel environment, such as Doppler shift and MIMO, is needed in the future.