## 3 Method of SBL-EM-QHLRT Modulation Identification  
### 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 \(-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}+\varepsilon \tag{7}\]  
The likelihood function of \(t\) can be expressed as:  
\[p(t|\tilde{h},\sigma)=(\pi\sigma)^{-N*K(K-1)/2}\exp(-\left\|t-\Psi\tilde{h} \right\|^{2}/\sigma) \tag{8}\]  
The maximum a posteriori (MAP) of \(\tilde{h}\) can be expressed as:  
\[\tilde{h}_{MAP}=\operatorname*{arg\,max}_{\tilde{h}}\left[-\left\|t-\Psi\tilde {h}\right\|^{2}/\sigma+\ln p(\tilde{h})\right] \tag{9}\]  
Figure 1: Single Input Multi Output channel model. 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)=\sum_{i=0}^{KsL-1}CN(0,\alpha_{i}^{-1}) \tag{10}\]  
where \(\alpha=[\alpha_{1},\alpha_{2},\ldots,\alpha_{i},\ldots,\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) \tag{11}\]  
The first term on the right of Equation (11) is posterior distribution of \(\tilde{h}\). According to the Bayesian rule:  
\[p(\tilde{h}|\alpha,\sigma,t)=p(t|\tilde{h},\sigma)p(\tilde{h}|\alpha)/p(t| \alpha,\sigma)=(2\pi)^{-ML-1}|\Sigma|^{-1/2}\exp\left\{-\frac{1}{2}\big{(} \tilde{h}-\mu\big{)}^{T}\Sigma^{-1}\big{(}\tilde{h}-\mu\big{)}\right\} \tag{12}\]  
where \(\sum=\big{(}A+\sigma^{-1}\Psi^{T}\Psi\big{)}^{-1}\), \(A=diag(\alpha_{1},\alpha_{2},\ldots,\alpha_{i},\ldots,\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(\alpha)p(\sigma)/p(t) \tag{13}\]  
Assuming \(\alpha\) and \(\sigma\) obey the uniform distribution and we can achieve the estimation of \(\alpha\) and \(\sigma\) by maximizing \(p(t|\alpha,\sigma)\):  
\[ \begin{split} p(t|\alpha,\sigma)&=\int p(t|\tilde{h },\sigma)p(\tilde{h}|\sigma)d\tilde{h}\\ =&(2\pi)^{-NK(K-1)/2}\Big{|}\sigma I+\Psi A^{-1} \Psi^{T}\Big{|}^{-1/2}\\ &\times\left\{-\frac{1}{2}t^{-T}\Big{(}\sigma I+\Psi A^{-1}\Psi^ {T}\Big{)}^{-1}t\right\}\end{split}  \tag{14}\]  
The derivatives of \(p(t|\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}\underset{ii}{\sum}\right)/\mu_{i}^{2} \tag{15}\]  
\[\sigma^{new}=\left\|t-\Psi\mu\right\|^{2}/\left(NK(K-1)/2-\underset{i}{\sum}(1 -\alpha_{i}\underset{ii}{\sum})\right) \tag{16}\]  
where \(\underset{ii}{\sum}\) is the \(ith\) diagonal element of \(\sum\), \(\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 \(\sum\) 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.