I - In the previous two posts we discussed MSK performance in an AWGN channel, first presenting the MATLAB/OCTAVE Code for one sample per symbol case [Post 1], and then extending it to the more general case of multiple samples per symbol [Post 2]. This helps us visualize the underlying beauty of Continuous Phase Modulation (CPM) which reduces out of band energy and consequently lowers Adjacent Channel Interference (ACI). We also briefly touched upon the case of MSK in Rayleigh fading, but did not go into the details. So here we take a deeper dive.

II - On the face of it its quite simple. Just as we add an Independent and Identically distributed (IID) Complex Gaussian Random Variable to the signal we can do the same for the fading channel (the only difference being that the channel is multiplicative and noise is additive). So the received signal would look something like this y=h.*s+n where 'h' and 'n' are defined as given below ('Ns' is the samples/symbol and 'Nbits' is the total number of bits or symbols):

n=sigma*(randn(1,Ns*Nbits)+i*randn(1,Ns*Nbits))

h=(1/sqrt(2))*(randn(1,Ns*Nbits)+i*randn(1,Ns*Nbits))

III - But there is a caveat. If lets say we have four samples per symbol and we multiple the signal vector 's' with channel coefficient vector 'h', with 'h' having four IID values then we are providing a diversity advantage and we will be experiencing a gain in the performance like we do in Maximal Ratio Combining (assuming we multiply the received signal with the conjugate of the channel before correlation with the basis function). Indeed this is validated through Monte Carlo simulation where the BER continuously improves as the number of samples per symbol increases.

IV - So what is the solution? One solution is to assume that the channel is Quasi Static i.e. it is constant over the duration of a single symbol and then changes in the next symbol (without knowledge of prior value or in other words without memory). In our simulation we considered number of samples per symbol to be four and the channel to be static over 100 symbols (or bits assuming 1 bit/symbol). In our simulations we found the simulated BER to exactly match the theoretical BER in Rayleigh channel which is given as 0.5*(1-sqrt(EbNo/(EbNo+1))).

## OCTAVE Code for MSK in Rayleigh Fading

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MINIMUM SHIFT KEYING %
% GENERALIZED MATRIX BASED %
% BER OF MSK IN RAYLEIGH FADING %
% www.raymaps.com %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
Ns=4;
Nbits=1e6;
Nstat=100;
Eb=Ns;
EbNodB=10;
% MODULATION
bits_in=round(rand(1,Nbits));
differential_bits=abs(diff([0, bits_in]));
bipolar_symbols=differential_bits*2-1;
symbols_oversampled=(ones(Ns,1))*bipolar_symbols;
symbols_oversampled=(1/Ns)*transpose(symbols_oversampled(:));
cumulative_phase=(pi/2)*cumsum(symbols_oversampled);
tx_signal=exp(i*(cumulative_phase));
% CHANNEL (NOISE AND FADING)
EbNo=10^(EbNodB/10);
sigma=sqrt(Eb/(2*EbNo));
AWGN_noise=randn(1,Ns*Nbits)+i*randn(1,Ns*Nbits);
ray_chan=(1/sqrt(2))*(randn(1,Nbits/Nstat)+i*randn(1,Nbits/Nstat));
ray_chan_mat=ones(Nstat*Ns,1)*ray_chan;
ray_chan_vec=transpose(ray_chan_mat(:));
rx_signal=ray_chan_vec.*tx_signal+sigma*AWGN_noise;
% DEMODULATION (EQUALIZATION, CORRELATION, DECISION)
rx_signal=rx_signal.*conj(ray_chan_vec);
sin_waveform=sin(pi/(2*Ns):pi/(2*Ns):pi);
rx_real=[real(rx_signal(Ns+1:end)), zeros(1,Ns)];
rx_imag=[imag(rx_signal(1:end))];
rx_I_metric=sin_waveform*reshape(rx_real, 2*Ns, Nbits/2);
rx_Q_metric=sin_waveform*reshape(rx_imag, 2*Ns, Nbits/2);
n=1:Nbits/2;
rx_I_metric=(-1).^(n+1).*rx_I_metric;
rx_Q_metric=(-1).^(n+1).*rx_Q_metric;
rx_combined=[rx_Q_metric; rx_I_metric];
rx_combined=transpose(rx_combined(:));
demodulated_symbols=sign(rx_combined);
bits_out=demodulated_symbols*0.5+0.5;
%BER CALCULATION
ber_theoretical=0.5*(1-sqrt(EbNo/(EbNo+1)))
ber_simulated=sum(bits_out!=bits_in)/Nbits
```

## Simulation Results

## Note:

- The code above also works for Ns=1 i.e. for one sample per symbol. Which as discussed earlier is just like doing BPSK simulation and does not require the channel to be Quasi Static.
- It must also be noted that same results can be obtained by using a different method of equalization, division by the channel instead of multiplication with the complex conjugate.
- Please note that we perform differential encoding before modulation at the transmitter side. This results in a very elegant coherent receiver architecture with BER performance same as BPSK [1].

[1] Digital Communication over Fading Channels by Marvin K. Simon, Mohamed-Slim Alouini

Can you please provide the MATLAB code for GMSK Bit Error Rate in Rayleigh Fading.

Thank you very much.

Good question Rahul,

This is quite simple. You just need to convolve the symbol stream with the Impulse Response of the Gaussian pulse shaping filter before passing it through the fading channel. Remember to use “same” as output length of the pulse shaping filter. The MATLAB function is named “conv”. Please let me know if something is still not clear.

YA

In your code, the noise standard deviation (sigma) is multiplied with sqrt(Eb), with Eb=samples per symbol.

But usually, the AWGN noise we considered is having 0 mean and sigma=sqrt(1/(2*10^(EbNodB/10)), Eb =1, i.e., unit energy.

Can you please explain why you are doing this?

EbNo=Eb/2*sigma^2

2*sigma^2=Eb/EbNo

sigma^2=Eb/(2*EbNo)

sigma=sqrt(Eb/(2*EbNo))

This is for the general case. If Eb=1 then noise standard deviation is equal to:

sigma=sqrt(1/(2*EbNo))

Remember, energy is nothing but square of the distance from the origin (typically we keep this as 1). When there are multiple samples per bit, the energy per bit (or energy per symbol in this case) is just equal to the number of samples per bit (or samples per symbol in this case).

But in your code also after oversampling the power of each bit is still unity because you have normalized it by Ns.

symbols_oversampled=(1/Ns)*transpose(symbols_oversampled(:));

cumulative_phase=(pi/2)*cumsum(symbols_oversampled);

tx_signal=exp(i*(cumulative_phase));

The normalization is required so that the phase of the exponential advances by pi/2 during each symbol period. But the energy of the exponential at each sampling instant is still unity. It is of the form a*exp(j*theta) with a=1. Look at the following line of code:

tx_signal=exp(i*(cumulative_phase));

Hope it makes sense now!

YA

Okay, make sense now. Thank you so much.

It means if I take a=1/sqrt(Ns), it is equivalent to what you did with the noise sigma.

Yes that would make Eb=1 and noise standard deviation would be reduced to sigma=sqrt(1/(2*EbNo)).