Before we delve deep into Minimum Shift Keying (MSK) and its performance in presence of co-channel interference the reader is advised to look at the following posts.

Post 1 – MSK BER performance in AWGN and flat fading environment when viewed as extension of BPSK

Post 2 – MSK Power Spectral Density and its BER performance in AWGN when viewed as a CPM

Post 3 – MSK BER Performance in AWGN and flat fading environment when viewed as a CPM

Co-channel interference is a phenomenon widely encountered in wireless communication systems and the main reason for that is frequency reuse, which allows the same frequency band to be used over and over again in geographically non-contiguous areas. GSM and other wireless communication systems, using MSK modulation, suffer from the same problem. This has been widely studied in the literature and interference rejection techniques have been proposed. The worst case is one where the power of both the signals (wanted signal and interference) is almost the same and there is no frequency or phase offset.

I – In the previous post we presented the mathematical model and code for BER calculation of a popular modulation scheme called MSK. However in the code we shared, we only considered one sample per symbol, which makes MSK look like BPSK. While BPSK symbols fall on the real axis, MSK symbols alternate between real and imaginary axes, progressing by π/2 phase during each symbol period. MSK signal thus has memory and this can help in demodulation using advanced techniques such as Viterbi Algorithm.

I - Minimum Shift Keying (MSK) is a type of Continuous Phase Modulation (CPM) that has been used in many wireless communication systems. To be more precise it is Continuous Phase Frequency Shift Keying (CPFSK) with two frequencies f1 and f2. The frequency separation between the two tones is the minimum allowable while maintaining orthogonality and is equal to half the bit rate (or symbol rate, as both are the same). The frequency deviation is then given as Δf=Rb/4. The two tones have frequencies of fc±Δf where fc is the carrier frequency. MSK is sometimes also visualized as Offset QPSK (OQPSK) but we will not go into its details here.

When a wireless signal travels from a transmitter (Tx) to a receiver (Rx) it undergoes some changes. In simple terms the signal s(t) is scaled by a factor h(t) and noise n(t) is added at the receiver. Let’s take this discussion forward with a simple example. Suppose the Tx transmits one of two possible symbols, +1 or -1. In technical lingo this is called Binary Phase Shift Keying (BPSK). If the channel scaling factor is 0.1 we will either get a +0.1 or -0.1 at the Rx to which AWGN noise is added. The noise is random in nature (having a Gaussian distribution) but for simplicity we assume that it can have one of two values, +0.01 or -0.01.

When wireless signals travel from a transmitter to a receiver they do so after reflection, refraction, diffraction and scattering from the environment. Very rarely is there a direct line of sight (LOS) between the transmitter and receiver. Thus multiple time delayed copies of the signal reach the receiver that combine constructively and destructively. In a sense the channel acts as an FIR (finite impulse response) filter. Furthermore since the transmitter or receiver may be in motion the amplitude and phase of these replicas varies with time.

There are several methods to model the amplitude and phase of each of these components. We look at one method called the “Smiths Fading Simulator” which is based on Clark and Gans model. The simulator can be constructed using the following steps.

1. Define N the number of Gaussian RVs to be generated, fm the Doppler frequency in Hz, fs the sampling frequency in Hz, df the frequency spacing which is calculated as df=(2*fm)/(N-1) and M total number of samples in frequency domain which is calculated as M=(fs/df).
2. Generate two sequences of N/2 complex Gaussian random variables. These correspond to the frequency bins up to fm. Take the complex conjugate of these sequences to generate the N/2 complex Gaussian random variables for the negative frequency bins up to -fm.
3. Multiply the above complex Gaussian sequences g1 and g2 with square root of the Doppler Spectrum S generated from -fm to fm. Calculate the spectrum at -fm and +fm by using linear extrapolation.
4. Extend the above generated spectra from -fs/2 to +fs/2 by stuffing zeros from -fs/2 to -fm and fm to fs/2. Take the IFFT of the resulting spectra X and Y resulting in time domain signals x and y.
5. Add the absolute values of the resulting signals x and y in quadrature. Take the absolute value of this complex signal. This is the desired Rayleigh distributed envelope with the required temporal correlation.

The Matlab code for generating Rayleigh random sequence with a Doppler frequency of fm Hz is given below.

We have already seen in previous posts that the BER of BPSK increases significantly when the channel changes from a simple AWGN channel to a fading channel. One solution to this problem, that was proposed by Alamouti, was to use Transmit Diversity i.e. multiple transmit antennas transmit the information over multiple time slots increasing the likelihood of receiving the information. We have considered the simplest case of two transmit antennas and BPSK modulation (QPSK modulation would give the same BER with twice the throughput). Given below is the Python code for this, feel free to modify it and run it from the console given below.

We have previously calculated the bit error rate of BPSK in an AWGN channel, we now do the same for a Rayleigh fading channel. Remember that we have now shifted our focus from MATLAB to Python since its open and free to use. We are currently using Python-2 but intend to Python-3 once some integration issues with Trinket are sorted out.

We have previously looked at frequency domain fading simulators i.e. simulators that define the Doppler components in the frequency domain and then perform an IDFT to get the time domain signal. These simulators include Smith’s Simulator, Young’s Simulator and our very own Computationally Efficient Rayleigh Fading Simulator. Another technique that has been widely reported in the literature is Sum of Sinusoids Method. As the name suggests this method generates the Doppler components in the time domain and then sums them up to generate the time domain fading envelope. There are three parameters that define the properties of the generated signal.

1) Number of sinusoids – Higher the number better the properties of the generated signal but greater the computational complexity
2) Angle of arrival – This can be generated statistically or deterministically, spread from –pi to pi.
3) Phase of the arriving wave – This is uniformly distributed between –pi and pi.

The MATLAB code below gives three similar sum of sinusoids techniques for generating a Rayleigh faded envelope [1].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUM OF SINUSOIDS FADING SIMULATORS
% fd - Doppler frequency
% fs - Sampling frequency
% ts - Sampling period
% N - Number of sinusoids
%
% www.raymaps.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
fd=70;
fs=1000000;
ts=1/fs;
t=0:ts:1;
N=100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 1 - Clarke
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(1,length(t));
y=zeros(1,length(t));
for n=1:N;n
alpha=(rand-0.5)*2*pi;
phi=(rand-0.5)*2*pi;
x=x+randn*cos(2*pi*fd*t*cos(alpha)+phi);
y=y+randn*sin(2*pi*fd*t*cos(alpha)+phi);
end
z=(1/sqrt(N))*(x+1i*y);
r1=abs(z);
plot(t,10*log10(r1))
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 2 - Pop, Beaulieu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(1,length(t));
y=zeros(1,length(t));
for n=1:N;n
alpha=2*pi*n/N;
phi=(rand-0.5)*2*pi;
x=x+randn*cos(2*pi*fd*t*cos(alpha)+phi);
y=y+randn*sin(2*pi*fd*t*cos(alpha)+phi);
end
z=(1/sqrt(N))*(x+1i*y);
r2=abs(z);
plot(t,10*log10(r2),'r')
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 3 - Chengshan Xiao
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(1,length(t));
y=zeros(1,length(t));
for n=1:N;n
phi=(rand-0.5)*2*pi;
theta=(rand-0.5)*2*pi;
alpha=(2*pi*n+theta)/N;
x=x+randn*cos(2*pi*fd*t*cos(alpha)+phi);
y=y+randn*sin(2*pi*fd*t*cos(alpha)+phi);
end
z=(1/sqrt(N))*(x+1i*y);
r3=abs(z);
plot(t,10*log10(r3),'g')
hold off
xlabel('Time(sec)')
ylabel('Envelope(dB)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

All the three techniques given above are quite accurate in generating a Rayleigh faded envelope with the desired statistical properties. The accuracy of these techniques increases as the number of sinusoids goes to infinity (we have tested these techniques with up to 1000 sinusoids but realistically speaking even 100 sinusoids are enough). If we want to compare the three techniques in terms of the Level Crossing Rate (LCR) and Average Fade Duration (AFD) we can say that the first and third technique are a bit more accurate than the second technique. Therefore we can conclude that a statistically distributed angle of arrival is a better choice than a deterministically distributed angle of arrival. Also, if we look at the autocorrelation of the in-phase and quadrature components we see that for the first and third case we get a zero order Bessel function of the first kind whereas for the second case we get a somewhat different sequence which approximates the Bessel function with increasing accuracy as the number of sinusoids is increased.

The above figures show the theoretical Bessel function versus the autocorrelation of the real/imaginary part generated by method number two. The figure on the left considers 20 sinusoids whereas the figure on the right considers 40 sinusoids. As can be seen the accuracy of the autocorrelation sequence increases considerably by doubling the number of sinusoids. We can assume that for number of sinusoids exceeding 100 i.e. N=100 in the above code the generated autocorrelation sequence would be quite accurate.

[1] Chengshan Xiao, “Novel Sum-of-Sinusoids Simulation Models for Rayleigh and Rician Fading Channels,” IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 5, NO. 12, DECEMBER 2006.

In the previous posts we had discussed generation of a correlated Rayleigh fading sequence using Smith’s method [1] and Young’s modification of Smith’s method [2]. The main contribution of Young was that he proposed a mechanism where the number of IDFTs was reduced by half. This was achieved by first adding two length N IID zero mean Gaussian sequences filtered by the filter F[k] and then performing the IDFT on the resulting complex sequence.

This was different to Smith’s method where the IDFT was performed simultaneously on two branches and then the outputs of these branches were added in quadrature to achieve the desired sequence with Rayleigh distributed envelope and Uniformly distributed phase. Another problem with Smith’s method was that the outputs of the two arms after performing IDFTs was assumed to be real which is not always the case in implementation and depends upon the combination of Doppler frequency (fm) and length of Gaussian sequence (N).

Young’s technique is shown graphically in the above figure. Also shown is the definition of filter F[k] which depends upon N, fm and km (please note that the fm in the above equation is normalized by the sampling frequency fs). Here km = N*(fm/fs). We propose three modifications to Young’s technique which significantly reduces computation and at the same time maintains the statistical properties of the generated sequence. The modifications we propose are.

1. First modification has to do with the generated Gaussian sequence. It is observed that the filter F[k], at very high sampling rates, is mostly zero and there are very few points which have some non-zero value. So when we multiply the Gaussian sequence with the filter we mostly get zeros at the output. So we propose that the filter response in the frequency domain must be calculated first and the the Gaussian random sequence must be generated for only those points where the filter F[k] is non-zero e.g. for a sampling frequency of 7.68 MHz (standard sampling frequency for a BW of 5 MHz in LTE) and Doppler frequency of 70 Hz (corresponding to medium Doppler case in LTE) the filter F[k] has 99.9982% zeros in its frequency response and it would be a highly wasteful to calculate Gaussian RVs for all those values.

2. Secondly according to Clarke [3] the Doppler Spectrum measurements have “Marked disagreement at very low frequencies and at frequencies in the region of the sharp cut-off associated with the maximum Doppler frequency shift. At very low frequencies the spectral energy is always observed to be higher than that predicted by theory”. He goes on to add “The reason for this is that neither theoretical model takes into account the large scale variations in total energy which result from the changing topography between transmitter and mobile receiver”. This suggests that the Classical Doppler Spectrum might not be the best choice under all scenarios. This has also been noted in [4] where a flat fading model is evaluated in terms of its Level Crossing Rate and Average Fade Duration. Such a flat spectrum is especially suited to indoor scenarios as noted in [5] and [6].

We propose a filter that gives equal weight to all the frequencies up to the maximum Doppler frequency. So our filter is a box-type filter which applies a constant scaling factor to all the frequencies in the pass-band and zeros out all the frequencies in the stop-band. So in fact the Gaussian sequence that is generated in the in-phase arm may directly be added with the Gaussian sequence from the other arm without applying the frequency domain filter and then IDFT of the complex sequence is taken. We will look into the deviation from ideality that this causes later.

3. The third modification that we propose is in the implementation of the IDFT. Here again we take into consideration that the complex sequence being fed to the IDFT is filled with zeros (as we noted earlier 99.9982% zeros for 7.68 MHz and even more for higher frequencies) so we can avoid a lot multiplications and summations. The IDFT is defined below and also given is our modification to it.

Further improvement in computation time is achieved by implementing the above as a matrix multiplication. The matrix multiplication is implemented as H*X where H is the IDFT coefficients of size N x 2(km+1) and X is a vector of size 2(km+1) x 1 upon which the IDFT has to be performed.

Now let us look at the output sequence generated by using the above techniques. We consider the case of Medium Doppler Frequency of 70 Hz (EVA channel) as defined by LTE specifications. Sampling frequency is fixed at 10 kHz giving a normalized Doppler frequency of 0.007. This was done due to limitation of memory on the machine. The author also experimented with a sampling frequency of 7.68 MHz but this did not yield enough samples for statistically accurate results. We did use a sampling frequency of 7.68 MHz for our bit error rate simulation which is shown in the end.

It is observed for fm=70 Hz the envelope of the output sequence using the proposed technique matches quite well with the envelope of the output sequence generated by the ideal filter proposed by Young. Also the phase and envelope of the sequence generated using the proposed technique has the desired distribution. Some of the other metrics that we can look at are the level crossing rate (LCR) and average fade duration (AFD) as well as the Auto Correlation of the real and imaginary parts of the complex sequence generated.

Parameter

Young

Modified

LCR (ideal)

48.1086

48.1086

LCR (sim)

48.1506

39.4348

AFD (ideal)

0.0018

0.0018

AFD (sim)

0.0018

0.0022

If we look at the results for LCR and AFD we see that the simulated results match reasonably well with the results predicted by theory. These results correspond to 100 snapshots of the fading sequence. It was important to take the average of several snapshots as results varied with each simulation run. Sometimes Young’s technique produced more accurate results while at other times the proposed technique was better. Again the limitations of computer memory and processing power dictated the length of the sequence that could be generated.

In general Young’s technique produced better results than our proposed technique. It was found that product of LCR and AFD for both cases matched quite well with the theoretical value. So the total time spent in a fade state per second was equal in both the cases. In the proposed method the duration of a single fading event was higher, whereas the number of fading events per second was lower. This can be attributed to the fact that in our proposed technique higher weighting is given to lower frequency components and the fading sequence is smoothed out by these low frequency components. One technique to overcome this is spectral broadening as suggested by [4] but this is not the subject here and we postpone its discussion to another article.

Auto Correlation of Real Part fm=70Hz

The Auto Correlation of the real and imaginary parts of the generated sequences are also calculated for a Doppler frequency of 70 Hz. It is found that the Auto Correlation sequence for the two techniques matches quite well. However, the Auto Correlation sequence deviates from the theoretical value as calculated the by Bessel function of the first kind and zero order. Since we have used a flat spectral mask the Auto Correlation function resembles the sinc(x) function which is the same as zeroth order Spherical Bessel Function of the first kind (which is related to 1/2 order Bessel Function of the first kind).

It was found that when the Rayleigh fading sequence is generated by the program provided in Young’s thesis the shape of the Auto Correlation function depends upon the sampling frequency. At a normalized Doppler frequency of 0.05 and N=2^16 Young’s technique produces quite accurate results. We also measured the mean squared error (MSE) between the two Auto Correlations sequences and found it to be a function of the Normalized Doppler Frequency. It was found that as the Normalized Doppler Frequency was increased from 0.00007 to 0.0007 the MSE error dropped from 0.0277 to 0.0041. The relationship between the Normalized Doppler Frequency and MSE, for a fixed sequence length, seems to be resembling an exponential function. For more accurate results, at higher sampling frequencies, the number of samples would have to be increased considerably. In fact it was found that if the variable km (km=N*fm/fs) is maintained at around 20 the error between the two correlation sequences is less than 1% for all possible sampling rates.

We also compared the bit error rate (BER) performance of different QAM modulation schemes using both the techniques for fading envelope generation and found these to be matching quite well. A single tap was used which results in a flat fading channel. This is a simplistic channel model but it gives us some confidence that the proposed approach does have the desired statistical properties. A good test of a temporally correlated Rayleigh fading sequence is to test it on a system that implements interleavers and channel coders whose performance strongly depends upon factors such as the LCR and AFD e.g. a certain forward error correction (FEC) code might work well in high LCR and low AFD as this distributes out errors in different blocks and allows the code to correct them. In simulations done so far (not shown here) we have found that for a 1/2 rate convolutional encoder with Hard Viterbi Decoding the BER for the two schemes matches quite well. In general the results for correlated fading are much worse than uncorrelated fading.

In future we would also like to evaluate the bit error rate (BER) performance of an M-QAM OFDM system with Frequency Selective Rayleigh fading as described by the LTE fading channels EPA, EVA and ETU. This is probably a good scenario to compare the accuracy of the two techniques used to generate Rayleigh fading sequences above. One challenge in this regard is that the LTE channel taps are described in increments of 10 nsec whereas the LTE signal sampling rate is defined on a different scale (minimum Ts=32.5521 nsec corresponding to a sampling rate of 30.72 Msps). So we would have to do sample rate conversion to implement a time varying frequency selective Rayleigh fading channel.

[1] John I. Smith, “A Computer Generated Multipath Fading Simulation for Mobile Radio”, IEEE Transactions on Vehicular Technology, vol VT-24, No. 3, August 1975.

[2] David J. Young and Norman C. Beaulieu, “The Generation of Correlated Rayleigh Random Variates by Inverse Discrete Fourier Transform”, IEEE Transactions on Communications vol. 48 no. 7 July 2000.

[3] R. H. Clarke, A Statistical Theory of Mobile Radio Reception”, Bell Systems Technical Journal 47 (6), pp 957–1000, July 1968.

[4] Rosmansyah, Y.; Saunders, S.R.; Sweeney, P.; Tafazolli, R., “Equivalence of flat and classical Doppler sample generators,” Electronics Letters , vol.37, no.4, pp.243,244, 15 Feb 2001.

[5] JTC (Joint Technical Committee T1 RIP1.4 and TIA TR46.3.3/TR45.4.4 on Wireless Access): “Draft final report on RF channel characterization”. Paper no. JTC(AIR)/94.01.17-238R4, January 17, 1994.

[6] ETSI (European Telecommunications Standards Institute): “Universal mobile telecommunications system (UMTS); selection procedures for the choice of radio transmission technologies of the UMTS (UMTS 30.03 version 3.2.0)”. Technical Report, TR 101 112 V3.2.0 (1998-04), http://www.etsi.org, 1998.

It is sometimes important to know the relationship between various distributions. This can be useful if there is a function available for one distribution and it can be used to derive other distributions. In the context of Wireless Communications it is important to know the relationship between the Uniform, Gaussian and Rayleigh distribution.

According to Central Limit Theorem the sum of a large number of independent and identically distributed random variables has a Gaussian distribution. This is used to model the amplitude of the in-phase and quadrature components of a wireless signal. Shown below is the model for the received signal which has been modulated by the Gaussian channel coefficients g1 and g2.

r=g1*a1*cos(2*pi*fc*t)+g2*a2*sin(2*pi*fc*t)

The envelope of this signal (sqrt(g1^2+g2^2)) as a Rayleigh distribution. Now if you only had a function for Uniform Distribution you can generate Rayleigh Distribution using the following routine.

clear all
close all
M=10000;
N=100;
for n=1:M;
x1=rand(1,N)-0.5;
x2=rand(1,N)-0.5;
y1=mean(x1);
y2=mean(x2);
z(n)=sqrt(y1^2+y2^2);
end
hist(z,20)

Note: Here a1 and a2 can be considered constants (at least during the symbol duration) and its really g1 and g2 that are varying.