Phased Fractional Low-Order Moment-Based Doppler Shift Estimation in the Presence of Interference Signals and Impulsive Noise

Abstract

Doppler shift estimation continues to be a critical challenge of utmost significance in both theoretical research and practical engineering applications. Many innovators have crafted solutions specific to this issue, with notable contributions across various signals and scenarios. Given that cyclostationary signals are prevalent in both artificial and natural phenomena, we propose a novel framework based on the phased fractional lower-order moment (PFLOM) for estimating Doppler shift in mixed cyclostationary signals. During the estimation process, a more realistic impulse noise model is examined in contrast to the ideal Gaussian noise typically assumed in conventional methods. This approach is meticulously derived through a series of detailed steps in line with cyclostationary signal processing and PFLOM principles. Furthermore, an extensive simulation has been conducted to validate the efficacy and robustness of our proposed method. It is anticipated that the concept and method presented here could be applied more broadly due to its solid theoretical underpinnings.

1. Introduction

Cyclostationary signals (processes), both in communication and mechanical systems, present a captivating subset of non-stationary signals characterized by periodic variations in their statistical properties, setting them apart from traditional stationary or non-stationary signals [1,2,3,4,5]. This periodicity extends beyond artificial signals to encompass natural phenomena like astronomical, hydrological, and oceanographic signals, influenced by Earth’s rotational and orbital cycles. Moreover, even in biological contexts such as electrocardiograms and electroencephalograms, periodic phenomena arise from the rhythmic patterns of the heartbeat and blood circulation. Given the broad spectrum of cyclostationary signals and their relevance across diverse domains, the relative research holds significant importance [6,7,8,9,10,11].

Doppler shift estimation holds immense significance across diverse fields, spanning communications, radar, medical imaging, and astronomy [12,13,14,15,16,17]. In communications, the accurate estimation of Doppler shifts ensures the reliability and efficiency of systems, particularly in scenarios involving high-speed mobility, enabling the dynamic adaptation of modulation schemes and mitigation of fading and interference. Radar systems rely on Doppler shift estimation for target detection, tracking, and velocity measurement, which are crucial for applications such as air traffic control, weather monitoring, and military defense. In medical imaging, Doppler ultrasound techniques utilize Doppler shift estimation to visualize blood flow dynamics and detect abnormalities, aiding diagnoses in cardiology, vascular medicine, and obstetrics. In astronomy, Doppler shift estimation allows astronomers to study the motion and properties of celestial objects, facilitating the discovery and characterization of explanatory systems and advancing our understanding of the universe. Overall, Doppler shift estimation serves as a fundamental tool for analyzing relative motion, providing insights into dynamic systems and phenomena essential for technological advancement and scientific discovery across a broad spectrum of disciplines.

To achieve Doppler shift estimation for mixed cyclostationary signals in impulsive noise is a challenging task. On the one hand, cyclostationary signals interfere severely with each other when their amplitudes are equal or similar. On the other hand, the traditional second- and higher-order-based methods fail when the background noise does not obey Gaussian distribution [18,19]. Gaussian noise is an ideal model designed to simplify the algorithm derivation since its characteristic exponent is a fixed constant of 2 [20,21]. This paper is proposed to complete the above-mentioned tough task through a suitable method.

To address the first problem, the 3D cyclic spectrum based on cyclostationary signal processing technology is employed [22,23]. The spectrum is built in both the frequency domain and cyclic frequency domain. The frequency domain is used to estimate Doppler shift, while the cyclic frequency domain is used to estimate the cyclic frequency. Both domains are used together to separate the cyclostationary signals. To address the second problem, PFLOM technology is employed. Nowadays, there are several effective technologies that can handle impulsive noise and Gaussian noise, including PFLOM [24,25], correntropy-like methods [26,27,28], and hyperbolic function-based methods [29,30,31,32,33]. However, correntropy-like methods are passed considering the necessary phase information used in the frequency domain. Hyperbolic function-based methods are passed considering the computational complexity caused by the exponent functions.

Through serious thinking about the technologies of cyclostationary signal processing and PFLOM, a novel PFLOM-based method is proposed in this paper to complete robust Doppler shift estimation for mixed cyclostationary signals in impulsive noise. The fractional order p in PFLOM is adjusted in a suitable interval to realize robustness against the impulsive noise modeled by alpha-stable distribution. Alpha-stable distribution is an extension of Gaussian distribution in theory according to the central limit theorem (CLT) [34,35] and the generalized central limit theorem (GCLT) [36,37]. It has a stronger generalization in various scenarios, which corresponds to the wide range of cyclostationary signals. In this paper, the essential issue comes from the needs of various applications, and the whole procedure is processed based on the mathematical derivations step by step.

The remainder of this paper is organized as follows. In Section 2, cyclostationary signal processing and impulsive noise modeled by alpha-stable distribution are introduced briefly. In Section 3, the methodology including a mathematical model for the receiving end, a basic idea for the framework, and an explanation of choosing PFLOM is studied theoretically, where the mathematical derivations are provided step by step. In Section 4, the simulation including one single experiment and Monte Carlo experiments is carried out to demonstrate the superiority of the proposed method, and four remarks are given from the results of the simulation. In Section 5, the main contribution and related reflection of this work are concluded.

2. Related Works

2.1. Cyclostationary Signal Processing in Ideal and Real Scenarios

Given a non-stationary random signal $X\left(t\right)$, its autocorrelation on time t and delay $\tau$ is defined by [6,8]

$$R_X(t,\tau )\triangleq \mathrm{E}\left[X\left(t\right)X^{*}(t+\tau )\right]=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^N{x}_n\left(t\right)x_n^{*}(t+\tau ),$$

(1)

where $X\left(t\right)$ and $X(t+\tau )$ are both random variables on fixed time t. In fact, the correlation describes the relationship between the random variables and has had nothing to do with the time functions until now. $x_n\left(t\right)$ denotes the n-th sample function from $X\left(t\right)$. Performing the statistical averaging $\mathrm{E}[·]$ on a random signal $X\left(t\right)$ yields a time function $R_X(t,\tau )$. In order to emphasize the fixed times, there is also some literature that uses time $t_1$ and $t_2$ to replace t and $\tau$, and the autocorrelation is expressed as follows [6,8].

$$R_X(t_1,t_2)\triangleq \mathrm{E}\left[X\left(t_1\right)X^{*}\left(t_2\right)\right]=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^N{x}_n\left(t_1\right)x_n^{*}\left(t_2\right),$$

(2)

where $t_1=t$ and $t_2=t+\tau$.

In practice, it is impossible to obtain an infinite number of sample functions. Therefore, time averaging is commonly used in place of statistical averaging under strict ergodicity. Further, if $R_X(t,\tau )$ is either a period function or an almost period function regarding t, then we say that the non-stationary random signal $X\left(t\right)$ has a second-order cyclostationarity. And $R_X(t,\tau )$ can be expanded by Fourier series (FS) [6,8].

$$R_X(t,\tau )=\sum _{\xi }{V}_X(\xi ,\tau ){\mathrm{e}}^{\mathrm{j}2\pi \xi t}$$

(3)

where $V_X(\xi ,\tau )$ denotes a Fourier-series coefficient and is also named cyclic autocorrelation. $\xi$ denotes Fourier frequencies and is also named cyclic frequencies. Based on the theory of FS, $V_X(\xi ,\tau )$ can be computed by [6,8]

$$V_X(\xi ,\tau )={\displaystyle \frac{1}{T_0}}{\int }_{T_0}{R}_X(t,\tau ){\mathrm{e}}^{-\mathrm{j}2\pi \xi t}\mathrm{d}t$$

(4)

where $T_0$ denotes the fundamental period of $R_X(t,\tau )$. At last, applying Fourier transform (FT), not FS, on $V_X(\xi ,\tau )$ on $\tau$, the spectrum of cyclic autocorrelation, also named the cyclic spectrum, is obtained [6,8].

$$Z_X(\xi ,f)={\int }_{{\Lambda }_2}{V}_X(\xi ,\tau ){\mathrm{e}}^{-\mathrm{j}2\pi f\tau }\mathrm{d}\tau$$

(5)

where ${\Lambda }_2$ denotes the search range for the delay $\tau$. Since FT can work for both periodic and non-periodic signals, the process of FS can be replaced by FT considering the real situations of background noises in $X\left(t\right)$ and $R_X(t,\tau )$ almost being a period function. In digital signal processing, FT can be simply replaced by discrete Fourier transform (DFT) through sampling in time, frequency, and cyclic frequency domains.

2.2. Proper Model for Impulsive Noise

Gaussian noise is used as a common model since it is observed in all electronic equipment, and in many other aspects of nature.

For impulsive noise, it is often modeled by generalized Gaussian distribution, the Gaussian mixture model, and alpha-stable distribution in various studies. Compared with the generalized Gaussian distribution and Gaussian mixture model, the most significant feature of alpha-stable distribution is that it comes from the generalized central limit theorem, while the original Gaussian distribution comes from the central limit theorem. In practice, the noise obeying alpha-stable distribution (heavy-tailed distribution) is rare [38,39,40]. But considering that Gaussian distribution is a subcase of alpha-distribution and it constitutes a rich class, such as Cauchy distribution, Levy distribution, and so on, we prefer alpha-stable distribution to model impulsive noise.

Definition 1.

Given two independent and identically distributed (i.d.d.) variables $X_1,X_2\sim X$ and $Y=a{X}_1+b{X}_2,a,b>0$, if $Y\sim cX+d,c>0,d\in \mathbb{R}$, then X is alpha-stable distributed.

Alpha-stable distribution does not have a closed form of probability density function (PDF). Therefore, its characteristic function (CF) is commonly used in the study, which is defined by [20,21]

$${\phi }_X\left(t\right)=exp\left\{\mathrm{j}\delta t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\left[1-\mathrm{j}\beta \mathrm{sgn}\left(t\right)\varphi (t,\alpha )\right]\right\},$$

(6)

where t is a real number, $exp(·)$ denotes the exponential function, $\mathrm{sgn}(·)$ denotes the sign function, and

$$\varphi (t,\alpha )=\left\{\begin{array}{}\mathrm{(7a)}& \tan \left(\frac{\pi \alpha }{2}\right)\hfill & \mathrm{for}\alpha \ne 1\mathrm{(7b)}& -\frac{2}{\pi }\ln |t|\hfill & \mathrm{for}\alpha =1.\end{array}\right.$$

$\alpha \in (0,2]$ denotes the characteristic exponent describing the tails of alpha-stable distribution. The smaller $\alpha$ is, the heavier the tails are. $\beta \in [-1,1]$, $\gamma \in (0,\infty )$, and $\delta \in (-\infty ,\infty )$ denote skewness, scale, and location, respectively. Given Gaussian distribution $\mathcal{N}(\mu ,{\sigma }^2)$, $\gamma$ and $\delta$ are much similar to $\sigma$ and $\mu$, respectively. According to the CF, the PDF of alpha-stable distribution can be calculated through Fourier transform as follows [20,21]:

$$f_X\left(x\right)={\int }_{\infty }{\phi }_X\left(t\right){\mathrm{e}}^{-\mathrm{j}tx}\mathrm{d}t={\int }_{\infty }exp\left\{\mathrm{j}(\delta -x)t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\left[1-\mathrm{j}\beta \mathrm{sgn}\left(t\right)\varphi (t,\alpha )\right]\right\}\mathrm{d}t.$$

(8)

From (6) and (8), both the PDF and CF of alpha-stable distribution are complicated. Without a loss of generality, $\beta$ is sometimes set to zero to simplify the two. In this case, alpha-stable distribution degenerates into symmetric $\alpha$-stable (S$\alpha$S) distribution. The PDF and CF of S$\alpha$S distribution are given by [20,21]

$$\begin{array}{cc}\hfill {\tilde{\phi }}_X\left(t\right)& =exp\left\{\mathrm{j}\delta t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\right\}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\tilde{f}}_X\left(x\right)& ={\int }_{\infty }exp\left\{\mathrm{j}(\delta -x)t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\right\}\mathrm{d}t.\hfill \end{array}$$

(10)

3. Methodology

3.1. Mathematical Model for the Receiving End

The scenario presented in our work is the elimination of time delay between two receivers in uncooperative communication. The information about the transmitter is unknown, which is a bit different from the scenario for the carrier-frequency offset (CFO) estimation in cooperative communication. The mathematical model of the received signals is given by

$$r\left(t\right)=\sum _{i=1}^N{a}_i{s}_i(t-d_i){\mathrm{e}}^{\mathrm{j}2\pi D_{i}t}+n\left(t\right)$$

(11)

where $a_i$, $s_i\left(t\right)$, $d_i$, $D_i$, and $n\left(t\right)$ denote the decreasing coefficient, sending signal, time delay, Doppler shift, and additive noise, respectively. It can be observed that when $N\ge 2$, there are multiple sources that interfere with one another. If $s_i\left(t\right)$ is designated as the signal of interest (SOI), then $s_j\left(t\right),j\ne i$ are classified as interference signals. When $s_i\left(t\right)$ is a cyclostationary signal, the cyclic frequency inside and the corresponding Doppler shift can be estimated through the second-order statistics. Commonly, the second-order statistics of $s_i\left(t\right)$ are obtained by the second-order statistics of $r\left(t\right)$.

$$\begin{array}{cc}\hfill & R_r(t,\tau )=\mathrm{E}\left[r\left(t\right)r^{*}(t+\tau )\right]\hfill \\ \hfill & =\mathrm{E}\left\{\left[\sum _{i=1}^N{a}_i{s}_i(t-d_i){\mathrm{e}}^{\mathrm{j}2\pi D_{i}t}+n\left(t\right)\right]\times \left[\sum _{j=1}^N{a}_j{s}_j^{*}(t+\tau -d_j){\mathrm{e}}^{-\mathrm{j}2\pi D_j(t+\tau )}+n^{*}(t+\tau )\right]\right\}\hfill \\ \hfill & =\mathrm{E}\left[\sum _{i=1}^N\sum _{j=1}^N{s}_i(t-d_i)s_j^{*}(t+\tau -d_j){\mathrm{e}}^{\mathrm{j}2\pi (D_{i}t-D_{j}t-D_j\tau )}\right]+E\left[n\left(t\right)n^{*}(t+\tau )\right]\hfill \\ \hfill & +\mathrm{E}\left[\sum _{i=1}^N{s}_i(t-d_i){\mathrm{e}}^{\mathrm{j}2\pi D_{i}t}{n}^{*}(t+\tau )\right]+\mathrm{E}\left[\sum _{j=1}^N{s}_j^{*}(t+\tau -d_j){\mathrm{e}}^{-\mathrm{j}2\pi D_j(t+\tau )}n\left(t\right)\right]\hfill \end{array}$$

(12)

Commonly, the signal and noise are usually assumed to be independent of each other. We have

$$\begin{array}{cc}\hfill & \forall i=1,2,\cdots ,N,\mathrm{E}\left[s_i(t-d_i){\mathrm{e}}^{\mathrm{j}2\pi D_{i}t}{n}^{*}(t+\tau )\right]=0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & \forall j=1,2,\cdots ,N,\mathrm{E}\left[s_j^{*}(t+\tau -d_j){\mathrm{e}}^{-\mathrm{j}2\pi D_j(t+\tau )}n\left(t\right)\right]=0.\hfill \end{array}$$

(14)

Further, the signals are independent of each other. We have

$$\begin{array}{cc}\hfill & \forall i\ne j,\mathrm{E}\left[s_i(t-d_i)s_j^{*}(t+\tau -d_j){\mathrm{e}}^{\mathrm{j}2\pi (D_{i}t-D_{j}t-D_j\tau )}\right]=0.\hfill \end{array}$$

(15)

Based on (13)–(15), $R_r(t,\tau )$ can be simplified as follows:

$$\begin{array}{c}\hfill R_r(t,\tau )=\mathrm{E}\left[\sum _{i=1}^N{a}_i^2{s}_i(t-d_i)s_i^{*}(t-d_i+\tau ){\mathrm{e}}^{-\mathrm{j}2\pi D_i\tau }\right]+\mathrm{E}\left[n\left(t\right)n^{*}(t+\tau )\right].\end{array}$$

(16)

In theory, the second-order statistics of $n\left(t\right)$ do not exist due to it obeying alpha-stable distribution. In practice, the second-order statistics of $n\left(t\right)$ exist because of the amplitude limiter at the receiving end. However, the periodic characterization of $R_r(t,\tau )$ is disrupted more or less, and the methodology based on it must be severely degraded. The reason is that the periodic characterization of $R_r(t,\tau )$ is disrupted due to the term of $\mathbf{E}\left[n\left(t\right)n^{*}(t+\tau )\right]$ when the additive noise $n\left(t\right)$ obeys alpha-stable distribution. More precisely, the second-order statistics of the random variable obeying alpha-stable distribution do not exist. In this case, we must use an efficient alternative to replace the traditional $R_r(t,\tau )$ and complete the estimation in a workable framework.

3.2. Basic Idea for the Framework

According to the above discussion, we need an alternate form to replace $R_r(t,\tau )$ to guarantee the existence of $\mathrm{E}\left[n\left(t\right)n^{*}(t+\tau )\right]$. $R_r(t,\tau )$ is computed by $r\left(t\right)$, and we employ the fractional lower-order statistics (FLOS) to modify $r\left(t\right)$. Among various FLOS, PFLOM is relatively good in this scenario. PFLOM is based on the fractional lower-order operator, and the fractional lower-order operator of the complex number x is defined by [18]

$$x^{\langle p\rangle }={\left|x\right|}^{p-1}x,x^{-\langle p\rangle }={\left|x\right|}^{p-1}{x}^{*},p\in [0,1].$$

(17)

Here, p denotes the p-th order. For the real number, (17) is simplified as [18]

$$x^{\langle p\rangle }=x^{-\langle p\rangle }={\left|x\right|}^p\mathrm{sgn}\left(x\right),p\in [0,1].$$

(18)

In general, impulsive noise can be handled when p tends to 0; Gaussian noise can be handled when p tends to 1. Especially, when p equals 0 and 1, $x^{\langle p\rangle }$ represents the phase and the number itself, respectively. For example, when $x=25$ and $p=0.5$, we have $x^{\langle p\rangle }=5$, and we can see that the amplitude of one sampling point of the receiving signal is compressed. Furthermore, the modified correlation, i.e., the PFLOM of $r\left(t\right)$ and $r(t+\tau )$, is given by

$$\begin{array}{cc}\hfill & R_r^{\langle p\rangle }(t,\tau )\hfill \\ \hfill & =\mathrm{E}\left[r^{\langle p\rangle }\left(t\right)r^{-\langle p\rangle }(t+\tau )\right]\hfill \\ \hfill & =\mathrm{E}\left[\sum _{i=1}^N{a}_i^{2p}{s}_i^{\langle p\rangle }(t-d_i)s_i^{-\langle p\rangle }(t-d_i+\tau ){\mathrm{e}}^{-\mathrm{j}2\pi D_i\tau }\right]\hfill \\ \hfill & =\mathrm{E}\left[{\mathrm{e}}^{-\mathrm{j}2\pi D_i\tau }\sum _{i=1}^N{a}_i^{2p}|s_i(t-d_i)s_i(t-d_i+\tau ){|}^{p-1}{s}_i(t-d_i)s_i^{*}(t-d_i+\tau )\right]\hfill \\ \hfill & +\mathrm{E}\left[n^{\langle p\rangle }\left(t\right)n^{-\langle p\rangle }(t+\tau )\right]\hfill \end{array}$$

(19)

From (19), it can be seen that the operation of ${(·)}^{\langle p\rangle }$ suppresses the amplitude and maintains the phase meanwhile. Further, as $f\left(x\right)={\left|x\right|}^p$ is a monotone function of the amplitude with respect to p, we have the following characteristic of the periodicity:

$$R_r(t,\tau )=R_r(t+T_1,\tau )\Rightarrow R_r^{\langle p\rangle }(t,\tau )=R_r^{\langle p\rangle }(t+T_1,\tau )$$

(20)

Applying FT on $R_r^{\langle p\rangle }(t,\tau )$ regarding t, we can obtain $V_r^{\langle p\rangle }(\xi ,\tau )$.

$$V_r^{\langle p\rangle }(\xi ,\tau )={\int }_{{\Lambda }_1}{R}_r^{\langle p\rangle }(t,\tau ){\mathrm{e}}^{-2\pi \xi t}\mathrm{d}t,$$

(21)

where ${\Lambda }_1$ denotes the search range for the time t. Applying FT on on $V_r^{\langle p\rangle }(\xi ,\tau )$ regarding $\tau$, we can obtain $Z_r^{\langle p\rangle }(\xi ,f)$, which can be used to jointly estimate the cyclic frequency ${\xi }_i$ and Doppler shift $D_i$ of $s_i\left(t\right)$ through peak finding.

$$Z_r^{\langle p\rangle }(\xi ,f)={\int }_{{\Lambda }_2}{V}_r^{\langle p\rangle }(\xi ,\tau ){\mathrm{e}}^{-2\pi f\tau }\mathrm{d}\tau ,$$

(22)

Based on $|Z_r^{\langle p\rangle }(\xi ,f)|$, the cyclic frequency ${\xi }_i$ and Doppler shift $D_i$ are jointly estimated by peak finding in the 3D spectrum $|Z_r^{\langle p\rangle }(\xi ,f)|$.

In our scenario, there are two receivers, and the received signals are $r\left(t\right)$ and $r(t+\tau )$. The time delay $\tau$ comes from the distance between the two receivers. When the number of statistically independent cyclostationary signals is one ($N=1$), there is only one peak in the positive domain of $|Z_r^{\langle p\rangle }(\xi ,f)|$. When the number of statistically independent cyclostationary signals is three ($N=3$), there are two peaks in the positive domain of $|Z_r^{\langle p\rangle }(\xi ,f)|$. Of course, with the increasing number of statistically independent cyclostationary signals, the cost of the two-dimensional search is expensive.

3.3. Explanation of Choosing PFLOM

Besides PFLOM, there are some methods that can effectively deal with the additive noises obeying alpha-stable distribution, such as correntropy-based methods and hyperbolic function-based methods. The correntropy of $r\left(t\right)$ and $r(t+\tau )$ is defined by [41]

$$C_r(t,\tau )=\mathrm{E}\left[K_{\sigma }\left(r\left(t\right),r(t+\tau )\right)\right]={\displaystyle \frac{1}{\sqrt{2\pi \sigma }}}\mathrm{E}\left[exp\left({\displaystyle \frac{{-\parallel r\left(t\right)-r(t+\tau )\parallel }_2^2}{2{\sigma }^2}}\right)\right],$$

(23)

where $K_{\sigma }(·,·)$ and $\sigma$ denote the Gaussian kernel and kernel size, respectively. No matter what the values of $n\left(t\right)$ and $n(t+\tau )$ are, $C_r(t,\tau )$ is bounded such that

$$\forall n\left(t\right),n(t+\tau )\in \mathbb{C},0\le C_r(t,\tau )\le {\displaystyle \frac{1}{\sqrt{2\pi \sigma }}}.$$

(24)

One problem of correntropy-based methods is that their results are real numbers without the frequency and phase information. Therefore, correntropy-based methods are hardly used for Doppler shift estimation. Besides, various kernels in correntropy-based methods introduce harmonics of higher orders. For instance, the exponent function in the Gaussian kernel can be expanded by

$$\begin{array}{cc}\hfill exp\left({\displaystyle \frac{{-\parallel r\left(t\right)-r(t+\tau )\parallel }_2^2}{2{\sigma }^2}}\right)& =exp\left({\displaystyle \frac{-{\parallel r\left(t\right)\parallel }_2^2-{\parallel r(t+\tau )\parallel }_2^2}{2{\sigma }^2}}\right)exp\left({\displaystyle \frac{2\langle r\left(t\right),r(t+\tau )\rangle }{2{\sigma }^2}}\right)\hfill \\ \hfill & =exp\left({\displaystyle \frac{-{\parallel r\left(t\right)\parallel }_2^2-{\parallel r(t+\tau )\parallel }_2^2}{2{\sigma }^2}}\right)\sum _{n=0}^{\infty }{\displaystyle \frac{{\langle r\left(t\right),r(t+\tau )\rangle }^n}{{\sigma }^{2n}n!}}.\hfill \end{array}$$

(25)

$\langle r\left(t\right),r(t+\tau )\rangle$ denotes the inner-product of two complex numbers $r\left(t\right)$ and $r(t+\tau )$. According to the analysis above, the harmonics of $s_i\left(t\right)$ and $s_j\left(t\right)$ operated by the kernel trick will affect each other. Thus, none or very little of an introduction of harmonics is important for multi-signal separation and parameter estimation. The hyperbolic tangent correlation of $r\left(t\right)$ and $r(t+\tau )$ is defined by

$$\begin{array}{cc}\hfill T_r(t,\tau )& =\mathrm{E}\left[tanh\left({\displaystyle \frac{|r\left(t\right)r^{*}(t+\tau )|}{H}}\right){\displaystyle \frac{r\left(t\right)r^{*}(t+\tau )}{|r\left(t\right)r^{*}(t+\tau )|}}\right]\hfill \\ \hfill & =\mathrm{E}\left[\sum _{n=1}^{\infty }{\displaystyle \frac{2^{2n}(2^{2n}-1)B_{2n}}{\left(2n\right)!}}{\left({\displaystyle \frac{|r\left(t\right)r^{*}(t+\tau )|}{H}}\right)}^{2n-1}{\displaystyle \frac{r\left(t\right)r^{*}(t+\tau )}{|r\left(t\right)r^{*}(t+\tau )|}}\right]\hfill \end{array}$$

(26)

where H and $B_n$ denote scale parameter and Bernoulli Numbers, respectively. Further, Bernoulli Numbers are computed by

$$\sum _{i=0}^n{B}_i{C}_{n+1}^i=0,B_0=1$$

(27)

According to the definition of the hyperbolic tangent correlation, the hyperbolic tangent correlation introduces odd harmonics. In addition, the hyperbolic tangent correlation has large computational complexity and is not suitable for use in algorithms for searching. Compared to the two methods above, PFLOM has the advantage of no harmonics and low computational complexity.

4. Simulation

4.1. Signal Generation and Parameter Settings

At the receiving end, there are two signals interfering with each other, i.e., $s_1\left(t\right)$ and $s_2\left(t\right)$. They are binary phase-shift keying (BPSK) signals with identical baud rates, mapping modes, and power, but different carrier frequencies, time delays, and Doppler shifts. Note that the proposed method can work for the same cyclic frequency related to carrier frequency. Here, we use a different carrier frequency to minimize the problem of occlusion in Figure 1. The additive noise $n\left(t\right)$ is modeled by S$\alpha$S distribution. Since the second-order statistics of $n\left(t\right)$ do not exist, we used the generalized power to replace the traditional one.

$$P_n^{\prime }={\gamma }^{\alpha },\beta =0,\delta =0.$$

(28)

Note that $\beta$ and $\delta$ must be set to zero. If $\beta$ and $\delta$ are not both zero, it is difficult to obtain randomly alpha-stable noise with the same generalized power. Consequently, the generalized signal-to-noise ratios, GSNR $(s_1/n)$ and GSNR $(s_2/n)$, are defined by

$$\begin{array}{cc}\hfill \mathrm{GSNR}(s_1/n)& =10lg\left({\displaystyle \frac{P_{s_1}}{P_n^{\prime }}}\right)\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \mathrm{GSNR}(s_2/n)& =10lg\left({\displaystyle \frac{P_{s_2}}{P_n^{\prime }}}\right)\hfill \end{array}$$

(30)

Similarly, GSNR $(a_1{s}_1/n)$ and GSNR $(a_2{s}_2/n)$ are defined by

$$\begin{array}{cc}\hfill \mathrm{GSNR}(a_1{s}_1/n)& =10lg\left({\displaystyle \frac{a_1^2{P}_{s_1}}{P_n^{\prime }}}\right)=20lg{a}_1+\mathrm{GSNR}(s_1/n)\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mathrm{GSNR}(a_2{s}_2/n)& =10lg\left({\displaystyle \frac{a_2^2{P}_{s_2}}{P_n^{\prime }}}\right)=20lg{a}_2+\mathrm{GSNR}(s_2/n)\hfill \end{array}$$

(32)

The parameters used in the simulation are listed in

Table 1. Parameter settings in the simulation.

Parameter	Symbol	Value
Baud rate	$B_1,B_2$	10 Baud, 15 Baud
Carrier frequency	$f_{\mathrm{c}}$	60 Hz, 30 Hz
Sampling frequency	$f_{\mathrm{s}}$	1200 Hz
Sampling Interval	$T_{\mathrm{s}}$	${1200}^{-1}$ s
Amplitude	$A_1,A_2$	10
Magnitude mismatch	$a_1,a_2$	$0.8$, $0.9$
Time delay	$d_1,d_2$	$13T_{\mathrm{s}}$, $20T_{\mathrm{s}}$
Doppler shift	$D_1,D_2$	$0.03f_{\mathrm{s}}$, $-0.04f_{\mathrm{s}}$
Fractional order	p	$0.7$
Generalized signal-to-noise ratio	GSNR	10 dB
Characteristic exponent	$\alpha$	$1.6$
Scale	$\gamma$	calculated by $\alpha$ and GSNR

.

4.2. One Single Experiment

To explain the proposed method deeply and clearly, one single experiment is carried out under the condition of $\alpha =1.6$ and $\mathrm{GSNR}=10\mathrm{dB}$. The preset parameter, i.e., the fractional order p, is set to $0.7$. Since $p<\alpha /2$, the proposed method is effective in the joint estimation on cyclic frequencies and Doppler shifts based on the characteristic of PFLOM dealing with alpha-stable noise, which is shown in in Figure 1. The four subfigures provide different views of $|Z_r^{\langle p\rangle }(\xi ,f)|$.

In Figure 1a, the highest peak $(0,-24,1)$ has no information about cyclic frequency, so it cannot be used for the joint estimation on cyclic frequency and Doppler shift. Consequently, the other four peaks are taken into consideration to achieve the estimation. $(-120,36,0.707322)$ and $(120,36,0.707322)$ are symmetric about the plane $\xi =0$. Similarly, $(-60,-48,0.634462)$ and $(60,-48,0.634462)$ are symmetric. The domain of cyclic frequency is $(-\infty ,\infty )$, and the domain of Doppler shift is also $(-\infty ,\infty )$. In practice, the range of Doppler shifts cannot be infinite, and we set it to be $[-600,600]$ instead. We could have ${\xi }_1=120$, $D_1=36$ for $s_1\left(t\right)$ and ${\xi }_2=60$, $D_2=-48$ for $s_2\left(t\right)$. Figure 1b shows the projection of the 3D graph onto the cyclic frequency domain, and it can be seen that the amplitudes of the five mentioned peaks are much higher than the others. The amplitude is directly affected by the corresponding magnitude mismatches $a_i$. The greater $a_i$ is, the higher the amplitude is. Figure 1c,d show the 2D top view of 3D graph. Figure 1c is with the labeling included as Figure 1a,b. Figure 1d is without any labeling to help show where the peaks are located.

4.3. Monte Carlo Experiments

According to the result and understanding from the single experiment above, one hundred Monte Carlo experiments are carried out to evaluate the average performance of the proposed PFLOM-based method compared to the traditional correlation-based method. The three main preset parameters—$\alpha$, p, and GSNR—are set as follows. $\alpha =1.0:0.1:2.0$, $p=0.2:0.2:0.8$, and $\mathrm{GSNR}=0:5:10$, which is different from the above single experiment. To demonstrate the significance of the peaks, the peak-to-average power ratio (PAPR) is used as the evaluation criterion shown in Figure 2. The definition of the PAPR for $s_i\left(t\right)$ is given by

$$\mathrm{PAPR}\left(s_i\right)=20lg\left({\displaystyle \frac{|Z_r^{\langle p\rangle }({\xi }_i,f_i)|}{\overline{|Z_r^{\langle p\rangle }(\xi ,f)|}}}\right)$$

(33)

where $|Z_r^{\langle p\rangle }({\xi }_i,f_i)|$ and $\overline{|Z_r^{\langle p\rangle }(\xi ,f)|}$ denote the amplitude of the peak and the mean value of the amplitudes, respectively.

In Figure 2, it can be seen that the traditional correlation-based method does well in Gaussian noise, and these experiment results are in perfect agreement with the classical theory, i.e., the second-order-based methods can effectively address or remove Gaussian noise. However, with the characteristic exponent $\alpha$ decreasing, the correlation-based method has severe performance degradation. Overall, the proposed PFLOM-based method with different fractional order p performs similarly well as the traditional one in Gaussian noise but has strong robustness against the decreasing of $\alpha$ (impulsive noise). The smaller the value of p, the more robust the proposed method is. These experiment results are also in perfect agreement that p should be smaller than or equal to half of $\alpha$$(p\le \alpha /2)$ to guarantee the effectiveness of the PFLOM.

In the six subfigures, the PFLOM-based method $(p=0.2)$ shows the best robustness. However, it does not mean a tiny value of p is optimal for all situations. In the localized zoomed-in images of Figure 2a,b, the PAPRs of the PFLOM-based method $(p=0.2)$ are $46.9370$ and $46.1015$, respectively, when $\alpha =2$, which are acceptably good results. Regrettably, they are both the worst results in their respective comparisons. With the increasing of GSNR, the PFLOM-based method $(p=0.4)$ demonstrates good overall performance compared to the PFLOM-based method $(p=0.2)$. In Figure 2c, the PFLOM-based method $(p=0.2)$ performs better in $\alpha \in [1.0,1.55]$, while the PFLOM-based method $(p=0.4)$ performs better in the other. In Figure 2c,d, the PFLOM-based method $(p=0.2)$ performs better in $\alpha \in [1.0,1.55]$, while the PFLOM-based method $(p=0.4)$ performs better in $\alpha \in [1.6,2.0]$. In Figure 2e,f, the PFLOM-based method $(p=0.2)$ performs better in $\alpha \in [1.0,1.35]$, while the PFLOM-based method $(p=0.4)$ performs better in $\alpha \in [1.4,2.0]$

Based on the simulation and comparison, there are four basic conclusions that can be drawn.

• Both GSNR and the characteristic exponent $\alpha$ are the crucial factors that affect the performance of the method.
• The traditional correlation-based method does well in Gaussian noise $(\alpha =2)$, but its performance is severely degraded under impulse noise $(\alpha <2)$.
• The proposed PFLOM-based method has strong robustness in impulsive noise when the fractional order inside satisfies $(p\le \alpha /2)$.
• One tiny value of p is acceptable, but it may not be optimal. GSNR should be taken into consideration to set the optimal value of p.

5. Conclusions

In this paper, we propose a novel framework that integrates the principles of cyclostationary signal processing with the PFLOM technique to achieve robust Doppler shift estimation in impulsive noise. By utilizing the estimated cyclic frequency, a fundamental characteristic of cyclostationary signals, our method effectively distinguishes between mixed cyclostationary signals. Through comprehensive Monte Carlo simulations, we have found that while the correlation-based method, which relies on the second-order statistical analysis, performs well under Gaussian noise, it is inadequate in the presence of impulsive noise. In contrast, our PFLOM-based method exhibits robust performance in both Gaussian and impulsive noises, with the additional advantage of a straightforward fractional order p setting, which enhances its practicality and potential for wider application across various fields.

We maintain a positive outlook regarding the proposed PFLOM-based method, which is expected to not only overcome the limitations inherent in current techniques for Doppler shift estimation but also to facilitate its application in various additional fields. The method’s robustness against impulsive noise, coupled with its straightforward implementation, positions it as a promising candidate for improving the accuracy and reliability of signal processing in critical applications, including radar systems, wireless communication, and remote sensing. As we persist in our exploration and refinement of this method, we anticipate its full potential will be realized, resulting in advancements and innovations within the realm of signal processing technologies.