Blind receiver design-based multi-dimensional stable distribution model using the EM algorithm

Abstract

Multi-antenna receivers are a key technology for modern communication systems. Signal attenuation in the very low frequency (VLF) channel is a serious problem. In addition to high signal attenuation, the resulting noise on the VLF channel is non-Gaussian. Hence, to analyze and address the issue of the normal operation of the multi-antenna receiver in the VLF channel, we must study the signal detection problem in a multi-dimensional non-Gaussian fading channel. Motivated by the existing blind receiver in an underwater submarine single-antenna receiver, we study the signal detection and estimation algorithm under the fading channel for a multi-dimensional non-Gaussian noise model. In this study, we propose a blind receiver based on the expectation-maximization (EM) algorithm. The proposed blind receiver can reduce non-Gaussian noise. In addition, we propose a nonlinear receiver that can accurately receive the transmitted signal over the high-attenuation VLF communication systems. Numerical and simulation results over uncorrelated and correlated non-Gaussian channels confirm that the design of the proposed blind receiver is close to optimal, with low computation complexity.

Introduction

Very low frequency (VLF) of 3–30 kHz is a sufficient frequency band for underwater communications and navigation. ¹ VLF is mainly used in submarine communication, and it has a relatively large skin depth to diffuse in the conductive water. ² VLF is used in one-way communication between a base station and a submarine. VLF is also used in remote sensing applications, ³ to collect continuous observation signals for the coverage areas. Such remote sensing observations can be performed using several models, such as the Long-Wave Propagation Capability (LWPC) program and Modefinder. ³ The application of such schemes is based on the ionospheric state and properties of the analyzed disturbances that affect the possibility of using approximation or empirical formulas. In general, the 3–30 kHz VLF has been used in early wireless telegraphy, navigation beacons and time signals, geophysical and atmospheric measurement, and military communications such as submarine communications.

In the early (1905–1925) era of wireless telegraphy, the VLF band was used for achieving long-distance communication. The era of telegraphy was based on high-power networks of low and very low frequency. Such networks were designed to use radiotelephony by employing amplitude modulation and single-sideband modulation within the band starting at 20 kHz. Due to the long communication range provided by the VLF, as well as its stable phase characteristics, the VLF band has been used for long-range hyperbolic radio navigation systems used by ships and aircrafts to determine their geographical position by comparing the phase of radio waves received from fixed VLF navigation beacon transmitters. The worldwide Omega system used frequencies from 10 to 14 kHz. In geophysical and atmospheric measurements, the VLF band is employed by geophysicists for investigating long-range lightning position as well as phenomena such as the aurora. Measurements of whistlers are used to deduce the physical properties of the magnetosphere. Geophysicists use the VLF-electromagnetic receivers to evaluate conduction near the Earth’s surface. VLF transmitters are used in military submarine applications. Due to their long-range and high dependability, it is expected that in a global war, VLF communications will be less interrupted by explosions than higher frequency communications. Because it can infiltrate ocean water, VLF is used by the navy to communicate with submarines near the surface.

In this study, we propose a blind receiver design based on a multi-dimensional stable distribution model using the Expectation-Maximization (EM) algorithm for submarine communications. The concept is that the blind receiver will reduce the complexity of the underwater submarine receiver under uncorrelated noise.

Multi-antenna receiving technology^4–10 can improve the receiver performance, particularly in VLF communication channels. The multi-antenna receiver in a submarine can overcome the single-antenna’s poor receiving capability. In addition, to achieve the best communication reception, the receiver often needs to estimate the channel parameters in advance, and then, the signal can be demodulated by equalization using the estimated channel parameters, particularly in the case of non-Gaussian noise. One estimation method is blind receiving technology, which obtains an optimal signal through the use of training sequences to estimate the channel and reduce the interference. However, due to its high computational complexity and significant hardware requirements, blind receiver detection still requires much research.

In multi-antenna receiving technology, the noise is often assumed to be a multi-dimensional Gaussian model, which is acceptable in the high-frequency channel. However, in the VLF channel, the noise distribution has the characteristics of a heavy tail and is significantly non-Gaussian.^11–15 Such non-Gaussian noise consideration is due to the influence of atmospheric thunder and lightning. To achieve optimal reception in submarine multi-antenna communications, it is necessary to study multi-antenna receiving technology in a non-Gaussian noise environment.

The multi-dimensional mixed Gaussian model and the multi-dimensional stable distribution model are commonly used as multi-dimensional non-Gaussian noise models. They can be regarded as weakened forms of Class A and Class B in the multi-dimensional Middleton noise model, ¹⁶ respectively. For signal detection under these two types of models, studies^17–19 have often focused on the structure of the receiver under known parameters. This type of receiver often requires a reserved training sequence to estimate the channel parameters, but in practice, particularly in the VLF channel, this method theoretically cannot achieve good performance due to the low communication rate and the fluctuation of channel parameters. Therefore, the channel parameters estimated using the training sequence may not be suited for subsequent signal estimation. Hence, blind non-Gaussian receivers can be incorporated in such communication systems. The blind non-Gaussian receivers can jointly estimate channel parameters and signals, with performance closer to the actual needs.

Ilow and Hatzinakos ²⁰ explained the relationship between the spherically isotropous symmetrical alpha-stable model ²¹ and the radio frequency interference created by a Poisson distributed field of interferers. In this research, the authors presume that every receiver is bordered by the same place of effective interferers and ignored the split receiver. The spherically isotropic alpha stable model is obtained under both uniform and non-uniform interference distribution that incorporates pathloss, lognormal shadowing, and Rayleigh fading. Unfortunately, this research is restricted to baseband signaling and ignored any correlation between the interferer and the antenna receiver. Delaney ²² proposed a new model as an extension of class A distribution, which was found to be useful in analyzing multiple-input multiple-output (MIMO) receivers. ²³ Unfortunately, the multi-antenna receiver proposed in previous studies is less than ideal, as it is closely limited to two antenna receivers, ²⁴ and it requires numerical dependency among those receivers.

Herein, considering that the receiver is usually a narrowband receiver in practice, we propose a blind receiver with a multi-dimensional stable distribution model based on the EM algorithm. The proposed receiver considers correlated channels as well as uncorrelated channels. To deal with the detection problem under correlated noise channels, we propose a multi-dimensional blind receiver under the distributed noise of a correlated channel. To improve the receiving model under multi-dimensional correlated noise channels, we also derive an EM algorithm with blind signal detection in uncorrelated noise channels. Simulation experiments show that the proposed blind detector has fast iteration convergence with robustness, and its low bit error rate approximates the bit error rate performance of the optimal receiver with known parameters.

The main contributions in this study can be summarized as follows:

1. This study proposes a blind receiver with a multi-dimensional stable distribution model. To achieve acceptable performance in a VLF multi-antenna receiver, a multi-dimensional non-Gaussian noise fading channel is studied under correlated and uncorrelated channels.

2. In this study, the EM algorithm is derived specifically to improve the proposed blind receiver performance over the VLF channel. The blind receiver problem is simplified by introducing additional data or hidden parameters.

3. Simulation experiments show that the blind detector has low computation complexity with robustness, and its bit error rate performance is nearly the same as that of the optimal receiver under the condition of known parameters.

The rest of this paper is structured as follows: Section “System model and formulation of the problem” presents the system model. Section “Proposed method” presents the proposed EM algorithm under uncorrelated and correlated channels. Section “Simulation and discussion” discusses the simulation. Section “Conclusions” presents the authors’ conclusions.

System model and formulation of the problem

In this section, we will discuss the multi-antenna receiver over correlated and uncorrelated channels, as well as the operation and performance of the multi-antenna receiver over these channels.

System model under the uncorrelated channel

We considered a multi-antenna receiver equipped with $R$ receiving antennas to receive the transmitted signals transmitted over a Rayleigh flat slow fading channel. The noise in such a channel is considered an additive non-Gaussian noise. Because it is a slow fading channel, it can be assumed that the channel fading coefficient is constant in each processing window. Let $A_1...,A_J$ represent the constellation diagram of the transmitter symbols, and the transmitted symbols are independently sampled $M$ times; so, we have ${\mathit{s}}_i={[s_{i1},\cdots ,s_{\mathrm{iM}}]}^T$ , $s_{\mathrm{ij}}\in \{A_1,...A_J\}$ , $i=1,...,N,j=1,...M,$ where $N$ is the total number of received symbols by the system. Assume the received complex signal is ${\mathbf{x}}_{\mathrm{ij}}={[x_{\mathrm{ij}}^1,\cdots ,x_{\mathrm{ij}}^R]}^T,$ which can be expressed as follows:

$${\mathbf{x}}_{\mathrm{ij}}=\mathit{a}{s}_{\mathrm{ij}}+{\mathit{n}}_{\mathrm{ij}},i=1,\cdots ,N,j=1,\cdots ,M,$$

(1)

where $\mathit{a}={[a_1,\cdots ,a_R]}^T$ is the channel fading coefficient, $s_{\mathrm{ij}}$ is the single sample values, and ${\mathit{n}}_{\mathrm{ij}}={[n_{\mathrm{ij}}^1,\cdots ,n_{\mathrm{ij}}^R]}^T$ is the multi-dimensional non-Gaussian noise. In this study, let ${\mathit{n}}_{\mathrm{ij}}$ be an independently distributed variable that obeys the complex sub-Gaussian distribution. Then, according to the definition of the complex Gaussian distribution, ${\mathit{n}}_{\mathrm{ij}}$ can be rewritten as

$${\mathit{n}}_{\mathrm{ij}}=v_{\mathrm{ij}}^{1/2}{\mathit{G}}_{\mathit{ij}},$$

(2)

where ${\mathit{G}}_{\mathit{ij}}$ is the complex Gaussian vector, which obeys the multi-dimensional complex Gaussian distribution, and can be written as

$${\mathit{G}}_{\mathrm{ij}}~\mathcal{C}\mathcal{N}(0,\sum ),$$

(3)

and $v_{\mathrm{ij}}$ is a stable distribution variable of $\alpha /2$ that can be written as

$$v_{\mathrm{ij}}~f_{\alpha /2}({\left(\cos \frac{\pi \alpha }{4}\right)}^{2/\alpha },1,0),$$

(4)

Therefore, we have

$${\mathit{n}}_{\mathrm{ij}}~\mathcal{C}\mathcal{N}(0,v_{\mathrm{ij}}\sum ).$$

(5)

According to the definition of sub-Gaussian distribution, each dimension of the random variable subject to sub-Gaussian distribution has the same parameter $\alpha$ . Therefore, the estimation of parameter $\alpha$ can be reduced to the estimation $\alpha$ under the one-dimensional stable distribution. Because $\alpha$ value changes slowly in practice, the value of $\alpha$ in each dimension can be estimated by a fast estimation algorithm ²⁵ prior to the signal transmission. Each estimated value is stored as ${\alpha }_i$ , and the final estimated value of $\alpha$ can be obtained by averaging storage ${\alpha }_i$ .

System model under the correlated channel

The correlated channel assumptions in this study are as follows. First, we assume that there is a received antenna in a multi-antenna receiver, and the $R$ irradiation areas between these antennas may overlap; so, they are not necessarily independent. Second, signals are transmitted over a Rayleigh flat slow fading channel. Because it is a slow fading channel, it can be assumed that the channel fading coefficient is constant in each processing window. Third, the channel noise is considered as an additive non-Gaussian noise. Fourth, in addition to the possible correlation between each data dimension of the sampling vector in received noise, the sampled noise vectors may also correlate with the transmission of each symbol. These assumptions were also considered in previous research.^26–29

Let $A_1,\cdots ,A_J$ represent the constellation diagram of the transmitter symbol. Assume that the transmitted symbol is independently sampled $M$ times and that its transmitted signal vector ${\mathit{A}}_j={[A_{1,j},\cdots ,A_{\mathrm{mj}}]}^T,j=1,\cdots ,J$ , has a dimension of $M\times 1$ . Suppose that the received signal vector ${\mathit{x}}_{\mathrm{ir}}={[{\mathit{x}}_{\mathrm{ir},1},\cdots ,{\mathit{x}}_{\mathrm{ir},M}]}^T$ is received by r-th receiver antenna in each interval transmitted symbol. This can be written as

$${\mathit{x}}_{\mathrm{ir}}=a_r{\mathit{S}}_i+{\mathit{n}}_{\mathrm{ir}},i=1,\cdots ,N,r=1,\cdots ,R,$$

(6)

where ${\mathit{S}}_i$ is the transmitted signal vector, and its values are $\{A_1,\cdots ,A_J\}$ . The $a_r$ represents the complex attenuation coefficient of r-th receiver antenna, and ${\mathit{n}}_{\mathrm{ir}}$ is $M\times 1$ dimension noise vector, where ${\mathit{n}}_{\mathrm{ir}}={[n_{\mathrm{ir},1},\cdots ,n_{\mathrm{ir},M}]}^T$ . By combining the observation equations from $R$ received antennas, we can obtain

$${\mathit{X}}_i=\mathit{a}\otimes {\mathit{S}}_i+{\mathit{n}}_i,n=1,\cdots ,N,$$

(7)

where ${\mathit{X}}_i={[{\mathit{x}}_{i1}^T,\cdots ,{\mathit{x}}_{\mathrm{iR}}^T]}^T$ , $\mathit{a}={[{\mathit{a}}_1,\cdots ,{\mathit{a}}_R]}^T$ , and ${\mathit{n}}_i={[{\mathit{n}}_{i1}^T,\cdots ,{\mathit{n}}_{\mathrm{iR}}^T]}^T$ . The ⊗ is the Kronecker product operation. In additive noise component ${\mathit{n}}_i$ , there are correlated components with RM× 1 dimension. In this study, we assume that these correlated components obey the complex sub-Gaussian distribution with RM× 1 dimension as

$${\mathit{n}}_i=v_i^{1/2}{\mathit{G}}_i,$$

(8)

where ${\mathit{G}}_i$ is the complex Gaussian vector with RM× 1 dimension, and this vector follows the multi-dimensional complex Gaussian distribution as

$${\mathit{G}}_i~\mathcal{C}\mathcal{N}(0,\sum ),$$

(9)

and $v_i$ is the stable distribution variable of $\alpha /2$ and can be represented as equation (4). Finally, the correlated noise over the correlated channel can be written as

$${\mathit{n}}_i~\mathcal{C}\mathcal{N}(0,v_i\sum ),v_i~f_{\alpha /2}({\left(\cos \frac{\pi \alpha }{4}\right)}^{2/\alpha },1,0).$$

(10)

Proposed method

The EM algorithm^30–32 is a common maximum likelihood estimation method. It is very effective in estimating the potential probability distribution under incomplete data. It is often very difficult to use the Newton method and other methods,^33–35 but the EM algorithm is intuitive. The EM algorithm consists of two steps at each iteration: solving the expectation (E-Step) and solving the maximum value (M-Step). First, there are indeed missing values in data due to limited observation conditions. Second, because it is difficult to directly solve the likelihood function, additional data or hidden parameters are introduced to simplify the problem. This article focuses on the second EM category.

EM algorithm under uncorrelated channels

Overall EM steps under uncorrelated channels

Consider $\mathit{X}=\{{\mathbf{x}}_{\mathrm{ij}}|i=1,\cdots ,N,j=1,\cdots ,M\}$ as incomplete or missing data; if the assistant variable $v_{\mathrm{ij}}$ is introduced, we can consider $\mathbf{Y}=\{({\mathbf{x}}_{\mathrm{ij}},v_{\mathrm{ij}})|i=1,\cdots ,N,j=1,\cdots M\}$ as completed data. In this study, we assume a transmitted symbol sample set $\mathit{S}=\{{\mathit{s}}_i|i=1,\cdots ,N\}$ with $\Phi =\{\mathit{a},\sum ,\mathit{S}\}$ estimated parameters. Suppose that the posterior probability of $\mathbf{Y}$ is

$$\begin{array}{c}f(\mathbf{Y}|\mathbf{\Phi })=\underset{i=1}{\overset{N}{\Pi }}\underset{j=1}{\overset{M}{\Pi }}\frac{1}{{\pi }^R|v_{\mathrm{ij}}\sum |}\hfill \\ \exp \{-v_{\mathrm{ij}}^{-1}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})\}P\left(v_{\mathrm{ij}}\right).\hfill \end{array}$$

(11)

By evaluating the logarithm of equation (11), we obtain

$$\ln f\left(\mathbf{Y}|\mathbf{\Phi }\right)\propto -{\sum }_{i=1}^L{\sum }_{j=1}^M\ln |\sum |-\sum _{i=1}^L\sum _{j=1}^M{v}_{\mathrm{ij}}^{-1}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}}).$$

(12)

The EM algorithm maximizes $\mathcal{Q}(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }})$ to obtain a maximum likelihood estimation in two steps:

E-Step: $\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathbf{\Phi }}^{\left(t\right)},\mathit{X}\}$ ;

M-Step: solve $\mathbf{\Phi }=\arg \underset{\mathbf{\Phi }}{max}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})$ and let ${\mathbf{\Phi }}^{(t+1)}={\mathbf{\Phi }}^{\left(t\right)}$ ,

where $\mathcal{Q}(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }})$ is given by

$$\mathcal{Q}\left(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }}\right)=E\left\{\ln f\left(\mathbf{Y}|\mathbf{\Phi }\right)|\stackrel{\prime }{\mathbf{\Phi }},\mathit{X}\right\}\propto -{\sum }_{i=1}^L{\sum }_{j=1}^M\ln |\sum |-\sum _{i=1}^L\sum _{j=1}^{M}E\left(v_{\mathrm{ij}}^{-1}\right){({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}}).$$

(13)

In the M-Step, it is necessary to solve $(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})$ to obtain the maximum value of $\mathbf{\Phi }$ . However, if $\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})$ directly solves the derivation of $\mathbf{\Phi }$ , the computation complexity will be very high. Therefore, to reduce the computation complexity of the traditional EM, in this study, the proposed EM algorithm is based on a variant of the space-altering generalized expectation-maximization (SAGE) algorithm.^36,37 In the SAGE-based EM algorithm, one of the variables is solved, while the remaining variables are constant, thereby reducing the complexity. A possible step-by-step iteration of the SAGE-EM algorithm is

1. E-Step: Solve $\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathbf{\Phi }}^{\left(t\right)},\mathit{X}\}$ ;

2. M-Step: Fixed ${\mathit{a}}^{\left(t\right)}$ and ${\mathit{s}}_i^{\left(t\right)}$ to solve $\sum _{}^{\left(t+1\right)}=\arg \underset{\sum }{\max }\mathcal{Q}\left(\sum |{\sum }^{\left(t\right)}\right)$ ;

3. E-Step: $\mathcal{Q}(\mathbf{\Phi }|{\mathit{a}}^{\left(t\right)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathit{a}}^{\left(t\right)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)},\mathit{X}\}$ ;

4. M-Step: Fixed $\sum ^{(t+1)}$ and ${\mathit{s}}_i^{\left(t\right)}$ to solve ${\mathit{a}}^{(t+1)}=\arg \underset{a}{max}\mathcal{Q}(\mathit{a}|{\mathit{a}}^{\left(t\right)})$ ;

5. E-Step: Solve $\mathcal{Q}(\mathbf{\Phi }|{\mathit{a}}^{(t+1)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathit{a}}^{(t+1)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)},\mathit{X}\}$ ;

6. M-Step: Fixed ${\mathit{a}}^{(t+1)}$ and $\sum ^{(t+1)}$ to solve ${\mathit{s}}_i^{(t+1)}=\arg \underset{s_i}{max}\mathcal{Q}({\mathit{s}}_i|{\mathit{s}}_i^{\left(t\right)})$ .

E-Step

The main task of the E-Step EM algorithm is calculated by $E(v_{\mathrm{ij}}^{-1})$ . The posterior probability $v_{\mathrm{ij}}$ can be calculated as follows:

$$\begin{array}{c}P(v_{\mathrm{ij}}|\cdots )\propto \hfill \\ \frac{1}{{\pi }^R|v_{\mathrm{ij}}\sum |}\exp \{-v_{\mathrm{ij}}^{-1}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})\}\hfill \\ f_{\alpha /2}(v_{\mathrm{ij}};\cos {\left(\frac{\pi \alpha }{4}\right)}^{2/\alpha },1,0).\hfill \end{array}$$

(14)

The expression in equation (14) is complicated, and the stable distribution does not have the probability density of the closed-form; so, it is difficult to obtain $E\left(v_{\mathrm{ij}}^{-1}\right)$ by analytical form. However, we can obtain its expectation as follows:

$$E\left(v_{\mathrm{ij}}^{-1}\right)\approx \frac{1}{L}\sum _{l=0}^L{v}_{\mathrm{ij}}^{\left(l\right)},$$

(15)

where $v_{\mathrm{ij}}^{\left(l\right)}$ are the sampling particles obeying the formula in equation (14), and $L$ is the total number of particles. Hence, a series of particles $v_{\mathrm{ij}}^{\left(l\right)}$ can be generated by the Monte Carlo method, and then, the expectation $E\left(v_{\mathrm{ij}}^{-1}\right)$ can be calculated as an average.

Because the $v_{\mathrm{ij}}$ distribution is not common, it is difficult to directly generate the random numbers $v_{\mathrm{ij}}$ by sampling. Thus, rejection sampling is considered in this study. For the R-dimensional complex Gaussian distribution, equation (14) is upper-bounded by

$$\begin{array}{c}\frac{1}{{\pi }^R|v_{\mathrm{ij}}\sum |}\exp \{-v_{\mathrm{ij}}^{-1}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})\}\hfill \\ \le u_b={|\frac{\pi }{M}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H\sum ^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})\sum |}^{-1}\exp \{-1/R\}.\hfill \end{array}$$

(16)

Therefore, $v_{\mathrm{ij}}$ can be generated by rejection sampling, where the steps of generating random numbers in each $v_{\mathrm{ij}}$ can be as follows:

1. Sample at $\alpha /2$ with stable distribution and generate ${\lambda }_{\mathrm{ij}}^{\left(l\right)}$ as a stable distribution variable of $\alpha /2$ , as shown in equation (4).

2. Sample based on a uniform distribution $u^{\left(l\right)}~U(0,u_b)$ and then generate $u^{\left(l\right)}$ ;

3. If $u^{\left(l\right)}\le \mathcal{C}\mathcal{N}({\mathbf{x}}_{\mathrm{ij}},\mathit{a}{s}_{\mathrm{ij}},v_{\mathrm{ij}}\sum )$ , accept ${\lambda }_{\mathrm{ij}}^{\left(l\right)}$ ; otherwise, go back to step (1).

However, rejection sampling has one drawback: when the probability of rejection is too large, random number generation will be prohibitively slow. Therefore, importance sampling is used here. For $E\left(v_{\mathrm{ij}}^{-1}\right)$ , $P\left(v_{\mathrm{ij}}\right)$ is decomposed into $P\left(v_{\mathrm{ij}}\right)=h\left(v_{\mathrm{ij}}\right)f_{\alpha /2}\left(v_{\mathrm{ij}}\right)$ , and we obtain

$$\begin{array}{c}E\left(v_{\mathrm{ij}}^{-1}\right)=\int v_{\mathrm{ij}}^{-1}P\left(v_{\mathrm{ij}}\right)d{v}_{\mathrm{ij}}\hfill \\ =\int v_{\mathrm{ij}}^{-1}h\left(v_{\mathrm{ij}}\right)f_{\alpha /2}\left(v_{\mathrm{ij}}\right)d{v}_{\mathrm{ij}}\approx \frac{1}{W_t}\sum _{l=1}^L{\left(q_{\mathrm{ij}}^{\left(l\right)}\right)}^{-1}w\left(q_{\mathrm{ij}}^{\left(l\right)}\right),\hfill \end{array}$$

(17)

where $q_{\mathrm{ij}}^{\left(l\right)}$ are random particles generated according to the proposed distribution ${\lambda }_{\mathrm{ij}}~f_{\alpha /2}(\cos {\left(\frac{\pi \alpha }{4}\right)}^{2/\alpha },1,0)$ , and $W_t=\sum _{l=1}^{L}w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ , where $w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)=h\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ . Because it will be normalized at the end, $w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ can be in the same proportion as $h\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ . This feature is very convenient in practical applications because $h\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ is often difficult to obtain in practice, but it is easy to remove the constant term formula. For example, the proportional equation of $h\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ in this study can be $\mathcal{C}\mathcal{N}({\mathbf{x}}_{\mathrm{ij}};\mathit{a}{s}_{\mathrm{ij}};v_{\mathrm{ij}}\sum )$ , that is, $w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)=({\mathbf{x}}_{\mathrm{ij}};\mathit{a}{s}_{\mathrm{ij}};v_{\mathrm{ij}}\sum )$ . The following is the importance sampling algorithm:

1. Generate $q_{\mathrm{ij}}^{\left(l\right)}$ from the proposed distribution $q_{\mathrm{ij}}~f_{\alpha /2}(\cos {\left(\frac{\cos \pi \alpha }{4}\right)}^{2/\alpha },1,0)$ , where $l=1,\cdots ,L$ ;

2. Calculate the weight of importance $w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)$ ;

3. Calculate $E\left(v_{\mathrm{ij}}^{-1}\right)$ according to equation (17).

M-Step

The M-Step is used to estimate the value of ∑, $\mathit{a}$ , and ${\mathit{s}}_i$ . To facilitate the symbolic representation, the old value $\mathbf{\Phi }$ and the new value $\stackrel{\prime }{\mathbf{\Phi }}$ are not distinguished in equation (18).

To estimate the value of the Hermitian matrix ∑, we need to solve ${\sum }^{†}=\arg \underset{\sum }{max}\mathcal{Q}(\sum |\sum ^{\prime })$ , that is, $\frac{\partial \mathcal{Q}}{\partial {\sum }^{-1}}=0$ . According to the derivation of the matrix:

$$\frac{\partial \log \left|{\sum }^{-1}\right|}{\partial {\sum }^{-1}}=\sum ^T,\frac{\partial x^H{\sum }^{-1}y}{\partial {\sum }^{-1}}={\mathbf{x}}^{*}{\mathit{y}}^T.$$

(18)

Then, the new estimation ∑ is given by

$${\sum }^{†}=\frac{\sum _{i=1}^N\sum _{j=1}^{M}E\left\{{\mathit{v}}_{\mathrm{ij}}^{-1}\right\}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}}){({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H}{\mathrm{NM}}.$$

(19)

To estimate the channel fading coefficient $\mathit{a}$ , we need to solve ${\mathit{a}}^{†}=\arg \underset{a}{max}\mathcal{Q}(\mathit{a}|\stackrel{\prime }{\mathit{a}})$ , that is, $\frac{\partial \mathcal{Q}}{\partial \mathit{a}}=0$ . Note that the latest estimate $\sum ^{†}$ uses ∑ at this time. According to the chain rule,

$$\begin{array}{c}\frac{\partial {({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H{\sum }^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}{\partial \mathit{a}}\hfill \\ =\frac{\partial {({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^T\partial {({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H{\sum }^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}{\partial \mathit{a}\partial ({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}.\hfill \end{array}$$

(20)

According to the nature of the conjugate gradient,

$$\frac{\partial {\mathbf{x}}^H\mathit{Ax}}{\partial \mathbf{x}}={\mathit{A}}^T{\mathbf{x}}^{*}.$$

(21)

Therefore, we have

$$\begin{array}{c}\frac{\partial {({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H{\sum }^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}{\partial \mathit{a}}\hfill \\ =-s_{\mathrm{ij}}^T{\left({\sum }^{-1}\right)}^T{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^{*}.\hfill \end{array}$$

(22)

In the above formula, * represents a conjugate operator. Therefore, the new estimate of $\mathit{a}$ is obtained as follows:

$${\mathit{a}}^{†}=\frac{\left(\sum _{i=1}^N\sum _{j=1}^{M}E\left\{v_{\mathrm{ij}}^{-1}\right\}{s}_{\mathrm{ij}}^H{\sum }^{-1}{\mathbf{x}}_{\mathrm{ij}}\right)}{\left(\sum _{i=1}^N\sum _{j=1}^{M}E\left\{v_{\mathrm{ij}}^{-1}\right\}{s}_{\mathrm{ij}}^H{\mathit{s}}_{\mathrm{ij}}{\sum }^{-1}\right)}.$$

(23)

For $s_{\mathrm{ij}}$ , we only need to find the appropriate value of $s_{\mathrm{ij}}$ to maximize $\mathcal{Q}$ . Therefore, we have

$$\begin{array}{c}{s}_{\mathrm{ij}}^{†}=\arg \underset{s_{\mathrm{ij}}}{max}\left\{\sum _{j=1}^{M}E\left\{v_{\mathrm{ij}}^{-1}\right\}{({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})}^H{\sum }^{-1}({\mathbf{x}}_{\mathrm{ij}}-\mathit{a}{s}_{\mathrm{ij}})\right\},\hfill \\ s_{\mathrm{ij}}\in \{A_1,\cdots ,A_J\}.\hfill \end{array}$$

(24)

Proposed EM algorithm under correlated channels

Overall EM steps under correlated channels

The received data $\mathbf{X}=\{{\mathbf{X}}_i|i=1,\cdots ,N\}$ in equation (7) are regarded as a type of incomplete data. If the variable $v_i$ is introduced, then $\mathbf{Y}=\{({\mathbf{X}}_i,v_i)|i=1,\cdots ,N\}$ can be regarded as complete data. The posterior probability distribution of $\mathbf{Y}$ is

$$f\left(\mathbf{Y}|\mathbf{\Phi }\right)={\prod }_{i=1}^N\frac{1}{{\pi }^{RM}\left|v_i{\sum }^{}\right|}\exp \{-v_i^{-1}{({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H\sum ^{-1}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)\}P\left(v_i\right),$$

(25)

where $\mathbf{\Phi }=\{\mathit{a},\sum ,{\mathit{S}}_i\}$ , and the logarithmic form of equation (25) is

$$\begin{array}{c}\ln f(\mathbf{Y}|\mathbf{\Phi })\propto -\sum _{i=1}^N\hfill \\ \ln |\sum |-\sum _{i=1}^N{v}_i^{-1}{({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H\sum _{}^{-1}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i).\hfill \end{array}$$

(26)

The EM algorithm maximizes $\mathcal{Q}(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }})$ to obtain the maximum likelihood estimation, where $\mathcal{Q}(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }})$ is

$$\begin{array}{c}\mathcal{Q}(\mathbf{\Phi }|\stackrel{\prime }{\mathbf{\Phi }})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|\stackrel{\prime }{\mathbf{\Phi }},\mathbf{X}\}\propto \sum _{i=1}^N\hfill \\ \ln |\sum |-\sum _{i=1}^{N}E\left(v_i^{-1}\right){({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H\sum _{}^{-1}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i).\hfill \end{array}$$

(27)

The main steps are:

E-Step: $\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathbf{\Phi }}^{\left(t\right)},\mathbf{X}\}$ ;

M-Step: Solve $\mathbf{\Phi }=\arg \underset{\mathbf{\Phi }}{max}\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})$ and let ${\mathbf{\Phi }}^{(t+1)}=\mathbf{\Phi }$ .

Because $\mathbf{\Phi }$ has multiple parameters, we use the SAGE algorithm to iterate computations. The t-th iteration calculation is as follows:

1. E-Step: Solve $\mathcal{Q}(\mathbf{\Phi }|{\mathbf{\Phi }}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathbf{\Phi }}^{\left(t\right)},\mathbf{X}\}$ ;

2. M-Step: Fixed ${\mathit{a}}^{\left(t\right)}$ and ${\mathit{S}}_i^{\left(t\right)}$ , and solve $\sum _{}^{(t+1)}=\arg \underset{\sum }{max}\mathcal{Q}(\sum |{\sum }^{\left(t\right)})$ ;

3. E-Step: Solve $\mathcal{Q}(\mathbf{\Phi }|{\mathit{a}}^{\left(t\right)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathit{a}}^{\left(t\right)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)},\mathbf{X}\}$ ;

4. M-Step: Fixed $\sum ^{(t+1)}$ and ${\mathit{S}}^{\left(t\right)}$ , and solve ${\mathit{a}}^{(t+1)}=\arg \underset{\mathit{a}}{max}\mathcal{Q}(\mathit{a}|{\mathit{a}}^{\left(t\right)})$ ;

5. E-Step: Solve $\mathcal{Q}(\mathbf{\Phi }|{\mathit{a}}^{(t+1)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)})=E\{\ln f(\mathbf{Y}|\mathbf{\Phi })|{\mathit{a}}^{(t+1)},{\sum }^{(t+1)},{\mathit{S}}^{\left(t\right)},\mathbf{X}\}$ ;

6. M-Step: Fixed $\sum ^{(t+1)}$ and ${\mathit{a}}^{(t+1)}$ , and solve ${\mathit{S}}^{(t+1)}=\arg \underset{s_i}{max}\mathcal{Q}({\mathit{S}}_i|{\mathit{S}}_i^{\left(t\right)})$ .

E-Step

The calculated E-Step is similar to that in Section “EM algorithm under uncorrelated channels.” However, because the system model condition in this section is under the correlated noise channel, for the sake of completeness, the calculation of $E\left\{v_i^{-1}\right\}$ is simply derived again. It is not difficult to obtain the posterior probability of $v_i$ as

$$\begin{array}{c}P(v_i|\cdots )\propto \hfill \\ \frac{v_i^{-\mathrm{RM}}}{{\pi }^{\mathrm{RM}}|\sum |}\exp \{-v_i^{-1}{({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H\sum ^{-1}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)\}\hfill \\ f_{\alpha /2}(v_i;\cos {(\pi \alpha /4)}^{2/\alpha },1,0).\hfill \end{array}$$

(28)

Similarly, using the Monte Carlo method to solve $E\left(v_i^{-1}\right)$ , for expectations, $E\left(v_i^{-1}\right)\approx \frac{1}{L}\sum _{l=1}^L{v}_i^{\left(l\right)}$ . For $\mathrm{RM}\times 1$ multi-dimensional complex Gaussian distribution, there is an upper bound of

$$\begin{array}{c}\frac{v_i^{-\mathrm{RM}}}{{\pi }^{\mathrm{RM}}|\sum |}\exp \{-v_i^{-1}{\mathit{\omega }}^H{\sum }^{-1}\mathit{\omega }\}\hfill \\ \le \mathcal{K}\left(\mathit{\omega }\right)={|\frac{\pi }{\mathrm{RM}}{\mathit{\omega }}^H{\sum }^{-1}\mathit{\omega }\sum |}^{-1}\exp \{-\frac{1}{\mathrm{RM}}\}.\hfill \end{array}$$

(29)

Therefore, $v_i$ can usually be generated by rejection sampling, where the steps of generating random numbers each time are as follows:

1. Generate $v_i^{\left(l\right)}$ from the proposed distribution $v_i~f_{\alpha /2}(\cos {(\pi \alpha /4)}^{2/\alpha },1,0)$ ;

2. Generate $u^{\left(l\right)}$ samples from the uniform distribution of

$$u^{\left(l\right)}~U(0,\mathcal{K}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)).$$

(30)

• (3) If $u^{\left(l\right)}\ge \mathcal{C}\mathcal{N}({\mathbf{X}}_i;\mathit{a}\otimes {\mathit{S}}_i,v_i\sum )$ , return to the step (1); otherwise, accept $v_i^{\left(l\right)}$ .

Furthermore, when the rejection probability is too large, the rejection sampling algorithm is time-consuming. Therefore, the following importance sampling algorithm can also be used:

1. Generate $q_{\mathrm{ij}}^{\left(l\right)}$ , $l=1,\cdots ,L$ , from the proposed distribution $q_{\mathrm{ij}}~f_{\alpha /2}(\cos {(\pi \alpha /4)}^{2/\alpha },1,0)$ .

2. Calculate the weight of importance as $w\left(q_{\mathrm{ij}}^{\left(l\right)}\right)=\mathcal{C}\mathcal{N}({\mathbf{X}}_i;\mathit{a}\otimes {\mathit{S}}_i,q_{\mathrm{ij}}^{\left(l\right)}\sum ),l=1,\cdots ,L.$

3. Calculate $E\left(v_{\mathrm{ij}}^{-1}\right)$ as follows:

$$E\left(v_{\mathrm{ij}}^{-1}\right)\approx \frac{1}{W_t}\sum _{l=1}^L{\left(q_{\mathrm{ij}}^{\left(l\right)}\right)}^{-1}w\left(q_{\mathrm{ij}}^{\left(l\right)}\right),$$

(31)

where $W_t=\sum _{l=1}^L{\left(q_{\mathrm{ij}}^{\left(l\right)}\right)}^{-1}.$

M-Step

As in Section “M-Step,” the M-Step is used to estimate the value of ∑, $\mathit{a},$ and ${\mathit{s}}_i$ . The new estimated value ∑ is given by

$${\sum }^{†}=\frac{\sum _{i=1}^{N}E\left\{{\mathit{v}}_i^{-1}\right\}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i){({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H}{N}.$$

(32)

Let $\sum =\sum ^{†},$ and we can estimate the fading coefficient $\mathit{a}$ . We can also obtain $\partial \mathcal{Q}/\partial \mathit{a}=0$ and perform the following transformation:

$$\mathit{a}\otimes {\mathit{S}}_i=({\mathit{I}}_{R\times R}\otimes {\mathit{S}}_i)\mathit{a}={\mathbf{\Theta }}_{\mathit{i}}\mathit{a},$$

(33)

where ${\mathbf{I}}_{R\times R}$ is an $R\times R$ identity matrix. By combining the above formula and using the chain rule of differentiation, we can obtain

$$\frac{\partial {\left({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a}\right)}^H{\sum }_{}^{-1}\left({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a}\right)}{\partial \mathit{a}}=\frac{\partial {({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a})}^T}{\partial \mathit{a}}\frac{\partial {({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a})}^H\sum _{\}}^{-1}({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a})}{\partial ({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a})}.$$

(34)

We also obtain

$$\begin{array}{l}\frac{\partial {\left({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a}\right)}^H\sum _{}^{-1}\left({\mathbf{X}}_i-{\mathbf{\Theta }}_i\mathit{a}\right)}{\partial \mathit{a}}\\ =-{\mathbf{\Theta }}_i^{\mathit{T}}{\left({\sum }^{-1}\right)}^T{\left({\mathbf{X}}_{\mathit{i}}-{\mathbf{\Theta }}_{\mathit{i}}\mathit{a}\right)}^{\ast }.\end{array}$$

(35)

Furthermore, in order to solve according to the nature of the conjugate gradient, we can obtain the channel fading coefficient $\mathit{a}$ as follows:

$${\mathit{a}}^{†}=\frac{\left(\sum _{i=1}^{N}E\left\{v_i^{-1}\right\}{\mathbf{\Theta }}_i^H{\sum }^{-1}{\mathit{X}}_i\right)}{\left(\sum _{i=1}^{N}E\left\{v_i^{-1}\right\}{\mathbf{\Theta }}_i^H{\sum }^{-1}{\mathbf{\Theta }}_i\right)}.$$

(36)

Finally, ${\mathit{S}}_i,i=1,\cdots ,N,$ can be estimated by using $\sum =\sum ^{†}$ and $\mathit{a}={\mathit{a}}^{†}$ , as follows:

$$\begin{array}{c}{\mathit{S}}_i=\arg \underset{{\mathit{S}}_i}{max}\{-{({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)}^H{\sum }^{-1}({\mathbf{X}}_i-\mathit{a}\otimes {\mathit{S}}_i)\},\hfill \\ {\mathit{S}}_i\in \{{\mathit{A}}_1,\cdots ,{\mathit{A}}_J\}.\hfill \end{array}$$

(37)

Simulation and discussion

To test the performance of the receiver described in Section “Proposed method,” a simulation experiment was conducted over the correlated and uncorrelated channels. The underwater channel was considered uncorrelated in this simulation due to the underwater ambient noise and propagation features of the uncorrelated ambient noise. It was observed that, when $M=1$ , the derivation algorithm was the same under correlated and uncorrelated channels, which bolsters the theory explained in Section “Proposed method.” Therefore, to avoid the redundancy length caused by the simulation under the two channels, this study only considers the simulation under the $M=1$ scenario, and the simulation results under other $M$ values are similar. In the simulation experiment, we assume that the number of antennas in the multi-antenna receiver is $R=2$ , which is also common in the low-frequency receiver (such as the observation platform for global atmospheric noise established by Stanford University, which used the east-west and north-south dual-antenna receiving format; the submarine can use the electric antenna and the magnetic antenna to form the orthogonal antenna array^38–42). The transmitted signal is modulated by binary phase-shift keying (BPSK) or minimum-shift keying (MSK) and the phase ambiguity is overcome by differential coding and differential decoding. In the simulation, the initial value ∑ is calculated as follows: cut or limit the large pulse amplitude data to 20% in received data, that is, perform approximate Gaussian processing. Then, calculate its covariance matrix according to multi-dimensional Gaussian noise. After the approximate Gaussian operation, the signal is optimally detected under Gaussian noise. The initial value of channel coefficient $\mathit{a}$ is simply set as ${[0.1;0.1]}^T$ . Because the iterative algorithm in the E-step is needed to estimate $E\left(v_{\mathrm{ij}}^{-1}\right)$ using the Monte Carlo algorithm, the Monte Carlo algorithm needs to set up a Monte Carlo particle number $L$ . In the simulation, through many simulation tests, we find that, when $L$ increases to $L=80$ , the estimator’s performance improvement is more obvious, but the improvement is slight when $L$ increases above 80. Taking the estimator accuracy, complexity, and redundancy into consideration, $L$ is set to be $L=100$ in this study.

In the simulation, we first tested the iterative performance of the estimator. In the measured atmospheric noise, the value of $\alpha$ is generally between 1.6 and 2.0; so, the value of $\alpha$ is set at $\alpha =1.8$ . Let $\mathit{a}={[1+i;1+i]}^T$ , $\sum =[1,0.5+0.5i;0.5-0.5i,1]$ . Figures 1 and 2 show the typical iterative convergence performance diagram of parameters $\mathit{a}$ and ∑, respectively. The relative error ordination is $|\hat{\theta }-\theta |/\left|\theta \right|$ , where $\theta$ represents the true value, and $\hat{\theta }$ represents the iterative value.

Figure 1.

(a) Iteration convergence performance of a1 and (b) Iteration convergence performance of a2.

Figure 2.

Iteration convergence performance of Σ (we use Σ(1,1) and Σ(1,2) as examples): (a) iteration convergence performance of Σ(1,1) and (b) iteration convergence performance of Σ(1,2).

Based on Figures 1 and 2, the iterative convergence rate of the estimation period is rapid, and the convergence can be achieved within five iterations. In addition, we compared the iterative performance of the estimator under a different number of symbols. When increases were from 100 to 300 and from 300 to 500, and then for 100 times, the estimated results were averaged, as shown in Table 1. It can be seen that the parameter estimation results become increasingly accurate as the number of symbols increases. As shown in Table 2, the normalized standard deviation of estimated parameters is compared in 100 simulation experiments. The normalized standard deviation is defined as follows:

$$d\theta =\frac{\sqrt{\overline{{|\theta -\overline{\theta }|}^2}}}{\left|\overline{\theta }\right|}.$$

(38)

As shown in Tables 1 and 2, with the continuous increase of the symbols $N$ , the normalized standard deviation of each estimator parameter decreases continuously, as the performance of the parameter estimator increases. When $N=1000$ , the normalized standard deviation of each parameter estimator is less than 0.06, which can better meet the actual communication use. Although the above conclusions are obtained under specific parameter conditions, they are not limited to this condition; similar conclusions are obtained under other parameter conditions.

Table 1.

Mean of the estimates.

$N$	$\overline{a}$	$\overline{\Sigma }$
100	[1.0125 + 0.9998i;1.0108 + 0.9925i]T	[1.1160,0.4802 + 0.5475i;0.4802–0.5475i,1.1272]
300	[0.9891 + 1.0105i; 0.9858 + 0.9882i]T	[1.0704,0.5185 + 0.5398i;0.5185–0.5398i,1.0655]
500	[0.9990 + 1.0011i; 0.9990 + 0.9998i]T	[1.0397,0.5130 + 0.5241i;0.5130–0.5241i,1.0393]
1000	[0.9996 + 1.0013i; 0.9998 + 1.0043i]T	[1.0268,0.5047 + 0.5111i;0.5047–0.5111i,1.0360]

Table 2.

Normalized sample standard deviation of estimated parameters.

N	da1	da2	dΣ (1,1)	dΣ (1,2)	dΣ (2,1)	dΣ (2,2)
100	0.0718	0.0688	0.2897	0.6071	0.6071	0.6029
300	0.0396	0.0391	0.1242	0.1516	0.1516	0.1150
500	0.0344	0.0342	0.0848	0.1072	0.1072	0.0947
1000	0.0219	0.0237	0.0333	0.0504	0.0504	0.0502

The bit error rate (BER) evaluation for the proposed blind receiver is an important metric, and it was tested through the simulation channel. In this simulation experiment, we set the parameter condition as $\mathit{a}=[1+i;1+i]$ , $\sum =[{\sigma }^2,{\sigma }^2(\rho +\rho i);{\sigma }^2(\rho -\rho i),{\rho }^2]$ . To evaluate the normalized BER performance of the proposed blind receiver, we compared the blind receiver with the optimal receiver that uses the maximum likelihood estimation method under known parameters. As shown in Table 3, the comparison between the blind receiver and the optimal receiver under $\alpha =1.6$ and $\rho =0.2$ conditions was performed 100 times for each set of simulation experiments, in which the normalized standard deviation, $d(\mathrm{BE}{R}_{\mathrm{opt}})$ , and $d(\mathrm{BE}{R}_{\mathrm{EM}})$ were defined as in equation (38).

Table 3.

BER performance of optimum receivers and blind receivers with $\alpha =1.6$ and $\rho =0.2$ .

	$\overline{\mathrm{BE}}{R}_{\mathrm{opt}}$	$\overline{\mathrm{BE}}{R}_{\mathrm{em}}$	$d(\mathrm{BE}{R}_{\mathrm{opt}})$	$d(\mathrm{BE}{R}_{\mathrm{opt}})$
${\sigma }^2=0.1$	0.0034	0.0035	0.0284	0.0350
${\sigma }^2=0.2$	0.0060	0.0062	0.0286	0.0290
${\sigma }^2=0.3$	0.0088	0.0088	0.0289	0.0351
${\sigma }^2=0.4$	0.0112	0.0113	0.0125	0.0151
${\sigma }^2=0.5$	0.0138	0.0142	0.0212	0.0324
${\sigma }^2=0.6$	0.0163	0.0164	0.0337	0.0352
${\sigma }^2=0.7$	0.0189	0.0189	0.0199	0.0216
${\sigma }^2=0.8$	0.0220	0.0223	0.0189	0.0221
${\sigma }^2=0.9$	0.0248	0.0249	0.0252	0.0276
${\sigma }^2=1.0$	0.0275	0.0279	0.0260	0.0312

Table 3 lists 10 sets of simulation results for ${\sigma }^2=0.1-1$ .

It can be seen from Table 3 that the BER performance for the proposed blind receiver very closely approximates the BER performance of the optimal receiver. The normalized standard deviation is also small, which is mainly due to the robustness of the EM algorithm. Figure 3 shows the comparison of the performance for the optimal receiver and the blind receiver at ${\sigma }^2=0.2$ condition. Each simulation experiment was repeated 500 times, and the results were averaged. As the value of $\alpha$ is generally between 1.6 and 2.0 in practice, the BER performance curves of the two receivers are conducted at $\alpha =\{1.6;1.7;1.8;1.9\}$ .

Figure 3.

Performance of optimum receivers and blind receivers with $\rho =0.2$

As shown in Figure 3, the performance of the blind receiver approaches that of the optimal receiver under the different values of $\alpha$ . In addition, the performance of the optimal receiver will decrease as $\alpha$ decreases. The reason may be explained when the value of $\alpha$ is smaller, the more pronounced the non-Gaussian noise channel is and the greater the probability of high-amplitude pulses appearing is. That is, when the value of $\alpha$ is smaller, the channel quality is worse, resulting in poorer signal reception quality and higher BER.

Conclusions

Due to the possible correlation between the sampling noise vector and the noise vector during the transmission of symbols, we propose a blind receiver based on a multidimensional stable noise distribution channel. The proposed blind receiver addresses the problem of signal detection under a correlated noise channel. An EM algorithm is derived for a multidimensional correlated noise channel, in particular, for blind signal detection in uncorrelated noise channels. Simulation results show that the proposed blind detector has fast iterative convergence with robustness, and the BER performance of the proposed blind receiver can approximate the optimal receiver’s BER performance under given parameters.