Maximum Likelihood Sequence Detector for GMSK

1. Abstract

This report is about the maximum likelihood sequence detection (MLSD) of GMSK scheme. The receiver is done both in coherent and noncoherent means without using limiter discriminator structure or PAM representation of GMSK.

You can find sample MATLAB codes here.

2. Introduction

GMSK waveform can be modeled as both PAM and as a Trellis. Such behavior is beneficial since it enables different receiver structures [1, 2, 3]. A separate serial receiver that exploits the PAM representation of GMSK was examined previously. Both the coherent and non-coherent serial receivers were realized according to Laurent's PAM representation of GMSK and Kaleh's serial receiver [4]. In this short report, MLSD is employed as receiver due to GMSK's Trellis model.

In the following section, necessary mathematical derivations are done in order to write the GMSK waveform's phase as a state machine, ie. some outputs are generated according to some start state. After necessary outputs are generated, the state changes according to a Trellis that is specific to GMSK waveform. In fact, this Trellis structure is not specific to GMSK only but applies to all Continuous Phase Modulation (CPM) schemes.

In the third chapter the MLSD receiver is set so that this Trellis can be solved via Viterbi Algorithm. Unlike regular Viterbi algorithm that is used to solve convolutional coded sequences, this one uses a soft input structure to solve the GMSK Trellis. The regular MLSD receiver is a coherent one. The non-coherent MSLD (NSD) modification is also presented here.

Finally, BER curves are presented with further comments.

3. CPM Waveform as a Trellis

GMSK is a form of CPM. If CPM can be written as a Trellis, GMSK can be as well. CPM waveform has the following structure [1]:

$\displaystyle s(t, \underline{\alpha}) = \sqrt{\frac{E}{T}} exp\{j\psi (t, \underline{\alpha})\},$(1)

where $E$ is the symbol energy and $T$ is the symbol duration and $\underline{\alpha}$ is a sequence of bipolar NRZ symbols for GMSK. $\psi (\cdot)$ is the phase of the signal which also carries the information. Phase function has the following form [1, 5, 2, 6]:

$\displaystyle \psi (t, \underline{\alpha}) = 2\pi h\sum_{i=-\infty}^{n}\alpha_i q(t-iT), \quad nT \leq t < (n+1)T,$(2)

where $h$ is the modulation index and $h = 0.5$ for GMSK. $q(\cdot)$ is the phase shaping filter and is a Gaussian curve for GMSK case. $q(\cdot)$ is defined as follows [1]:

$\displaystyle q(t) = \begin{cases} 0, & t < 0, \\ \int_{0}^{t}f(\tau)d\tau, & 0 \leq t < LT,\\ 1/2, & t \geq LT \end{cases}$(3)

where $f(\cdot)$ is called the frequency pulse. Here, the integral of the frequency pulse might add up to more than $1/2$ when $t \ge LT$, then $q(t)$ is basically scaled in amplitude. For GMSK, this has Gaussian shape and is calculated via [1]:

$\displaystyle f(t) = \frac{1}{2T}\bigg[Q(2\pi B \frac{t-\frac{T}{2}}{\sqrt{ln(2)}}) - Q(2\pi B \frac{t+\frac{T}{2}}{\sqrt{ln(2)}})\bigg],$(4)

where $B = BT/T$, with $BT$ the so called time bandwidth product. $Q(\cdot)$ is the famous $Q$ function.

One important note with Eq. (2) before going further is that, it is defined for the time interval $nT \leq t < (n+1)T$ when the $n^{th}$ symbol arrives and the phase waveform in Eq. (2) has a transient shape only in $LT$ duration, which limits the transient curve to $L$ consecutive symbols within the $\underline{\alpha}$. From another perspective, this is due to the fact that $f(\cdot)$ is only defined in $(0, LT_{sym})$.

4. Partition of $\psi(t, \underline{\alpha})$

Phase shaping filter $q(t)$ has variations only within $t\in [0, LT]$. Looking at the Eq. (3), when $t < 0$, $q(t) =0$, which means $\alpha_i, i > n$ will have no effect on the value of $\psi(t, \underline{\alpha})$, due to $nT \leq t < (n+1)T$. On the contrary, $\alpha_i, i \leq n-L$ will have $\sum_{i=0}^{n-L} \alpha_i\frac{2\pi h}{2}$ cumulative effect on the value of $\psi(t, \underline{\alpha})$. Here, $i$ starts from 0 instead of $-\infty$ since we cannot have infinite amount of symbols. This comes from the way $\psi(t, \underline{\alpha})$ is defined. Lets now partition $\psi(t, \underline{\alpha})$ to see aforementioned effect better.

Say, $L = 3$, $N = 5$, symbols so, $i = 0,1,2,...,5$ and $n = 4$, ie. $n=4^{th}$ symbol just arrived and value of $\psi (t, \underline{\alpha})$ at $\quad 4T \leq t < 5T$ is of concern:

$\displaystyle \psi (t, \underline{\alpha}) = 2\pi h \sum_{i=0}^{4} \alpha_i q(t - iT), \quad 4T \leq t < 5T,$(5)

$\displaystyle \begin{aligned} \psi (t, \underline{\alpha}) = 2\pi h\bigg[ & \alpha_0 q(t) + \alpha_1 q(t - T) + \alpha_2 q(t - 2T) \\ & + \alpha_3 q(t - 3T) + \alpha_4 q(t - 4T)\bigg], \quad 4T \leq t < 5T, \end{aligned}$(6)

$\displaystyle \begin{aligned} \psi (t, \underline{\alpha}) = 2\pi h\bigg[ & \alpha_0 q(4T \leq t < 5T) + \alpha_1 q(3T \leq t < 4T) \\ & + \alpha_2 q(2T \leq t < 3T) + \alpha_3 q(T \leq t < 2T) + \alpha_4 q(0 \leq t < T)\bigg], \end{aligned}$(7)

First of all, $q(4T \leq t < 5T) = q(3T \leq t < 4T) = 1/2$ since $L = 3$ and (3) dictates $q(t) = 1/2$ for $t \geq LT$. This is the direct effect of the cumulative $\sum_{i=0}^{n-L} \alpha_i\frac{2\pi h}{2}$. Since $L = 3$ and $n = 4$ is chosen. Thus for this case, when $n = 4^{th}$ symbol arrives, $\alpha_0$ and $\alpha_1$ has a cumulative and constant effect on the value of $\psi (t, \underline{\alpha})$.

Furthermore, $i=5$ is said to have a similar effect, since the time interval of concern is $\quad 4T \leq t < 5T$, $q(\cdot)$ belonging to the $i=5^{th}$ symbol would be, $q(-T \leq t < 0) = 0$ due to the way $q(\cdot)$ is defined in Eq. (3). Thus, it is true that there are 6 symbols that have an impact on $\psi (t, \underline{\alpha})$. When $n=4^{th}$ symbol arrives, $\alpha_i, i \in [0, n-L]$ have constant cumulative effect whereas $\alpha_i, i \in [n+1, N]$ have no effect, since they are yet to be fed into the system.

Thus, $\psi (t, \underline{\alpha})$ is partitioned into two parts [2, 6]:

$\displaystyle \psi (t, \underline{\alpha}) = \theta(t, \underline{\alpha_n}) + \theta_{n-L},$(8)

where,

$\displaystyle \theta(t, \underline{\alpha_n}) = 2\pi h \sum_{i=n-L+1}^{n} \alpha_i q(t-iT),$(9)

and

$\displaystyle \theta_{n-L} = \sum_{i=-\infty}^{n-L} \alpha_i\frac{2\pi h}{2} = \pi h\sum_{i=-\infty}^{n-L} \alpha_i,$(10)

Here, $nT \leq t < (n+1)T$.

Eq. (9) is the transient phase that defines the phase curve and is dependent on the $L-1$ previous symbols that come before $n^{th}$ symbol. This resembles the idea of a start state that consists of $L-1$ previous symbols, an input at time $n$ and an output is calculated, $\psi (t, \underline{\alpha})$. Obviously, when the $n^{th}$ symbol comes in, it causes a transition from the start state to an end state. Which awfully sounds like a Trellis.

However, there is also the cumulative value of $\theta_{n-L}$ in Eq. (10) that needs to be taken care. When $n^{th}$ symbol arrive, $\alpha_i, i \in [n-L+1, n]$ is used to calculate the transient state as per Eq. (9). The previous $\alpha_i, i \in [0, n-L]$ are used to calculate the cumulative state up until to the $n-L+1^{th}$ time. Which is some sort of a separate so called phase state, even though its not actually a state.

In order to visualize this, a small example is as follows; lets assume that $BT = 0.3$, $L = 3$, $sps = 20$. For this case, there are $2^{L-1} = 4$ transient states, excluding the aforementioned cumulative state right now. The states, inputs and corresponding outputs are then:

Figure 1. States, inputs and corresponding transient phase curves.These curves are added side by side with the addition of cumulative value calculated via Eq. (10). Let us now assume a sequence of symbols flowing:

$\displaystyle \underline{\alpha} = [\underset{\underset{n=0}{\downarrow}}{0} 1\, 1\, 0\, 1].$(11)

Here, the starting state of the GMSK scheme is $[0\, 0]$, since $L=3$. First symbol arrives at $n=0$.

1. First the symbols are converted into bipolar-NRZ such that, $0\rightarrow -1, \quad 1 \rightarrow 1$,
2. A sliding window of length $L$ travels along the data, hence at time $n=0$, the sub-sequence of concern is $\underline{\alpha}_{sub} = [-1\, -1\, -1]$,
3. The cumulative phase is calculated as per Eq. (10), for this case, at symbol $n=0$, cumulative phase is $0$ since there are no symbols $n-L = -3$ yet,
4. The necessary transient phase of this sub-sequence is chosen from Fig. 1,
5. Cumulative phase is added onto the transient phase to calculate the total phase for the time interval of concern.

These steps are done over and over until the data in the sequence finishes. Thus for the symbol sequence of concern, the phase values will be as follows:

Figure 2. Phase curve of the aforementioned $\underline{\alpha}$ sequence.
In Fig. 2, at time $n=1$, the sub-sequence is, $\underline{\alpha}_{sub} = [-1\, -1\, 1]$ and cumulative phase that is left over from time $n=0$ is $-\frac{\pi}{2}$. Thus the curve that corresponds to the sub-sequence is shifted by this value, which can be seen from sample number 21, beginning of the transient phase of the sub-sequence. At time $n=2$, a similar behavior is apparent, $\underline{\alpha}_{sub} = [-1\, 1\, 1]$ and the cumulative phase that is left over from time $n=0$ and $n=1$ is now $-\pi$, which shifts the transient curve of sub-sequence by this amount. This can be observed at sample 41 and so forth.

Thus the cumulative phase affects the Trellis of CPM scheme and enlarges it.

5. Coherent MLSD Receiver

Received signal is modeled as:

$\displaystyle r(t, \underline{\alpha}) = s(t, \underline{\alpha}) + n(t),$(12)

where, $n(t)$ is assumed to be AWGN.

The MLSD receiver is intended for finding an $\underline{\hat{\alpha}}$ sequence such that it is closest (in a sense) to the received noisy signal:

$\displaystyle \Lambda(\underline{\hat{\alpha}}) = -\int_{-\infty}^{\infty} |r(t) - s(t, \underline{\hat{\alpha}})|^2 dt$(13)

An $\underline{\hat{\alpha}}$ sequence that maximizes the Eq. (13) metric, is the most probable sequence that is transmitted. Maximizing the above sequence for CPM schemes is equal to maximizing correlation [1]:

$\displaystyle \lambda(\underline{\hat{\alpha}}) = Re\bigg[\int_{-\infty}^{\infty} r(t)s^*(t, \underline{\hat{\alpha}}) dt\bigg],$(14)

Due to the aforementioned Trellis-like structure of CPM signals, the correlation in Eq. (14) can be calculated via Viterbi Algorithm (VA).

Using Eq. (9) and (10), an L-tuple start state can be defined at time $n$ when $\alpha_n$ arrives,

$\displaystyle S_n = (\theta_{n-L}, \alpha_{n-L+1}, ..., \alpha_{n-2}, \alpha_{n-1}),$(15)

As previously shown, this start state depends on the $L-1$ $\alpha$ symbols that came just before $n^{th}$ symbol and the cumulative phase that accumulated from the beginning of the data transmission to the $(n-L)^{th}$ symbol. The corresponding end state after $n^{th}$ symbol is then [2, 6]:

$\displaystyle E_n = (\theta_{n-L+1}, \alpha_{n-L+2}, ..., \alpha_{n-1}, \alpha_{n}),$(16)

The transition from $S_n$ to $E_n$ is done by the following $L+1$-tuple:

$\displaystyle \sigma_n = (\theta_{n-L}, \alpha_{n-L+1}, ..., \alpha_{n-1}, \alpha_{n}).$(17)

Combination of these definitions enable a recursive metric update, the metric that accumulated up to $n^{th}$ symbol as follows:

$\displaystyle \lambda_{n+1}(\hat{E}_n) = \lambda_n(\hat{S}_n) + Re\bigg[e^{-j\hat{\theta}_{n-L}} z_n(\underline{\hat{\alpha}}_n)\bigg].$(18)

Here, $e^{-j\theta_{n-L}}$is the cumulative phase up until to $n^{th}$ symbol for the hypothetical sequence $\underline{\hat{\alpha}}_n$ and $z_n(\cdot)$ is the sampled matched filter output:

$\displaystyle z_n(\underline{\hat{\alpha}}_n) = \int_{nT}^{(n+1)T} r(t)s^*(t, \underline{\hat{\alpha}}_n)dt.$(19)

Writing the metric update as in Eq. (18) enables the opportunity to decrease the number of states from $\rho M^{L-1}$ to $M^{L-1}$ because of the way $e^{-j\theta_{n-L}}$ is calculated. Starting from the first time instant $n=0$, if the starting state is known a-priori and if $e^{-j\theta_{n-L}}$ is updated when the winner branch up to that time instant, then $e^{-j\theta_{n-L}}$ should have $M^{L-1}$ different values for each time step. Thus each state would have their own "winner" phases that accumulate as the maximum of branch metrics are chosen. Thus, it is not necessary to search for the best of all $\rho$ possible $e^{-j\theta_{n-L}}$ values, but calculating the necessary winner previous state is enough,

$\displaystyle \hat{\theta}_{n-L+1}(\hat{E}_n) = \hat{\theta}_{n-L}(\hat{S}_n) + \pi h\hat{\alpha}_{n-L+1}.$(20)

This approach decreases the total number of states to $M^{L-1}$ for the Viterbi Algorithm.

6. Noncoherent Receiver

In noncoherent case received signal is modeled as:

$\displaystyle r(t, \underline{\alpha}) = s(t, \underline{\alpha})e^{j\phi(t)} + n(t),$(21)

where $\phi(t)$ is some time varying phase. In slow fading channels, it can be assumed that $\phi(t) = \phi \sim U[0, 2\pi]$ during the interval in which data is transmitted, which is the case here as well. In such circumstances, the receiver structure also changes. Phase locked loops are common solutions to the random additional phase but they are also expensive.

Thus a noncoherent receiver which has to compensate for the random phase of the channel has to be developed. In the MLSD receiver, the recursive metric calculation in Eq. (14) is updated as follows:

$\displaystyle \lambda_{n+1}(\hat{E}_n) = \lambda_n(\hat{S}_n) + \bigg[Q^*_n(\hat{S}_n)e^{-j\hat{\theta}_{n-L}} z_n(\underline{\hat{\alpha}}_n)\bigg],$(22)

where $Q_n(\cdot)$ is the so called phase reference and is updated via a gradient descent like approach [7];

$\displaystyle Q_{n+1}(\hat{E}_n) = aQ_n(\hat{S}_n) + (1-a)\bigg[e^{-j\hat{\theta}_{n-L}}z_n(\hat{\alpha}_n)\bigg].$(23)

Here, $z_n(\hat{\alpha}_n)$ is calculated via Eq. (19), the matched filter output and $e^{-j\hat{\theta}_{n-L}}$ is the same phase update in the coherent receiver. The $a$ parameter is the so called forgetting factor and has the same purpose as in gradient descent approach. As $a \rightarrow 1$, the system becomes more open to the side effects of the random phase shift, thus becomes "more" coherent. In fact $a = 1$ is exactly the coherent receiver. As $a \rightarrow 0$ the opposite happens. Hence, $a$ controls, how much system depends on the current phase as opposed to the previous phase.

7. Example Problem for the Coherent Case

In this problem, it is very important to understand that, if the seed of the random noise is not given same as here, it may not be possible to reach the same results.

Lets assume a sample problem with GMSK, $BT = 0.3$, $L = 3$, $N = 4$ and starting state of $[0\, 0]$, and transmitted symbols $\underline{\alpha}$ sequence is as follows:

$\displaystyle \underline{\alpha} = [\underset{\underset{n=0}{\downarrow}}{0} 1\, 1\, 0\, 1].$(24)

Corresponding phase function of this symbol stream is in Fig. 2 as previously shown. The reduced state Trellis diagram has $2^{L-1} = 4$ states. In this example, system is simulated in MATLAB according to the equations in this report. The SNR is, $\frac{E_b}{N_0} = 8$ dB and oversampling ratio $sps = 20$. Noise seed is rng(40) and noise is added via awgn(.) function with 'measured' flag. Thus the noise code snippet becomes:

r = awgn(s, 8-10*log10(sps), 'measured');
The corresponding Trellis is in Fig. 3.

Figure 3. Trellis of the example problem.
In Fig. 3, the winner branch at the end of the process is marked with red. The solid lines correspond to hypothetical input of $-1$ and dashed lines correspond to a hypothetical input $1$. Numbers on the lines correspond to the cumulative metrics as calculated in Eq. (18). Cumulative phases that correspond to the phase update as per Eq. (20) can be seen under the state nodes.

Going backwards in the trellis, the winner path is in Fig. 4. This winner path corresponds to the data sequence of $\underline{\hat{\alpha}} = [0\, 1\, 1\, 0\, 1]$ which is in fact the transmitted data sequence.

Figure 4. Winner path.

8. Conclusion

The corresponding BER curves for various other GMSK receivers and some BPSK receivers is shown in Fig. 5. In the figure, the MSLD based receivers prove to be best among other types when GMSK is of concern even at low BT values such as $BT = 0.3$. Especially for the non-coherent case the performance of MSLD based receivers are better.

Figure 5. BER results of several different receivers, $BT = 0.3$.

References

[1] John B Anderson, Tor Aulin, and Carl-Erik Sundberg. Digital phase modulation. Springer Science & Business Media, 2013.
[2] Erik Samuel Perrins. Reduced complexity detection methods for continuous phase modulation. 2005.
[3] P¨ar Moqvist et al. Serially concatenated systems: An iterative decoding approach with application to continuous phase modulation. Lic. Eng. thesis, Chalmers University of Technology, Gothenburg, Sweden, 1999.
[4] Musa Tunç Arslan. PAM Representation of GMSK and Serial Receivers. Technical report, Meteksan Defence Inc., 2018.
[5] Subhashini Gupta, Vikas Bhatia, and LC Mangal. An efficient fpga implementation of gmsk (bt= 0.3) transceiver with non coherent sequence detection for tactical v/uhf waveforms. In Communications (NCC), 2012 National Conference on, pages 1–5. IEEE, 2012.
[6] Afzal Syed. Comparison of noncoherent detectors for soqpsk and gmsk in phase noise channels. International Foundation for Telemetering, 2007.
[7] Lutz Lampe, Robert Schober, and Mani Jain. Noncoherent sequence detection receiver for bluetooth systems. IEEE Journal on Selected Areas in Communications, 23(9):1718–1727, 2005.