Mathematical Analysis Power Spectrum of M-ary MSK and Detection with Optimum Maximum Likelihood

Abstract

In this paper, the power spectral density(PSD) for Multilevel Minimum Shift Keyed signal with modulation index h = 1/2 (M-ary MSK) are derived using the mathematical method of the Markov Chain model. At first, according to an essential requirement of the phase continuity characteristics of MSK signals, a complete model of the whole process of signal generation is built. Then, the derivations for autocorrelation functions are carried out precisely. After that, we verified the correctness and accuracy of the theoretical derivation by comparing the derived results with numerical simulations using MATLAB. We also divided the spectrum into four components according to the derivation. By analyzing these figures in the graphic, each component determines the characteristics of the spectrum. It is vital for enhanced spectral characteristics. To more visually represent the energy concentration of the main flap and the roll-down speed of the side flap, the specific out-of-band power of M-ary MSK is given. OMLCD(Optimum Maximum Likelihood Coherent Detection) of M-ary MSK is adopted to compare the signal received with prepared in advance in a code element T to go for the best. And M-ary MSK BER(Bit Error Rate) is compared with the same ary PSK (Phase Shift Keying) with M=2,4,6,8. The results show the detection method could improve performance by increasing the length of L(memory inherent) in the phase continuity.

1. Introduction

In this paper, we analyzed the M-ary MSK communication system in two aspects: spectral properties of M-ary MSK (Multilevel Minimum Shift Keyed) signal with modulation index h=1/2 and optimum performance of optimum maximum likelihood coherent detection in BER. We derive the PSD expression using the mathematical method of the Markov Chain model to make the process more streamlined. The results can help us search the key points specific to a particular formula that affects the spectral characteristics and facilitates spectrum modification. The simulation of M-ary MSK in BER can be used to assess the stability of communication system transmissions. The MSK signals are very widely used in the VLF(Very Low Frequency) communication field. The benefits of the modulated signal possess a consistent envelope and a high bandwidth efficiency. Chen Xingbo [1] changed MSK’s carrier to LFM (Linear Frequency Modulation) signal, MSK-LFM. He analyzed the energy leakage of MSK-LFM based on the spectrum in expression and simulation. Dong [2] demonstrated a receiver using PLL(Phase-Locked Loop) to get a high-speed transmission with a low BER. Wang [3] derived the BER of MSK modulation based on the gamma-gamma channel model with coherent detection and analyzed the BER performance of MSK in weak, medium, and strong turbulence. This continuous phase limitation affects the signal with two crucial aspects. At first, with the transient effects, lessened spectral bandwidth will get better. Secondly, the continuous phase will improve demodulation performance by making decisions with several symbols rather than one symbol only [4]. M-ary MSK also has the advantages of MSK apply for Very Low Frequency to ensure no excessive power during phase shift. It can not only increase the communication rate but also improve spectral efficiency. Le Nguyen Binh applied the M-ary MSK modulation scheme for long-haul optically amplified transmission systems. According to the numerical result, a significant increase is achieved in transmission rates [5]. For a variety of MSK-derived communication systems, M-ary MSK is used to improve spectral efficiency. For example, Tao analyzed IRZ-MSK in a high-speed optical fiber transmission system from [6], research can be continued on multi-decimal.

However, using multi-decimal symbols instead of binary symbols within a code element can result in a significant increase in transmission rate [7]. While improving transmission efficiency, analyzing the spectral characteristics has become an essential part. Anderson and Salz [8] considered getting the PSD of M-ary CPFSK(Continuous Phase Frequency Shift Keying) by directly calculating the Fourier transform over a long time interval signal waveform. In addition, we derive the theoretical expressions of M-ary MSK with h=1/2 and present the compositions of M-ary MSK, which is vital for analyzing the spectrum structure and use in communication system design for improving the performance [9]. Ekanayake [10] derived the spectral properties of M-ary CPFSK with h=1/M by extending reference [9], which the method developed for binary digital FM(Frequency-Modulated) gives the PSD using a wave train made up of sine-wave segments. Obviously, as the binary number increases, the modulation index also rises, reducing the phase difference between adjacent characters, thus leading to a higher chance of demodulation errors. Biglieri used the representation of Rimoldi [11] to calculate a closed-form expression for the power spectral density of full-response M-ary CPM with modulation index J/M by straightforward computations, J an integer. The direct derivation process can be tedious and challenging to understand than building a model.

M-ary MSK compared to MSK, the transmission efficiency improves significantly. Considering the characteristics of MSK signals, the spectrum of M-ary MSK is analyzed by the Markov Chain model at the same time. Reasons for MSK modulation process can be built with Markov Chain process for spectral analytical analysis are as follows:

1. The signal sent by the MSK source is completely statistically independent, which corresponds to the Markov Chain with no after-effects.

2. With the fixed modulation index 1/2 and periodicity of trigonometric functions in phase, the number of states is inevitable, corresponding to a Markov Chain state space.

3. The transfer between states can accordingly be expressed in terms of the transfer probabilities of Markov Chains [12, 13, 14].

Simon and Divsalar [15] have discussed the performance of coherent detection for MPSK and give the structure of making decisions on two successive symbols. However, Park, Lee, and Feher [16], following the CPM(Continuous Phase Modulation) characteristics, proposed an optimum ML(Maximum-Likelihood) coherent receiver for unique CPM signals, while Andersson, Svensson [17] considered a noncoherent detector with an environment of white Gaussian noise, and they proposed a decoding algorithm for these convolutional coding CPM schemes. In 1992, Jung, Baier proposed a new novel signal generation concept for CPM with h=1/2 based on real-valued impulses. They derived the expressions in mathematical for both the real impulses and the bandpass CPM signals[18].

2. Signal Expression

MSK is a special CPFSK with h = 1/2. We can express the i-th phase of M-ary MSK signal as [19]:

\(\Phi_{i}=\frac{d_{i} \pi h(t-(i-1) T)}{T}+\pi h \sum_{j=1}^{i-1} d_{j}+\theta \quad(i-1) T       (1)

So the signal can be simplified and written as:

\(r(t)=A \cos \left(\Phi_{i}\right) \quad(i-1) T       (2)

where

A The amplitude of the signal

h Modulation index

d_i Digital information signal

T The whole symbol interval

Θ Starting phase

It is clear that when \(i=1, \sum_{j=1}^{i-1} d_{j}=0\). We set the modulation index h =1/2 , and the amplitude A is normalized. In each extended interval, one value of numeric information is assumed from the set \(\{\pm 1, \pm 3, \pm 5, \ldots, \pm(\mathrm{M}-1)\}\). According to (1), we can see clearly that the difference with PSK modulation is that the current phase is related to the current input information and the phase of the last moment to ensure the phase continuity. All possible phase trajectories for all values of the phase term Φ_i are shown in Fig. 1. For 4-ary MSK modulation, it is assumed that the starting phase θ is known, so it is generally considered to be 0, and there is no loss in generality. Therefore we set θ = 0 and d_i = ±1,±3 . When the next phase with digital information input increases over 2π or decreases under −2π , it can always be equated to one of the existing state trajectories because of the periodicity of the cosine function.

Fig. 1. Phase trellis of 4-ary MSK

3. The Derivation of M-ary MSK Autocorrelation Function in Mathematical Expressions

As mentioned in the Introduction, we apply the Markov Chain Process Strategy to derive the analytical expression of the spectrum of M-ary MSK. At first, we evaluate the autocorrelation functions of M-ary MSK. Then we get the theoretical spectrum expressions. Markov Chain Strategy is based on references [13] [20].

3.1 Markov Chain Process Strategy

We assumed X(t) = X_R(t) + jX_I(t) is random M-ary MSK signal because within one symbol interval T_b digital information takes one of the set values \(\left\{\pm D_{1}, \pm D_{2}, \pm D_{3}, \ldots, \pm D_{d}\right\}=\{\pm 1, \pm 3, \pm 5, \ldots, \pm(\mathrm{M}-1)\}\) with equal probability. So that X_R(t) and X_I(t) have 4M types of states in every symbol interval, Those states are:

\(\begin{gathered} \sin \Phi\left(\pm D_{1}\right), \sin \Phi\left(\pm D_{2}\right), \sin \Phi\left(\pm D_{3}\right), \ldots, \sin \Phi\left(\pm D_{d}\right) \\ \cos \Phi\left(\pm D_{1}\right), \cos \Phi\left(\pm D_{2}\right), \cos \Phi\left(\pm D_{3}\right), \ldots, \cos \Phi\left(\pm D_{d}\right) \end{gathered}\)

Because of the symmetry of the I and Q channels, we only take the actual envelope of the I channel signal for analysis. The types of X_R(t) states of each symbol interval are decimated to 2M :

\(\sin \Phi\left(D_{d}\right), \sin \Phi\left(D_{d}+1 / 2 \pi\right), \sin \Phi\left(D_{d}+\pi\right), \sin \Phi\left(D_{d}+3 / 2 \pi\right)\)

That is : 

\(\sin \Phi\left(D_{d}\right), \cos \Phi\left(D_{d}\right),-\cos \Phi\left(D_{d}\right),-\sin \Phi\left(D_{d}\right)\)       (3)

The states are shown in equation (3) is arranged from 1 to d by D_d and gather each different state groups numbering the states in order states 1 to state 2M .

Because the signals sent by the sources are statistically independent and equal in probability. For each of the 2M states, the starting state probability is

\(P(i)=\frac{1}{N}=\frac{1}{2 M}, i \in[1,2 M]\)       (4)

N is the number of starting states. Under a different number of transfers m , p(j|i,m) represents the probability of transition from state i to state j . Put them into the matrix to get its state transfer matrix \({P_{M}}^{m}\):

a) In the case m = 0 :

\(P_{M}{ }^{0}=\overbrace{\left(\begin{array}{lll} 1 & & \\ & \ddots & \\ & & 1 \end{array}\right)}^{2 M}\)       (5)

p(j|i,m) corresponds to the element \(p^{(m)}_{i,j}\) in the matrix \(P_{M}^{m}\). The superscript m represents the number of transfer steps, and the subscript M denotes the M-ary.

b) In the case m =1:

When a one-step transfer is performed, the transfer probability of each state can be obtained from the phase initiation, termination characteristics of each state signal, and the phase continuity of the signal. We can get the state transfer matrix at m =1. The following Fig. 2 shows the 8-ary MSK state figure, from which the general rule can be obtained.

\(\begin{gathered} P_{M}{ }^{1}=\frac{M}{2}\{\overbrace{\left(\left(\begin{array}{ccccc} A & A & \cdots & A & A \\ B & B & \cdots & B & B \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ A & A & \cdots & A & A \\ B & B & \cdots & B & B \end{array}\right)\right.}^{M / 2} \\ A=\left(\begin{array}{cccc} 0 & a & 0 & 0 \\ b & 0 & 0 & b \\ b & 0 & 0 & b \\ 0 & 0 & a & 0 \end{array}\right) B=\left(\begin{array}{cccc} 0 & 0 & a & 0 \\ b & 0 & 0 & b \\ b & 0 & 0 & b \\ 0 & a & 0 & 0 \end{array}\right) \end{gathered}\)       (6)

Fig. 2. Different states of 8-ary MSK

The values of a and b depend on the number of a or b in each row, and the sum of probability each row should be equal to 1.

c) In the case m = 2 :

\(P_{M}{ }^{2}=\left(P_{M}{ }^{1}\right)^{2}=\frac{M}{2}\{\overbrace{\left(\begin{array}{cccc} \frac{M}{4}\left(A^{2}+\mathrm{AB}\right) & \cdots & \frac{M}{4}\left(A^{2}+\mathrm{AB}\right) \\ \frac{M}{4}\left(\mathrm{~B}^{2}+\mathrm{AB}\right) & \cdots & \frac{M}{4}\left(\mathrm{~B}^{2}+\mathrm{AB}\right) \\ \vdots & \ddots & \vdots \\ \frac{M}{4}\left(A^{2}+\mathrm{AB}\right) & \cdots & \frac{M}{4}\left(A^{2}+\mathrm{AB}\right) \\ \frac{M}{4}\left(\mathrm{~B}^{2}+\mathrm{AB}\right) & \cdots & \frac{M}{4}\left(\mathrm{~B}^{2}+\mathrm{AB}\right) \end{array}\right)}\\ A^{2}=B^{2}=A B=\left(\begin{array}{cccc} a b & 0 & 0 & a b \\ 0 & a b & a b & 0 \\ 0 & a b & a b & 0 \\ a b & 0 & 0 & a b \end{array}\right)\)       (7)

In Fig. 2, take 8-ary MSK as an example. We can see that the first to fourth columns of the figure correspond to each of the four types of states in (3). The fact that state 1 ends magnitude one means that its next state must start withmagnitude one, the various states in the second column. It is easy to find the matrix of states arranged in this pattern, summarize to (6) and (7). To better understand matrix composition, examples of MSK, 4-ary MSK and 8-ary MSK are listed in Appendix 1.

Based on knowledge of stochastic processes, the autocorrelation function is given as follow expression from [21]:

\(R(\tau)=E\left[X_{I}(t+\tau) X_{I}(t)\right]=R\left(m T+\tau^{\prime}\right)\)

The above equation is the statistical average of the two random variables corresponding to the values taken by the stochastic process XI at two moments separated by τ . Where τ = mT + τ′ since X_I is cyclostationary stochastic processes. In terms of the transition relationship between states and the corresponding transition probability, R(τ) can be expressed as:

\(\begin{aligned} &R(\tau)=R\left(m T+\tau^{\prime}\right)=\sum_{i=1}^{N} \sum_{j=1}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)+\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) \mathrm{s}_{j i}{ }^{*}\left(\tau^{\prime}\right)\right] \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N}\left\{P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right]+P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) \mathrm{s}_{j i}{ }^{*}\left(\tau^{\prime}\right)\right]\right\} \end{aligned}\)       (8)

where

\(\begin{aligned} &S_{i j}\left(T-\tau^{\prime}\right)=\frac{1}{T} \int_{0}^{T-\tau^{\prime}} S_{i}(u) s_{j}^{*}\left(u+\tau^{\prime}\right) d u \\ &S_{j i}^{*}\left(\tau^{\prime}\right)=\frac{1}{T} \int_{0}^{\tau^{\prime}} S_{j}^{*}(u) s_{i}\left(u+T-\tau^{\prime}\right) d u \end{aligned}\)

3.2 Case 1: Where 𝟎 ≤ 𝐭 ≤ 𝐓 and 𝐦 = 𝟎

Based on the above Markov Chain Process Strategy, evaluate the autocorrelation functions of M-ary MSK. First, we calculate the autocorrelation function \({{R^{M}}_{m=o}} (\tau)\) in a fixed time interval [0,T]. When m = 0 , equation (8) can be expressed as follows.

\(\begin{aligned} {{R^{M}}_{m=0}}(\tau) &={{R^{M}}_{m=0}}\left(m T+\tau^{\prime}\right)=\sum_{i=1}^{N} \sum_{j=1}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)+\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) \mathrm{s}_{j i}{ }^{*}\left(\tau^{\prime}\right)\right] \\ &=\sum_{i}^{N} \sum_{j}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right]+\sum_{i}^{N} \sum_{j}^{N} P(i)\left[p(j \mid i, m+1) s_{j i}{ }^{*}\left(\tau^{\prime}\right)\right] \\ &=R 1^{M}{ }_{m=0}(\tau)+R 2^{M}{ }_{m=0}(\tau) \end{aligned}\)

To facilitate the operation, the autocorrelation function \({R^{M}_{m=0}}(\tau)\) can be divided into the sum of two parameters \({R1^{M}_{m=0}}(\tau)\) and \({R2^{M}_{m=0}}(\tau)\) for the procedure. Before that, since the excessive number of terms in the expansion formula, patterns can only be found in the calculation process in 4-ary and 8-ary first.

\(\begin{aligned} {{R^{M}}_{m=0}}(\tau) &=\frac{1}{2 M}\left\{2 s_{11}+2 s_{22}+2 s_{55}+2 s_{66}+\frac{1}{2} s_{21}{ }^{*}-\frac{1}{2} s_{25}{ }^{*}+\frac{1}{2} s_{61}{ }^{*}-\frac{1}{2} s_{65}{ }^{*}\right\} \\ &+\frac{1}{2 M}\left\{\frac{1}{2} s_{10,1}{ }^{*}-\frac{1}{2} s_{10,5}{ }^{*}+\frac{1}{2} s_{14,1}{ }^{*}-\frac{1}{2} s_{14,5}{ }^{*}\right\} \\ &+\frac{1}{2 M}\left\{2 s_{99}+2 s_{10,10}+2 s_{13,13}+2 s_{14,14}+\frac{1}{2} s_{29}{ }^{*}-\frac{1}{2} s_{2,13}{ }^{*}+\frac{1}{2} s_{6,9}{ }^{*}-\frac{1}{2} s_{6,13}{ }^{*}\right\} \\ &+\frac{1}{2 M}\left\{\frac{1}{2} s_{10,9}{ }^{*}-\frac{1}{2} s_{10,13}{ }^{*}+\frac{1}{2} s_{14,9}{ }^{*}-\frac{1}{2} s_{14,13}{ }^{*}\right\} \\ &+\frac{1}{2 M}\{\cdots\} \end{aligned}\)

The partial expansion of \({{R^{M}}}_{m=0}(\tau)\) is shown above, with the first row being the result of a 4-ary calculation and the first four rows being the result of an 8-ary calculation, thus concluding that:

\(\begin{aligned} R 1^{M}{ }_{m=0}(\tau) &=\sum_{i}^{N} \sum_{j}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right]=\frac{1}{M}\left[s_{11}+s_{22}+s_{33}+\cdots+s_{M M}\right] \\ &=\frac{2}{M}\left[\frac{1}{\mathrm{~T}} \int_{0}^{\mathrm{T}-\tau^{\prime}} \sin \frac{n \pi}{2 \mathrm{~T}} u \sin \frac{n \pi}{2 \mathrm{~T}}\left(\mathrm{u}+\tau^{\prime}\right) \mathrm{du}\right] \\ &=\frac{2}{M}\left[\frac{1}{2}\left(1-\frac{\tau^{\prime}}{\mathrm{T}}\right) \cos \frac{\pi}{2 \mathrm{~T}} \tau^{\prime}+\frac{1}{2}\left(1-\frac{\tau^{\prime}}{\mathrm{T}}\right) \cos \frac{3 \pi}{2 \mathrm{~T}} \tau^{\prime}+\cdots+\frac{1}{2}\left(1-\frac{\tau^{\prime}}{\mathrm{T}}\right) \cos \frac{n \pi}{2 \mathrm{~T}} \tau^{\prime}\right] \\ &=\frac{2}{M} \sum_{n=1}^{M-1} \frac{1}{2}\left(1-\frac{\tau^{\prime}}{T}\right) \cos \frac{n \pi \tau^{\prime}}{2 T}, \quad n=1,3,5 \ldots \end{aligned}\)

The result of integration is different for different n , m terms with various positive and negative signs. \({{R2^{M}}}_{m=0}(\tau)\) is the result of the accumulation of the terms. What’s more, we found that differences in the initial phases of cos and sin within the integral lead to different generalizations of the integral result. Therefore, the \({{R2^{M}}}_{m=0}(\tau)\) is divided into two terms according to whether the n , m is the same or not. \(\Sigma^{M=1}_{n=1}\) denotes an odd summation of n from 1 to (M -1).

\(\begin{aligned} &R 2^{M}{ }_{m=0}(\tau)=\sum_{i}^{N} \sum_{j}^{N} P(i)\left[p(j \mid i, m+1) s_{j i}{ }^{*}\left(\tau^{\prime}\right)\right]=R 2_{-} \operatorname{{{same}^{M}}_{m=0}}(\tau)+R 2_{-} \operatorname{{dif}^{M}}_{m=0}(\tau) \\ &=\frac{2}{M^{2}} *\left\{\begin{array}{l} (-1)^{(n+1) / 2} \sum_{n}^{M-1}\left\{\frac{1}{T} \int_{0}^{\tau^{\prime}} \cos \frac{n \pi}{2 T} u \sin \frac{n \pi}{2 T}\left(u+T-\tau^{\prime}\right) d u\right\}, n=1,3,5 \ldots \\ (-1)^{(m+1) / 2} \sum_{n=1}^{M-1} \sum_{m=1}^{M-1} \frac{1}{T} \int_{0}^{\tau^{\prime}} \cos \frac{n \pi}{2 T} u \sin \frac{m \pi}{2 T}\left(u+T-\tau^{\prime}\right) d u, n \neq m \text { and } n, m=1,3,5 \ldots \end{array}\right. \\ &=\frac{2}{M^{2}} *\left\{\begin{array}{l} (-1)^{(n+3) / 2} \sum_{n=1}^{M-1}\left\{\frac{1}{2} \frac{\tau^{\prime}}{T} \sin \frac{\pi\left(T-\tau^{\prime}\right) n}{2 T}+\frac{\sin \left(\frac{1}{2} \pi n\right) \sin \frac{\pi \tau^{\prime} n}{2 T}}{\pi n}\right\}, n=1,3,5 \ldots \\ (-1)^{(m+1) / 2} \sum_{n=1}^{M-1} \sum_{m=1}^{M-1} 2 * \frac{n\left[\sin \left(\frac{1}{2} \pi m\right) \sin \frac{\pi \tau^{\prime} n}{2 T}\right]-m\left[\cos \frac{\pi\left(T-\tau^{\prime}\right) m}{2 T}\right]}{\pi\left[m^{2}-n^{2}\right]}, n \neq m \text { and } n, m=1,3,5 \ldots \end{array}\right. \end{aligned}\)       (9)

3.3 Case2: Where 𝐓 ≤ 𝐭 ≤ 𝟐T and 𝐦 = 𝟏

The calculation process is similar to that of Case 1. \({{R^{M}}}_{m=1}(\tau)\) is divided into \(\text{R1_}{{\text{same}^{M}}}_{m=1}(\tau)\)  and \(\text{R1_}{{\text{dif}^{M}}}_{m=1}(\tau)\) for the summation calculation.

\(\begin{aligned} {{R^{M}}_{m=1}}(\tau) &={{R^{M}}_{m=1}}(\tau)\left(m T+\tau^{\prime}\right)=\sum_{i=1}^{N} \sum_{j=1}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)+\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) {\mathrm{s}_{j i}}^{*}\left(\tau^{\prime}\right)\right] \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N}\left\{P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right]+P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) {\mathrm{s}_{j i}}^{*}\left(\tau^{\prime}\right)\right]\right\} \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N} P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right] \end{aligned}\)

Since the high symmetry and repetition of the matrix, the \(P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1){ \mathrm{s}_{j i}}^{*}\left(\tau^{\prime}\right)\right]\) term is eliminated during the calculation. Summing up the items gives:

\(\begin{aligned} &{R^{M}}_{m=1}(\tau)=R 1^{M}{ }_{m=1}(\tau) \\ &=R 1_{-} {\operatorname{same}^{M}}_{m=1}(\tau)+R 1_{-} \operatorname{dif}^{M}{ }_{m=1}(\tau) \\ &=\frac{2}{M^{2}} *\left\{\begin{array}{l} (-1)^{(n+1) / 2} \sum_{n}^{M}\left\{\frac{1}{2}\left(1-\frac{\tau^{\prime}}{T}\right) \sin \frac{\pi \tau^{\prime} n}{2 T}-\frac{\sin \left(\frac{1}{2} \pi n\right) \sin \frac{\pi\left(T-\tau^{\prime}\right) n}{2 T}}{n \pi}\right\}, \mathrm{n}=1,3,5 \ldots \\ (-1)^{(n+3) / 2} \sum_{n}^{M} \sum_{m}^{M} \frac{(n-m) \cos \frac{\pi\left[T m+\left(T-\tau^{\prime}\right) n\right]}{2 T}-2 n \cos \frac{\pi \tau^{\prime} m}{2 T}+(m+n) \cos \frac{\pi\left(T m+\left(\tau^{\prime}-T\right) n\right)}{2 T}, \ldots}{\left(m^{2}-n^{2}\right) \pi},n \neq m \text { and } n, m=1,3,5 \end{array}\right. \end{aligned}\)

Bring \(\tau^{\prime}=\left\{\begin{array}{lll} \tau & m=0 & (0 \leq \tau \leq \mathrm{T}) \\ \tau-\mathrm{T} & m=1 & (\mathrm{~T} \leq \tau \leq 2 \mathrm{~T}) \end{array}\right.\) into (10) to get the autocorrelation function, so the final M-ary MSK expression is:

\(R^{M}(\tau)={R^{M}}_{m=0}(\tau)+{R^{M}}_{m=1}(\tau)={R 1^{M}}_{m=0}(\tau)+{R 2^{M}}_{m=0}(\tau)+{R^{M}}_{m=1}(\tau)\)       (11)

Where T is a time interval. Because of the symmetry of the correlation function, the expressions in −2T ≤ t ≤ -T and −T ≤ t ≤ 0 are equal to equations T ≤ t ≤ 2 and 0 ≤ t ≤ T , respectively.

4. PSD function

Once the autocorrelation function R^M(τ) in section 3 has been evaluated. The PSD function can be derived directly by performing a Fourier transform on the autocorrelation function. Here is the derivation process:

\(\begin{aligned} S^{M} &=F\left\{R^{M}(\tau)\right\}=\int_{-\infty}^{\infty} R^{M}(\tau) e^{-j \omega \tau} d \tau=2 \int_{-\infty}^{\infty} R^{M}(\tau) \cos 2 \pi f \tau d \tau \\ &=2 \int_{-\infty}^{\infty} {R^{M}}_{m=0}(\tau) \cos 2 \pi f \tau d \tau+2 \int_{-\infty}^{\infty} {R^{M}}_{m=1}(\tau) \cos 2 \pi f \tau d \tau \\ &={S^{M}}_{m=0}+{S^{M}}_{m=1} \end{aligned}\)       (12)

Based on (12), the PSD function can easily be derived by direct calculation, and the detailed derivations are moved to Appendix 2. The following cases are solved for \({S^{M}}_{m=0}\) and \({S^{M}}_{m=1}\) .

4.1 Case 1: Where 𝟎 ≤ 𝐭 ≤ 𝐓 and 𝐦 = 𝟎, \({S^{M}_{m=0}}\)

\(\begin{gathered} {S^{M}}_{m=0}(f)=\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau+\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) {\mathrm{s}_{j i}}^{*}\left(\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau \\ =S 1^{M}{ }_{m=0}(f)+S 2_{-} {\operatorname{same}^{M}}_{m=0}(f)+S 2_{-} d i f^{M}{ }_{m=0}(f) \\ S 1^{M}{ }_{m=0}(f)=T\left\{\frac{\sin [\varphi(n)+\omega(f)]}{2[\varphi(n)+\omega(f)]}+\frac{\sin [\varphi(n)-\omega(f)]}{2[\varphi(n)-\omega(f)]}-\frac{\cos [\varphi(n)-\omega(f)]+[\varphi(n)-\omega(f)] \sin [\varphi(n)-\omega(f)]-1}{2[\varphi(n)-\omega(f)]^{2}}\right. \\ \left.-\frac{\cos [\varphi(n)+\omega(f)]+[\varphi(n)+\omega(f)] \sin [\varphi(n)+\omega(f)]-1}{2[\varphi(n)+\omega(f)]^{2}}\right\} \\ S 2_{-} {\operatorname{same}^{M}}_{m=0}(f)=4 T \frac{\varphi(n)\left[\cos \omega(f) * S_{2}+8 \omega(f) \sin \omega(f)\right]-\sin \varphi(n) * S_{1}(n)}{S_{2}(n)} \\ -4 T \frac{\sin \varphi(n)[\omega(f) \sin \omega(f) \sin \varphi(n)+\varphi(n) \cos \omega(f) \cos \varphi(n)-\varphi(n)]}{\varphi(n) S_{2}(n)} \\ S 2_{-} {d i f^{M}}_{m=0}(f)=8 T \frac{\varphi(n) \sin \varphi(m) *[\omega(f) \sin \omega(f) \sin \varphi(n)+\varphi(n) \cos \omega(f) \cos \varphi(n)-\varphi(n)]}{[\varphi(n)-\varphi(m)][\varphi(n)+\varphi(m)] S_{2}(n)} \\ +8 T \varphi(m) \frac{\varphi(m) \sin \varphi(m)-\omega(f) \sin \omega(f)}{[\varphi(n)-\varphi(m)][\varphi(n)+\varphi(m)] S_{2}(m)} \end{gathered}\)

 4.2 Case 2: Where T ≤ 𝐭 ≤ 2𝐓 and 𝐦 = 1, \({S^{M}_{m=1}}\)

\(\begin{aligned} {S^{M}}_{m=1} &=\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau \\ &=S 1_{-} {\text {same }^{M}}_{m=1}(f)+S 1_{-} {d i f^{M}}_{m=1}(f) \end{aligned}\)

\({S^{M}}_{m=1}\) si divided into \(S 1_{-} {\text {same }^{M}}_{m=1}(f)\) and \(S 1_{-} {\text {dif }^{M}}_{m=1}(f)\) for the summation calculation.

\(\begin{aligned} S_{-} {\operatorname{same}^{M}}_{m=1}(f) &=-2 T \frac{2[\varphi(n)+\omega(f)] \cos [\varphi(n)-2 \omega(f)]+2[\varphi(n)-\omega(f)] \cos [\varphi(n)+2 \omega(f)]-4 \varphi(n) \cos \omega(f)}{S_{2}(n)} \\ &+2 T\left\{\frac{2[\varphi(n)-\omega(f)] \cos [\varphi(n)-2 \omega(f)]-\sin [\varphi(n)-2 \omega(f)]}{4[\varphi(n)-\omega(f)]^{2}}\right.\\ &\left.+\frac{[\varphi(n)+\omega(f)] \cos [\varphi(n)+2 \omega(f)]-\sin [\varphi(n)+2 \omega(f)]}{4[\varphi(n)+\omega(f)]^{2}}-\frac{2 \varphi(n)\left[S_{2} \cos \omega(f)+8 \omega(f) \sin \omega(f)\right]}{S_{2}^{2}(n)}\right\} \\ &+2 T \sin \varphi(n)\left[\frac{\cos [\varphi(n)-\omega(f)]+\cos [\varphi(n)+\omega(f)]-2 \cos 2 \omega(f)}{S_{2}(n)}-\frac{2 \omega(f) \sin \omega(f) \sin \varphi(n)}{\varphi(n) S_{2}(n)}\right] \\ S 1_{-} {d i f^{M}}_{m=1}(f)=& \frac{T}{2}\left\{\frac{\sin [\varphi(m)+2 \omega(f)]-\sin [\varphi(n)+\varphi(m)+\omega(f)]}{[\varphi(n)+\varphi(m)][\varphi(n)-\omega(f)]}-\frac{\sin [2 \omega(f)-\varphi(m)]+\sin [\varphi(n)+\varphi(m)-\omega(f)]}{[\varphi(n)+\varphi(m)][\varphi(n)+\omega(f)]}\right\} \\ &-T \varphi(n) \frac{[\varphi(m)-\omega(f)] \sin [\varphi(m)+2 \omega(f)]+[\varphi(m)+\omega(f)] \sin [\varphi(m)-2 \omega(f)]+2 \omega(f) \sin \omega(f)}{\left[\varphi^{2}(m)-\varphi^{2}(n)\right]\left[\varphi^{2}(m)-\omega^{2}(f)\right]} \\ &+T\left\{\frac{\sin [\varphi(m)+2 \omega(f)]-\sin [\varphi(m)-\varphi(n)+\omega(f)]}{2[\varphi(m)-\varphi(n)][\varphi(n)+\omega(f)]}-\frac{\sin [\varphi(m)-\varphi(n)-\omega(f)]-\sin [\varphi(m)-2 \omega(f)]}{2[\varphi(m)-\varphi(n)][\varphi(n)-\omega(f)]}\right\} \end{aligned}\)

where

\(\begin{aligned} &S_{1}(n)=4\left[\varphi^{2}(n)+\omega^{2}(f)\right]=\pi^{2}\left(n^{2}+16 f^{2} T^{2}\right) \\ &S_{2}(n)=4\left[\varphi^{2}(n)-\omega^{2}(f)\right]=\pi^{2}\left(n^{2}-16 f^{2} T^{2}\right) \\ &S_{1}(m)=4\left[\varphi^{2}(m)+\omega^{2}(f)\right]=\pi^{2}\left(m^{2}+16 f^{2} T^{2}\right) \\ &S_{2}(m)=4\left[\varphi^{2}(m)-\omega^{2}(f)\right]=\pi^{2}\left(m^{2}-16 f^{2} T^{2}\right) \\ &\omega(f)=2 \pi f T \qquad \varphi(m)=m \pi / 2 \qquad \varphi(n)=n \pi / 2 \qquad n, m=1,3,5 \ldots \end{aligned}\)

5. Simulation for Verification and Analysis

The expression results are compared with those signal simulations using MATLAB to verify the accuracy of autocorrelation functions and PSDs. The result is shown in Fig. 3, Fig. 4, Fig. 5, and Fig. 6. We set simulation parameters as the length of the simulated input signal n =10⁵ , code element time interval T = 0.1s , bit time interval  T_b = T / log2M , sampling frequency fs = 10⁵Hz .

Fig. 3. Comparison of M-ary MSK autocorrelation function, M=2,4,8,16

Fig. 4. Comparison of theoretical expression curve and signal simulation curve of 8-ary MSK autocorrelation function

Fig. 5. Theoretical expression curve and signal simulation curve of 4-ary MSK power spectraldensity function

Fig. 6. Comparison of theoretical expression curve and signal simulation curve of 16-ary MSK power spectral-density function

5.1 autocorrelation function of M-ary MSK

In Fig. 3, it can be concluded that the larger the value of M , the more the autocorrelation function is concentrated from −T to T . As a result, the amplitude between T to 2T and −2T to −T is getting lower, closer to 0. We can still conclude that there is no discrepancy between the theoretical autocorrelation function expressions in (11) and obtained from the simulation. And to understand the subtle differences between the theoretical expressions and simulation, Fig. 4 is drawn to compare the details of the octal signal as an example.

In Fig. 4, we can see that the overall shape has not changed in most of the exact matches, but there are subtle differences in the magnitude. Although this property is verified with little doubt, one can reasonably conclude that both two have approximately the same. In the simulation process, it was found that as the length of the simulated input signal n increases, the difference between the two curves will become lower, and the simulated curve will gradually approach the theoretical curve. So we believe that the reason for the slight difference in the curves is the experimental error.

5.2 Power spectral-density function

As 8-ary MSK autocorrelation function, we take the spectrum of 4-ary and 8-ary MSK as an example to verify the accuracy. Fig. 5 and Fig. 6 show that the simulation results match highly with the theoretical derivation in the amplitude of the main lobe and the first side lobe, after which distortion occurs because of the cut-off effect of the Fourier transform of the signal simulation.

The number of bits increases within a code element interval. This leads to the dispersion of energy in the spectrum. So it is necessary to normalize the energy to one bit when making the comparison. The horizontal coordinate is set to fT_b .

The PSD for M = 2,4,8,16 is shown in Fig. 7. We can see that at M = 4 and 8 the amplitude spectrum becomes more concentrated, and the signal intensity increase more. As well known, the PSD of higher-level M rolls off much quicker than that M = 2.

Fig. 7. Normalised PSDs of M-ary MSK,M=2,4,8,16

When M = 16, the top of the lobe is flatter so that the bandwidth can be better utilized and the transmission efficiency becomes improved, but h is fixed equal to 0.5. The spectral main lobe bandwidth of M-ary MSK becomes wider substantially with the number of levels up.

In Fig. 8, the logarithmic form of the spectrum function can be seen. It is worth mentioning that the first side lobe of 4-ary MSK is lower than that of MSK, and the main lobe is higher, and the energy is more concentrated.

Fig. 8. Logarithmic form of normalized power spectral densities of M-ary MSK, M=2,4,8,16

5.3 Spectrum composition

Fig. 9 shows the composition of the spectrum of M-ary MSK where M = 2,4,8,16. Because for the case of m = 1, the spectrum with many higher-order components whose amplitude becomes extremely low and the proportion of the overall, so it is considered a variable for comprehensive analysis. For the case of m = 0, the proportion is more significant, and the subdivision is more detailed.

Fig. 9. Signal composition of MSK,M=2,4,8,16

As far as those figures, we can provide an analysis of the composition of the spectrum in detail. The first one we can see is \(S 1^{M}{ }_{m=0}(f)\) contributed the majority of the spectrum and is more pronounced when the level M increases. Second, \(S 2_{-} {\operatorname{dif}^{M}}_{m=0}(f)\) determine the shape of the top of the spectrum waveform, and the fluctuations become smaller with the increase of M, which is why the top of the lobe gets flattered within a specific bandwidth.

5.4 Fractional out-of-band power

It is also helpful to show the out-of-band power with M-ary MSK. It defines as [10]:

\(P_{o b}=1-\left\{\int_{-B}^{+B} S(f) d f / \int_{-\infty}^{+\infty} S(f) d f\right\}\)       (13)

In Fig. 10, P_ob is plotted against BT_b for M-ary MSK where M = 2,4,8,16. We can see 4-ary MSK has less out-of-band power compared to 2-ary MSK when BT_b > 0.71 , while  BTb > 0.84, 8-ary MSK has less than 2-ary MSK but upper than 4-ary MSK.

Fig. 10. Fractional out-of-band power of MSK, M=2,4,8,16

6. Optimum Maximum Likelihood Coherent Detection

The M-ary MSK detection in this paper uses coherent demodulation with maximum likelihood probability. In a code-time T, the phase state will increase or decrease 1/2πD_d with the input data d_i based on the starting phase θ . For the modulation index h = 1/2, there are 2*2π / (1/2π) = 8 possibilities for the phase at the beginning of a code element. Because of the symmetry of the cosine function, the beginning phase states are reduced by half to four. For each beginning phase, the phase change can be generated according to the input symbols, so there are k_max = 4M theoretical variations, denote as Φ_k . Then the waveform at the time i of each received codeword can be correlated with:

\(\begin{gathered} l_{1}=\int_{(i-1) T}^{i T} \cos \Phi_{1}(t) * r(t) d t \\ l_{2}=\int_{(i-1) T}^{i T} \cos \Phi_{2}(t) * r(t) d t \\ \quad \vdots \\ l_{4 M-1}=\int_{(i-1) T}^{i T} \cos \Phi_{4 M-1}(t) * r(t) d t \\ l_{4 M}=\int_{(i-1) T}^{i T} \cos \Phi_{4 M}(t) * r(t) d t \end{gathered}\)

Where l_k = max(l₁, l₂, … , l_4M-1, l_4M).

We get the maximum value l_K . The cosΦ_k(t) t is chosen as the best reception for this waveform at the time i . We can also reduce the BER of the received signal with coherent demodulation by increasing the inherent memory L. In this case, Φ_k expressed as Φ^L_k , the optimal receive memory inherent L corresponds to the observation time of the code element while satisfying the condition of phase continuity. The numbers of all possibilities of length L, k_max is then increased to 4M^L .

\(\begin{gathered} {l^{L}}_{1}=\int_{(i-1) T}^{(i+L-1) T} \cos {\Phi^{L}}_{1}(t)^{*} r(t) d t \\ {l^{L}}_{2}=\int_{(i-1) T}^{(i+L-1) T} \cos {\Phi^{L}}_{2}(t)^{*} r(t) d t \\ \vdots \\ {l^{L}}_{4 M^{2}-1}=\int_{(i-1) T}^{(i+L-1) T} \cos {\Phi^{L}}_{4 M^{2}-1}(t)^{*} r(t) d t \\ {l^{L}}_{4 M^{2}}=\int_{(i-1) T}^{(i+L-1) T} \cos {\Phi^{L}}_{4 M^{2}}(t)^{*} r(t) d t \\ {l^{L}}_{K}=\max \left(l^{L}{ }_{1}, l^{L}{ }_{2}, \cdots, {l^{L}}_{4 M^{L}-1}, {l^{L}}_{4 M^{L}}\right) ; \end{gathered}\)        (14)

Fig. 11 shows the Receiving process. The diagram on the left shows the reception of the first code element. We correlate the received signal code elements numbered T1 to T(L) with \(\cos \Phi_{1}^{L}, \cos \Phi_{2}^{L}, \ldots, \cos \Phi_{4 M^{L}-1}^{L} \cos \Phi_{4 M^{L}}^{L}\) to get \(l_{1}^{L}, l_{2}^{L}, \ldots, l_{4 M^{L}-1}^{L}, l_{4 M^{L}}^{L}\) following (14). The waveform \(\cos \Phi_{K}^{L}(t)\) corresponding to the largest \(l^{K}_{L}\) is selected as the best reception for this length L signal. Choose only the first one as the decoder for this reception afterward.

Fig. 11. Length L best match schematic

The right part of Fig. 11 shows the reception of the second code element, the range of detection is shifted to the right by one code element in length, and the matching method is similar to that of the first code element and continues.

We set: h = 1/2 , SNR = E_b / (N₀log₂M)

Sampling frequency fs = 20000

Symbol duration T = 0.001

Number of input symbols n = 100000

The theoretical BER of M-PSK is as follows:

BPSK theoretical  BER: \(P_{e}=1-\frac{1}{2} \operatorname{erfc} \sqrt{r}\)

QPSK theoretical BER:\(P_{e}=1-\left[1-\frac{1}{2} \operatorname{erfc} \sqrt{r / 2}\right]^{2}\)

MPSK theoretical BER:\(P_{e} \approx \operatorname{erfc}\left(\sqrt{2 r} \sin \frac{\pi}{2 M}\right)\)

The BER curves of 2-ary MSK, 4-ary MSK, and 8-ary MSK are obtained by MATLAB simulation, respectively in Fig. 12. When L = 1, compared with theoretical M-ary PSK under each of the same M, 2-ary MSK OMLCD only approximates theoretical values. But for M > 2, the overall BER of the correlate reception can be better than the theoretical M-ary PSK, and the larger M, the better it is. This reflects the relative superiority of OMLCD compared to theoretical PSK when arys increase.

Fig. 12. Probability of symbol error for M-ary MSK, M=2,4,8,16Fig. 12. Probability of symbol error for M-ary MSK, M=2,4,8,16

To add optimal memory inherent L, there is no doubt that it will significantly improve BER. Those four MSK average lift of about 4 dB when L = 3 compared with L = 1 at high SNR, for the case of increasing L from 1 to 2 BER improve about 3 dB. Specifically, 4-ary MSK with L = 3 can promote theoretical QPSK BER by 3.1 dB, and the coherent 16-ary MSK with L = 3 can outperform theoretical 16-ary PSK signal reaching up to 5.2 dB at high SNR.

According to the experiment, the case of L = 2 greatly reduces the probability of error than the case of L = 1 and the calculation time is at an acceptable level. Taking into account the calculation time and BER, L = 2 is preferable. If the computing time is sufficient, the BER in the case of L = 3 is further reduced.

7. Conclusions

This paper has built a Markov Chain model to evaluate the autocorrelation function and get the M-ary MSK by Fourier transform for M-ary MSK. And numerical signal simulations verify the correctness of the results in those figures. It has been shown that the composition of each level M spectrum analyzes the characteristics of the spectrum. The fractional out-of-band power of each M-ary MSK that M = 2,4,8,16 has been compared, which shows that when BT_𝑏 > 0.71 , 4-ary MSK has significantly less out-of-band power compared to another communication system. We adopt OMLCD to receive the signal and draw the BER curves at different SNRs. For inherent memory L, we can conclude that when L changes from 1 to 2, the BER improves about 3dB. When L = 3, the BER improves about 4 dB. With the increase of MSK level, except binary MSK, the BER is 4 dB higher than the theoretical PSK level.

Appendix 1

In appendix 1, examples of MSK, 4-ary MSK and 8-ary MSK matrix are listed. Example: MSK:

\(P_{M=2}{ }^{1}=A P_{M=2}{ }^{2}=A^{2}\)

4-ary MSK:

\(P_{M=4}{ }^{1}=\left(\begin{array}{ll} A & A \\ B & B \end{array}\right) P_{M=4}{ }^{2}=\left(\begin{array}{cc} A^{2}+A B & A^{2}+A B \\ A B+B^{2} & A B+B^{2} \end{array}\right)\)

8-ary MSK:

\(P_{M=8}{ }^{1}=\left(\begin{array}{llll} A & A & A & A \\ B & B & B & B \\ A & A & A & A \\ B & B & B & B \end{array}\right) \quad P_{M=8}{ }^{2}=2 *\left(\begin{array}{cccc} A^{2}+A B & A^{2}+A B & A^{2}+A B & A^{2}+A B \\ A B+B^{2} & A B+B^{2} & A B+B^{2} & A B+B^{2} \\ A^{2}+A B & A^{2}+A B & A^{2}+A B & A^{2}+A B \\ A B+B^{2} & A B+B^{2} & A B+B^{2} & A B+B^{2} \end{array}\right)\)

Appendix 2

In appendix 2 shows the details of deriving the PSD function.

4.1 Case 1: Where 𝟎𝟎 ≤ 𝐭𝐭 ≤ 𝐓𝐓 and 𝐦𝐦 = 𝟎𝟎, \(S_{m=0}^{M}\)

\(\begin{aligned} S_{m=0}^{M}(f) &=\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau+\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}+1) \mathrm{s}_{j i}^{*}\left(\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau \\ &=2 \int_{0}^{T_{b}} R 1^{M}{ }_{m=0}(\tau) \cos 2 \pi f \tau d \tau+2 \int_{0}^{T_{b}} R 2^{M}{ }_{m=0}(\tau) \cos 2 \pi f \tau d \tau \\ &=2 \int_{0}^{T_{b}} R 1^{M}{ }_{m=0}(\tau) \cos 2 \pi f \tau d \tau+2 \int_{0}^{T_{b}} R 2_{-} \operatorname{sam} e_{m=0}^{M}(\tau) \cos 2 \pi f \tau d \tau+2 \int_{0}^{T_{b}} R 2_{-} d i f^{M}{ }_{m=0}(\tau) \cos 2 \pi f \tau d \tau \\ &=S 1^{M}{ }_{m=0}(f)+S 2_{-} \operatorname{same}_{m=0}^{M}(f)+S 2_{-} d i f^{M}{ }_{m=0}(f) \end{aligned}\)

\(S_{m=0}^{M}(f)\) is divided into \(S 2_{-} \text {same }_{m=0}^{M}(f), S 1_{m=0}^{M}(f)\) and \(S 2_{-} \operatorname{dif}_{m=0}^{M}(f)\) for the summation calculation.

\(\begin{aligned} S 1^{M}{ }_{m=0}^{M}(f) &=2 \int_{0}^{T_{b}} \frac{1}{2}\left(1-\frac{\tau}{T_{b}}\right) \cos \frac{(2 n-1) \pi \tau}{2 T_{b}} \cos 2 \pi f \tau d \tau \\ &=\int_{0}^{T_{b}} \cos \frac{n \pi \tau}{2 T_{b}} \cos 2 \pi f \tau d \tau-\int_{0}^{T_{b}} \frac{\tau}{T_{b}} \cos \frac{n \pi \tau}{2 T_{b}} \cos 2 \pi f \tau d \tau \\ &=T\left\{\frac{\sin [\varphi(n)+\omega(f)]}{2[\varphi(n)+\omega(f)]}+\frac{\sin [\varphi(n)-\omega(f)]}{2[\varphi(n)-\omega(f)]}-\frac{\cos [\varphi(n)-\omega(f)]+[\varphi(n)-\omega(f)] \sin [\varphi(n)-\omega(f)]-1}{2[\varphi(n)-\omega(f)]^{2}}\right.\\ &\left.-\frac{\cos [\varphi(n)+\omega(f)]+[\varphi(n)+\omega(f)] \sin [\varphi(n)+\omega(f)]-1}{2[\varphi(n)+\omega(f)]^{2}}\right\} \end{aligned}\)

\(\begin{aligned} S 2_{-} \operatorname{sam} e_{m=0}^{M}(f) &=2 \int_{0}^{T_{b}} \frac{1}{2} \frac{\tau^{\prime}}{T_{b}} \sin \frac{\pi\left(T_{b}-\tau^{\prime}\right) n}{2 T_{b}} \cos 2 \pi f \tau d \tau+2 \int_{0}^{T_{b}} \frac{\sin \left(\frac{1}{2} \pi n\right) \sin \frac{\pi \tau^{\prime} n}{2 T_{b}}}{\pi n} \cos 2 \pi f \tau d \tau \\ &=4 T \frac{\varphi(n)\left[\cos \omega(f) * S_{2}+8 \omega(f) \sin \omega(f)\right]-\sin \varphi(n) * S_{1}(n)}{S_{2}^{2}(n)} \\ &-4 T \frac{\sin \varphi(n)[\omega(f) \sin \omega(f) \sin \varphi(n)+\varphi(n) \cos \omega(f) \cos \varphi(n)-\varphi(n)]}{\varphi(n) S_{2}(n)} \end{aligned}\)

\(\begin{aligned} S 2_{-} d i f^{M}{ }_{m=0}(f) &=4 \int_{0}^{T_{b}} \frac{n\left[\sin \left(\frac{1}{2} \pi m\right) \sin \frac{\pi \tau^{\prime} n}{2 T_{b}}\right]}{\pi\left[m^{2}-n^{2}\right]} \cos 2 \pi f \tau d \tau-4 \int_{0}^{T_{b}} \frac{m\left[\cos \frac{\pi\left(T_{b}-\tau^{\prime}\right) m}{2 T}\right]}{\pi\left[m^{2}-n^{2}\right]} \cos 2 \pi f \tau d \tau \\ &=8 T \frac{\varphi(n) \sin \varphi(m) *[\omega(f) \sin \omega(f) \sin \varphi(n)+\varphi(n) \cos \omega(f) \cos \varphi(n)-\varphi(n)]}{[\varphi(n)-\varphi(m)][\varphi(n)+\varphi(m)] S_{2}(n)} \\ &+8 T \varphi(m) \frac{\varphi(m) \sin \varphi(m)-\omega(f) \sin \omega(f)}{[\varphi(n)-\varphi(m)][\varphi(n)+\varphi(m)] S_{2}(m)} \end{aligned}\)

4.2 Case 2: Where 𝐓𝐓 ≤ 𝐭𝐭 ≤ 𝟐𝟐𝟐𝟐 and 𝐦𝐦 = 𝟏𝟏, \(S_{m=1}^{M}\)

\(\begin{aligned} S_{m=1}^{M} &=\sum_{i=1}^{N} \sum_{j=1}^{N} \int_{-\infty}^{\infty} 2 P(\mathrm{i})\left[\mathrm{p}(\mathrm{j} \mid \mathrm{i}, \mathrm{m}) \mathrm{s}_{i j}\left(\mathrm{~T}-\tau^{\prime}\right)\right] \cos 2 \pi f \tau d \tau \\ &=2 \int_{T_{b}}^{2 T_{b}} R 1_{-} s a m e_{m=1}^{M}(\tau) \cos 2 \pi f \tau d \tau+2 \int_{T_{b}}^{2 T_{b}} R 1_{-} d i f_{m=1}^{M}(\tau) \cos 2 \pi f \tau d \tau \\ &=S 1_{-} s a m e_{m=1}^{M}(f)+S 1_{-} d i f^{M}{ }_{m=1}(f) \end{aligned}\)

\(S_{m=1}^{M}\) is divided into \(S 1_{-} \operatorname{same}_{m=1}^{M}(f)\) and \(S_{1_{-}} d i f_{m=1}^{M}(f)\) for the summation calculation.

\(S 1_{-} \operatorname{same}_{m=1}^{M}(f)=2 \int_{T_{b}}^{2 T_{b}} \sin \frac{\pi\left(\tau-T_{b}\right) n}{2 T_{b}} \cos 2 \pi f \tau d \tau\)

\(\begin{aligned} &-\int_{T_{b}}^{2 T_{b}} \frac{\tau}{T_{b}} \sin \frac{\pi\left(\tau-T_{b}\right) n}{2 T_{b}} \cos 2 \pi f \tau d \tau-2 \int_{\tau_{b}}^{2 T_{b}} \frac{\sin \left(\frac{1}{2} \pi n\right) \sin \frac{\pi\left(2 T_{b}-\tau\right) n}{2 T_{b}}}{\pi n} \cos 2 \pi f \tau d \tau \\ &=-2 T \frac{2[\varphi(n)+\omega(f)] \cos [\varphi(n)-2 \omega(f)]+2[\varphi(n)-\omega(f)] \cos [\varphi(n)+2 \omega(f)]-4 \varphi(n) \cos \omega(f)}{S_{2}(n)} \\ &+2 T\left\{\frac{2[\varphi(n)-\omega(f)] \cos [\varphi(n)-2 \omega(f)]-\sin [\varphi(n)-2 \omega(f)]}{4[\varphi(n)-\omega(f)]^{2}}\right. \\ &\left.+\frac{[\varphi(n)+\omega(f)] \cos [\varphi(n)+2 \omega(f)]-\sin [\varphi(n)+2 \omega(f)]}{4[\varphi(n)+\omega(f)]^{2}}-\frac{2 \varphi(n)\left[S_{2} \cos \omega(f)+8 \omega(f) \sin \omega(f)\right]}{S_{2}^{2}(n)}\right\} \\ &+2 T \sin \varphi(n)\left[\frac{\cos [\varphi(n)-\omega(f)]+\cos [\varphi(n)+\omega(f)]-2 \cos 2 \omega(f)}{S_{2}(n)}-\frac{2 \omega(f) \sin \omega(f) \sin \varphi(n)}{\varphi(n) S_{2}(n)}\right] \end{aligned}\)

\(\begin{aligned} S 1_{-} d i f^{M}{ }_{m=1}(f) &=2 \int_{T_{b}}^{2 T_{b}} \frac{[n-m] \cos \frac{\pi\left[T_{b} m+(2 T-\tau) n\right]}{2 T_{b}}}{\pi\left[m^{2}-n^{2}\right]} \cos 2 \pi f \tau d \tau \\ &-2 \int_{T_{b}}^{2 T_{b}} \frac{2 n \cos \frac{\pi\left(\tau-T_{b}\right) m}{2 T_{b}}}{\pi\left[m^{2}-n^{2}\right]} \cos 2 \pi f \tau d \tau+2 \int_{T_{b}}^{2 T_{b}} \frac{[m+n] \cos \frac{\pi\left(T_{b} m+\left(\tau-2 T_{b}\right) n\right)}{2 T_{b}}}{\pi\left[m^{2}-n^{2}\right]} \cos 2 \pi f \tau d \tau \\ &=\frac{T}{2}\left\{\frac{\sin [\varphi(m)+2 \omega(f)]-\sin [\varphi(n)+\varphi(m)+\omega(f)]}{[\varphi(n)+\varphi(m)][\varphi(n)-\omega(f)]}-\frac{\sin [2 \omega(f)-\varphi(m)]+\sin [\varphi(n)+\varphi(m)-\omega(f)]}{[\varphi(n)+\varphi(m)][\varphi(n)+\omega(f)]}\right\} \\ &-T \varphi(n) \frac{[\varphi(m)-\omega(f)] \sin [\varphi(m)+2 \omega(f)]+[\varphi(m)+\omega(f)] \sin [\varphi(m)-2 \omega(f)]+2 \omega(f) \sin \omega(f)}{\left[\varphi^{2}(m)-\varphi^{2}(n)\right]\left[\varphi^{2}(m)-\omega^{2}(f)\right]} \\ &+T\left\{\frac{\sin [\varphi(m)+2 \omega(f)]-\sin [\varphi(m)-\varphi(n)+\omega(f)]}{2[\varphi(m)-\varphi(n)][\varphi(n)+\omega(f)]}-\frac{\sin [\varphi(m)-\varphi(n)-\omega(f)]-\sin [\varphi(m)-2 \omega(f)]}{2[\varphi(m)-\varphi(n)][\varphi(n)-\omega(f)]}\right\} \end{aligned}\)