1 Introduction

Since the Earth and Mars are far away from each other, there is an attenuation of signal and exists a delay of time for communication between them, both of which destroy the signal quality in receivers [1, 2]. On the one hand, the time of essential data transmission is limited in Mars exploration system, and on the other hand, the relative position and condition of orbiters vary quickly all the time, which makes it difficult to obtain the accurate and significant information for Mars exploration [3, 4]. Therefore, it is obvious that the information of signal-to-noise ratio (SNR) promotes the system to adjust the signal transmission rate automatically [5]. The existing SNR estimation methods can be classified into two groups, namely data-aided (DA) and non-data-aided (NDA) [6, 7] estimators, respectively, which apply to the cases that the transmitted signals are known to the receivers for DA estimators and the receivers of the NDA estimators have no information about the transmitted symbols, respectively [8]. Although DA estimation algorithms outperform NDA-based methods in accuracy and effectiveness, it is at the expense of lower transmission efficiency, which is extremely worthwhile for deep space exploration [9]. Therefore, NDA-based estimators are considered for the deep space scene. Recently, many new researches are also investigated [10–16].

Generally, NDA-based SNR estimators include signal-to-variation ratio (SVR) [17], second- and fourth-order moments (M2M4) [18], squared signal-to-noise variance (SNV), and subspace-based method [1]. SVR estimator, which can be applied to not only fading channels but also additive white Gaussian noise (AWGN) channels, is developed to monitor and evaluate the channel quality based on moment operations. It is noticeable that SVR estimator is not applicable to other forms of digital modulated signals except M-ary phase-shift keying (PSK) modulated ones. M2M4 estimator was first proposed to estimate the strength of carrier and noise in real channels, and we further extended it to complex domain. Due to its independence of the transmitted symbols, namely, it only requires the estimation of the second and fourth order related to the received symbols, M2M4 is also one of the NDA estimators. SNV estimator utilizes the first and second moment of the received signals, which is the output of the matched filter (MF). Finally, the subspace method is based on the singular value decomposition (SVD) theory.

This paper aims to develop a NDA SNR estimation method based on maximum likelihood (ML) method to improve the system performance of deep space model; at the same time, it can achieve a higher accuracy than other estimators, such as SVR, SNV, M2M4, and subspace-based method. The Cramer-Rao lower bound (CRLB) of this proposed NDA estimation method is also derived to verify the effectiveness of the proposed ML-based estimator.

The rest of this paper is organized by the following structure. In Section 2, we discuss the system model. The maximum likelihood method based on data-aided is explained in Section 3. Section 4 demonstrates the NDA SNR estimation algorithm through maximum likelihood method. Section 5 shows the experimental figures. Section 6 concludes the contributions made by this paper.

2 System model

We aim to find the best SNR estimator while consuming the least energy. Table 1 is the notation of variables for this paper. Generally, the statistical properties are averaged to generate the SNR estimate through a large amount of symbols. The discrete and complex binary phase-shift keying (BPSK) [19] signal constellations are adopted to generate the baseband-equivalent and band-limited transmit symbols in real AWGN channels as illustrated in Fig. 1. The length of the transmitted symbols, which is upsampled to L = 16 samples in each symbol, is denoted as N. The root raised-cosine (RRC) filter, whose rolloff R equals to 0.5, and number of tap coefficients Q equals to 65, is assumed [20]. The binary source symbols can be calculated by [21, 22]:

$$\begin{array}{@{}rcl@{}} s_{n} = ae^{j\alpha_{n}}, \end{array} $$
(1)
Fig. 1
figure 1

The system model

Full size image
Table 1 Notations of variables
Full size table

where a represents the signal amplitude whose probability is equal by obtaining values from {−A,A} and αn represents one of the two evenly spaced phases on a unit circle with \( n\in \left \lbrace 1,2, \dots, N\right \rbrace \). The upsampled information sequence can be shown as:

$$\begin{array}{@{}rcl@{}} c_{k} = \sum_{n}s_{n}\delta_{k,nN}, \end{array} $$
(2)

where δij denotes the Kronecker delta. The signal after being sampled and pulse-shaped can be given as:

$$\begin{array}{@{}rcl@{}} d_{k} = \sum_{n}s_{n}h_{k-nN}, \end{array} $$
(3)

where hk is the tap coefficients of the RRC filter with \( k\in \left \lbrace -\left (R-1\right)/2,\dots,-1,0,1,\dots, \left (R-1\right)/2 \right \rbrace \) and hk equals to zero for |k|>(R−1)/2. The received signal can be presented as:

$$\begin{array}{@{}rcl@{}} r_{k} = d_{k} + w_{k}, \end{array} $$
(4)

where wk denotes the complex and sampled AWGN with zero mean and variance of σ2. The output of the MF is expressed as:

$$\begin{array}{@{}rcl@{}} y_{k} = r_{k}\otimes h_{-k}^{\ast} = \sum_{l}h_{l}d_{k-l} + \sum_{l}h_{l}w_{k-l}, \end{array} $$
(5)

where ⊗ and ∗ denote the discrete convolution and complex conjugation, respectively, and \( h_{k} = h_{-k}^{\ast } \) can be attributed to the assumption that the impulse response of RRC is real and even symmetric. Finally, the downsampled symbols at the receiver can be expressed as:

$$\begin{array}{@{}rcl@{}} y_{n} = y_{k}\mid_{k=nN} = s_{n}f_{0} + z_{n}, \end{array} $$
(6)

where f0 denotes the maximum impulse response of the full raised-cosine, the samples of which can be written as:

$$\begin{array}{@{}rcl@{}} f_{k} = h_{-k}^{\ast}\otimes h_{k} = \sum_{l}h_{k-l}h_{l}, \end{array} $$
(7)

and the downsampled and filtered samples of noise zn can be expressed as:

$$\begin{array}{@{}rcl@{}} z_{n} = z_{k}\mid_{k=nN} = \sum_{l}h_{l}w_{k-l}\mid_{k=nN}, \end{array} $$
(8)

Subsequently, the SNR is deduced as:

$$\begin{array}{@{}rcl@{}} \rho = \frac{E\left\lbrace \left| s_{n}f_{0}\right|^{2}\right\rbrace }{var\left\lbrace z_{n}\right\rbrace }, \end{array} $$
(9)

where E{·} represents the expectation and var{·} represents the variance. The SNR can be independent of the channel by normalizing the corresponding square of tap coefficients of the RRC; subsequently, the SNR can be solely related to A and σ2, namely:

$$\begin{array}{@{}rcl@{}} \rho = \frac{A^{2}}{2\sigma^{2}}. \end{array} $$
(10)

3 Maximum likelihood estimation method based on NDA

In this section, we propose an NDA estimation method [23, 24] based on ML algorithm. The probability density function (PDF) of yn can be expressed as:

$$\begin{array}{@{}rcl@{}} P\left(y_{n}\right) = \frac{1}{2}\left\lbrace P_{+}\left(y_{n}\right) + P_{-}\left(y_{n}\right) \right\rbrace = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{y_{n}^{2}+A^{2}}{2\sigma^{2}}}\cosh\left(\frac{y_{n}A}{\sigma^{2}}\right), \end{array} $$
(11)

where \( P_{+}\left (y_{n}\right) = \frac {1}{\sqrt {2\pi }\sigma }e^{-\frac {\left (y_{n}-A\right)^{2}}{2\sigma ^{2}}}, P_{-}\left (y_{n}\right) = \frac {1}{\sqrt {2\pi }\sigma }e^{-\frac {\left (y_{n}+A\right)^{2}}{2\sigma ^{2}}} \), and \( \cosh \left (x\right) = \frac {1}{2}\left (e^{x} + e^{-x}\right) \).

The joint PDF of the signal vector at the receiver \( \left [ y_{1},y_{2},\dots,y_{n}\right ] \) is hence given by (due to independence):

$$\begin{array}{@{}rcl@{}} P_{N}\left(y_{1},y_{2},\dots,y_{n}\right) = \prod_{n=1}^{N}P\left(y_{n}\right) = \prod_{n=1}^{N}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{y_{n}^{2}+A^{2}}{2\sigma^{2}}}\cosh\left(\frac{y_{n}A}{\sigma^{2}}\right). \end{array} $$
(12)

The log-likelihood function is described as:

$$\begin{array}{@{}rcl@{}} L_{N}\left(y_{1},y_{2},\dots,y_{N}\mid A,\sigma^{2}\right) =&-&N\ln\left(\sqrt{2\pi}\sigma\right)- \frac{1}{2\sigma^{2}}\sum_{n=1}^{N}\left(y_{n}^{2} + A^{2}\right)\\ &+& \sum_{n=1}^{N}\ln\cosh\left(\frac{y_{n}A}{\sigma^{2}}\right). \end{array} $$
(13)

Setting \( \frac {{\partial {L_{N}}\left ({{y_{1}},{y_{2}}, \ldots,{y_{N}}} \right)}}{{\partial A}} = 0 \), we implicitly obtain the ML estimate of A as the solution to the equation:

$$\begin{array}{@{}rcl@{}} A = \frac{1}{N}\sum_{n=1}^{N}y_{n}\tanh\left(\frac{y_{n}A}{\sigma^{2}}\right). \end{array} $$
(14)

We consider an iterative algorithm to explore the most suitable value of amplitude that satisfies (13). According to the signal vector at the receiver with N samples of yn, it is obvious to define the function:

$$\begin{array}{@{}rcl@{}} F\left(x\right) = x - \frac{1}{N}\sum_{n=1}^{N}y_{n}\tanh\left(\frac{y_{n}x}{\sigma^{2}}\right). \end{array} $$
(15)

By solving the equation above, we can obtain the amplitude estimation of the maximum likelihood method, namely, F(x)=0 at \( x = \hat {A} \), where \( \hat {A} \) represents the estimate of A. We can determine the value of \( \hat {A} \) by the following steps:

Step 1: The received vector is normalized as \(\frac {1}{N}\sum _{n=1}^{N}y_{n}^{2} = 1\). Find the maximum and minimum values of Amin and Amax. At the same time, determine the total number of iteration I. Let A1=Amin and A2=Amax. Initialize i = 0.

Step 2: Compute \( A_{m} = \frac {A_{1}+A_{2}}{2} \), and therefore, \( \sigma _{m}^{2} = 1-A_{m}^{2} \).

Step 3: Compute \( F\left (A_{m}\right) = A_{m} - \frac {1}{N}\sum _{n=1}^{N}y_{n}\tanh \left (\frac {y_{n}A_{m}}{\sigma ^{2}}\right) \).

Step 4: If F(Am)>0, then let A2=Am; otherwise, let A1=Am.

Step 5: Let i = i+1, if i = I, output \( A_{m} = \frac {A_{1}+A_{2}}{2} \) as to the ML estimate of the amplitude and \(\frac {A_{m}^{2}}{1-A_{m}^{2}}\) as the estimated SNR. Otherwise, go back to step 2.

4 Performance comparisons

4.1 Cramer-Rao lower bound (CRLB)

We consider CRLB [25, 26] to evaluate whether the estimators work effectively or not. The definition of the SNR has been given in (10). We intend to estimate ρ by N observed samples of yn. The estimation task involves two parameters. The estimated parameter vector is denoted as:

$$\begin{array}{@{}rcl@{}} \theta = \left[ A \quad \sigma^{2}\right]^{T}. \end{array} $$
(16)

Since the estimation of the SNR is generally expressed in decibel scale, the following form of CRLB is adopted:

$$\begin{array}{@{}rcl@{}} g\left(\theta\right) = 10\log_{10}\left(\frac{A^{2}}{\sigma^{2}}\right), \end{array} $$
(17)
$$\begin{array}{@{}rcl@{}} \text{CRLB}\left(\rho\right) = \frac{\partial g\left(\theta \right)}{\partial \theta }K^{-1}\left(\theta\right)\frac{\partial g{{\left(\theta \right)}^{T}}}{\partial \theta }, \end{array} $$
(18)

where K(θ) is the Fisher information matrix, which can be shown as:

$$\begin{array}{@{}rcl@{}} K\left(\theta\right) = \left[ {\begin{array}{*{20}{c}} { - E\left({\frac{{\partial^{2}}\ln P\left({x;\theta} \right)}{\partial {A^{2}}}} \right)}&{ - E\left({\frac{{\partial^{2}}\ln P\left({x;\theta} \right)}{\partial A\partial {\sigma^{2}}}} \right)}\\ { - E\left({\frac{{\partial^{2}}\ln P\left({x;\theta} \right)}{\partial {\sigma^{2}}\partial A}} \right)}&{ - E\left({\frac{{\partial^{2}}\ln P\left({x;\theta} \right)}{\partial {\sigma^{2^{2}}}}} \right)} \end{array}} \right]. \end{array} $$
(19)

From (17), \( \frac {\partial g\left (\theta \right)}{\partial \theta } \) is determined as:

$$\begin{array}{@{}rcl@{}} \frac{\partial g\left(\theta \right)}{\partial \theta }{\mathrm{ = }}\left[ {\begin{array}{*{20}{c}} {\frac{20}{\ln \left({10} \right)A}}&{\frac{ - 10}{\ln \left({10} \right){\sigma^{2}}}} \end{array}} \right]. \end{array} $$
(20)

According to the N observed symbols, the Fisher information matrix is described as:

$$\begin{array}{@{}rcl@{}} K_{\mathrm{BPSK-NDA}}\left(\theta\right) = \frac{N}{\sigma^{4}} \left[ {\begin{array}{*{20}{c}} {{\sigma^{2}} - {\sigma^{2}}J\left(\rho \right)}&{AJ\left(\rho \right)}\\ {AJ\left(\rho \right)}&{1 - \frac{A^{2}}{\sigma^{2}}J\left(\rho \right)} \end{array}} \right], \end{array} $$
(21)

where \( J\left (\rho \right) = \frac {e^{- \rho }}{\sqrt {2\pi } }\int \limits _{- \infty }^{+ \infty } {\frac {{u^{2}}{e^{- \frac {u^{2}}{2}}}}{\cosh \left ({u\sqrt {2\rho } } \right)}} du \). By substituting (21) to (17), the CRLB of the BPSK signals can be expressed as:

$$\begin{array}{@{}rcl@{}} \mathrm{CRLB_{BPSK-NDA}}\left(\rho\right) = \frac{{100\left({\frac{2}{\rho} - J\left(\rho \right) + 1} \right)}}{{N{{\left({\ln \left({10} \right)} \right)}^{2}}\left({1 - J\left(\rho \right) - 2\rho J\left(\rho \right)} \right)}}{\left({dB} \right)^{2}}. \end{array} $$
(22)

4.2 Normalized mean square error (NMSE)

For the above iterative SNR estimator, we provide the required experimental results. To appropriately benchmark our proposed method, the results of SVR, M2M4, SNV, and subspace-based methods are also included. According to the system model given in Section 2, we obtain the NMSE of the aforementioned SNR estimators. The NMSE represents the deviation between estimated and true values and can be calculated as follows:

$$\begin{array}{@{}rcl@{}} \text{NMSE}\left\lbrace \hat{\rho}\right\rbrace = \frac{E\left\{ {{\left({\rho - \hat \rho} \right)}^{2}} \right\}}{\rho^{2}}, \end{array} $$
(23)

where \( \hat {\rho } \) describes the estimated SNR while ρ represents the true SNR value. The NMSE of each estimator is calculated by using the Monte Carlo method as follows:

$$\begin{array}{@{}rcl@{}} \text{NMSE}\left\lbrace \hat{\rho}\right\rbrace = \frac{1}{N}\sum\limits_{n = 1}^{N} {\frac{{{{\left({\rho - \hat \rho} \right)}^{2}}}}{{{\rho^{2}}}}}. \end{array} $$
(24)

The number of symbols N should be large enough to promise the accuracy and objectivity of the experimental results. In terms of the ML estimator, the NMSE under different simulation numbers is also shown to further compare the performance of the proposed estimator.

5 Simulation results

In this section, the Monte Carlo simulation results of the NMSE performances of the aforementioned estimators are presented. The CRLB performance, the ML estimator under several different simulation numbers, and the differences among perfect SNR, SNR under subspace-based method, and SNR under ML-based method are also displayed.

Figure 2 shows the NMSE performance among SNV, SVR, M2M4, subspace method, and ML-based estimator with the length of symbols N = 1000. From the simulation result, we can make a conclusion that SNV-based method performs worse than all the other estimators mentioned above. Subspace-based method outperforms SNV, SVR, and M2M4 in a reasonable SNR region. However, when compared to the proposed ML-based estimator, it suffers an SNR deficit under the same NMSE, which verifies the efficiency of the proposed ML-based method.

Fig. 2
figure 2

The NMSE comparison of several SNR estimators. The NMSE of SVR, SNV, M2M4, subspace method, and ML under N = 1000

Full size image

Figure 3 describes the CRLB of the proposed ML-based estimator under two cases, namely, N = 300 and N = 600, respectively. From the simulation result, we can find that the CRLB with N = 300 is higher than that of N = 600. The reason of this phenomenon can be easily understood by the essence that the larger simulation symbols lead to a higher estimation accuracy.

Fig. 3
figure 3

The CRLB of the proposed algorithm. ML under N = 300, ML under N = 600, CRLB under N = 300, and CRLB under N = 600

Full size image

Figure 4 demonstrates the NMSE performance under different lengths of simulation symbols, which are N = 100, N = 200, and N = 500, respectively. It is obvious that with the increase of simulation symbols N, the NMSE performance becomes better, which is consistent with the theoretical analysis. This result is useful for engineering to choose the appropriate length of symbols.

Fig. 4
figure 4

The NMSE under different length of simulation symbols. N = 100, N = 200, and N = 500, respectively

Full size image

Figure 5 shows the differences among the perfect SNR, SNR under subspace-based method, and SNR under ML-based method under N = 600. We can obtain from the curves that compared to the subspace-based method, the ML-based method achieves a higher estimation accuracy and the gap across the perfect SNR is smaller. Therefore, the proposed ML-based method performs better than all the other estimators mentioned above.

Fig. 5
figure 5

The difference comparison among the subspace method-based SNR, ML-based SNR, and perfect SNR under N = 600. Perfect SNR (N = 600), subspace-based method (N = 600), and ML-based method (N = 600)

Full size image

6 Conclusions

We proposed a novel NDA-based ML estimator, which is based on an iterative algorithm and achieves a higher SNR transmission accuracy compared to several traditional SNR estimators. The simulation results showed that the proposed new ML estimation method achieves a superior NMSE performance compared to the existing ones, which is significant for the researches of Mars exploration. The CRLB of the proposed ML-based method was also derived. Furthermore, the NMSE of the ML-based method under different simulation symbols was compared to further verify the effectiveness of the proposed algorithm.