Category Archives: Simulation

BER for BPSK-OFDM in Frequency Selective Channel

OFDM Tx-Rx Block Diagram

As the data rates supported by wireless networks continue to rise the bandwidth requirements also continue to increase (although spectral efficiency has also improved). Remember GSM technology which supported 125 channels of 200KHz each, which was further divided among eight users using TDMA. Move on to LTE where the channel bandwidth could be as high as 20MHz (1.4MHz, 3MHz, 5MHz, 10MHz, 15MHz and 20MHz are standardized).

This advancement poses a unique challenge referred to as frequency selective fading. This means that different parts of the signal spectrum would see a different channel (different amplitude and different phase offset). Look at this in the time domain where the larger bandwidth means shorter symbol period causing intersymbol interference (as time delayed copies of the signal overlap on arrival at the receiver).

The solution to this problem is OFDM that divides the wideband signal into smaller components each having a bandwidth of a few KHz. Each of these components experiences a flat channel. To make the task of equalization simple a cyclic prefix (CP) is added in the time domain to make the effect of fading channel appear as circular convolution. Thus simplifying the frequency domain equalization to a simple division operation.

Shown below is the Python code that calculates the bit error rate (BER) of BPSK-OFDM which is the same as simple BPSK in a Rayleigh flat fading channel. However there is a caveat. We have inserted a CP which means we are transmitting more energy than simple BPSK. To be exact we are transmitting 1.25 (160/128) times more energy. This means that if this excess energy is accounted for the performance of BPSK-OFDM would be 1dB (10*log10(1.25)) worse than simple BPSK in Rayleigh flat fading channel.

Note:

  1. Although we have shown the channel as a multiplicative effect in the figure above, this is only true for a single tap channel. For a multi-tap channel (such as the one used in the code above) the effect of the channel is that of a filter which performs convolution operation on the transmitted signal.
  2. We have used a baseband model in our simulation and the accompanying figure. In reality the transmitted signal is upconverted before transmission by the antennas.
  3.  The above model can be easily modified for any modulation scheme such as QPSK or 16-QAM. The main difference would be that the signal would have a both a real part and an imaginary part, much of the simulation would remain the same. This would be the subject of a future post. For a MATLAB implementation of 64-QAM OFDM see the following post (64-QAM OFDM).
  4. Serial to parallel and parallel to serial conversion shown in the above figure was not required as the simulation was done symbol by symbol (one OFDM symbol in the time domain represented 128 BPSK symbols in the frequency domain).
  5. The channel model in the above simulation is quasi-static i.e. it remains constant for one OFDM symbol but then rapidly changes for the next, without any memory.

Rayleigh Fading Envelope Generation – Python

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.

Alamouti – Transmit Diversity Scheme – Implemented in Python

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.

Implementation on Trinket

Implementation on REPL

Run Python Code from the Browser

Here is a piece of Python code that calculates Bit Error Rate (BER) of BPSK. The code is a bit slow at the moment, compared to MATLAB implementation, but this is work in progress and further optimizations would be carried out. We would like to point out that the main reason for this slower implementation is that a bit by bit error calculation is done, instead of a vectorial implementation. We already pointed out in our previous post that a “for loop” implemented in Python is not that efficient.

BPSK Bit Error Rate Calculation Using Python

Have you ever thought about how life would be without MATLAB. As it turns out there are free and open source options such as Python. We have so far restricted ourself to MATLAB in this blog but now we venture out to find out what are the other options. Given below is a most basic Pyhton code that calculates the Bit Error Rate of Binary Phase Shift Keying (BPSK). Compare this to our MATLAB implementation earlier [BPSK BER].

There are various IDEs available for writing your code but I have used Enthought Canopy Editor (32 bit) which is free to download and is also quite easy to use [download here]. So as it turns out that there is life beyond MATLAB. In fact there are several advantages of using Python over MATLAB which we will discuss later in another post. Lastly please note the indentation in the code below as there is no “end” statement in a for loop in Python.

from numpy import sqrt
from numpy.random import rand, randn
import matplotlib.pyplot as plt
  
N = 5000000
EbNodB_range = range(0,11)
itr = len(EbNodB_range)
ber = [None]*itr

for n in range (0, itr): 
 
    EbNodB = EbNodB_range[n]   
    EbNo=10.0**(EbNodB/10.0)
    x = 2 * (rand(N) >= 0.5) - 1
    noise_std = 1/sqrt(2*EbNo)
    y = x + noise_std * randn(N)
    y_d = 2 * (y >= 0) - 1
    errors = (x != y_d).sum()
    ber[n] = 1.0 * errors / N
    
    print "EbNodB:", EbNodB
    print "Error bits:", errors
    print "Error probability:", ber[n] 
        
plt.plot(EbNodB_range, ber, 'bo', EbNodB_range, ber, 'k')
plt.axis([0, 10, 1e-6, 0.1])
plt.xscale('linear')
plt.yscale('log')
plt.xlabel('EbNo(dB)')
plt.ylabel('BER')
plt.grid(True)
plt.title('BPSK Modulation')
plt.show()
BPSK Bit Error Rate
BPSK Bit Error Rate

MATLAB vs PYTHON A COMPARISON

 

Knife Edge Diffraction Model

Diffraction is a phenomenon where electromagnetic waves (such as light waves) bend around corners to reach places which are otherwise not reachable i.e. not in the line of sight. In technical jargon such regions are also called shadowed regions (the term again drawn from the physics of light). This phenomenon can be explained by Huygen’s principle which states that as a plane wave propagates in a particular direction each new point along the wavefront is a source of secondary waves. This can be understood by looking at the following figure. However one peculiarity of this principle is that it is unable to explain why the new point source transmits only in the forward direction.

Image result for diffraction

The electromagnetic field in the shadowed region can be calculated by combining vectorially the contributions of all of these secondary sources, which is not an easy task. Secondly, the geometry is usually much more complicated than shown in the above figure. For example consider a telecom tower transmitting electromagnetic waves from a rooftop and a pedestrian using a mobile phone at street level. The EM waves usually reach the receiver at street level after more than one diffraction (not to mention multiple reflections).  However, an approximation that works well in most cases is called knife edge diffraction, which assumes a single sharp edge (an edge with a thickness much smaller than the wavelength) separates the transmitter and receiver.

The path loss due to diffraction in the knife edge model is controlled by the Fresnel Diffraction Parameter which measures how deep the receiver is within the shadowed region. A negative value for the parameter shows that the obstruction is below the line of sight and if the value is below -1 there is hardly any loss. A value of 0 (zero) means that the transmitter, receiver and tip of the obstruction are all in line and the Electric Field Strength is reduced by half or the power is reduced to one fourth of the value without the obstruction i.e. a loss of 6dB.  As the value of the Fresnel Diffraction Parameter increases on the positive side the path loss rapidly increases reaching a value of 27 dB for a parameter value of 5. Sometimes the exact calculation is not needed and only an approximate calculation, as proposed by Lee in 1985, is sufficient.

Fresnel Diffraction Parameter (v) is defined as:

v=h*sqrt(2*(d1+d2)/(lambda*d1*d2))

where

d1 is the distance between the transmitter and the obstruction along the line of sight

d2 is the distance between the receiver and the obstruction along the line of sight

h is the height of the obstruction above the line of sight

and lambda is the wavelength

Diffraction Loss Using Knife-Edge Model

The MATLAB codes used to generate the above plots are given below (approximate method followed by the exact method). Feel free to use them in your simulations and if you have a question drop us a comment.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of the path loss based on the value of
% Fresnel Diffraction Parameter as proposed by Lee
% Lee W C Y Mobile Communications Engineering 1985
% Copyright www.raymaps.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all

v=-5:0.01:5;

for n=1:length(v)

    if v(n) <= -1
        G(n)=0;
    elseif v(n) <= 0
        G(n)=20*log10(0.5-0.62*v(n));
    elseif v(n) <= 1
        G(n)=20*log10(0.5*exp(-0.95*v(n)));
    elseif v(n) <= 2.4
        G(n)=20*log10(0.4-sqrt(0.1184-(0.38-0.1*v(n))^2));
    else
        G(n)=20*log10(0.225/v(n));
    end
    
end

plot(v, G, 'b')
xlabel('Fresnel Diffraction Parameter')
ylabel('Diffraction Loss (dB)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exact calculation of the path loss (in dB)
% based on Fresnel Diffraction Parameter (v)
% T S Rappaport Wireless Communications P&P
% Copyright www.raymaps.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all

v=-5:0.01:5;

for n=1:length(v)
    
    v_vector=v(n):0.01:v(n)+100;
    
    F(n)=((1+1i)/2)*sum(exp((-1i*pi*(v_vector).^2)/2));
    
end

F=abs(F)/(abs(F(1)));

plot(v, 20*log10(F),'r')
xlabel('Fresnel Diffraction Parameter')
ylabel('Diffraction Loss (dB)')

An interesting video explaining the phenomenon of diffraction.

Reflection vs Scattering in Ray-Tracing

In ray-tracing simulations one is faced with a complexity at points of intersection between rays and objects. How to handle the reflection, should it be specular or scattered. In this post we consider the two extreme cases; pure reflection vs scattering.

For this we consider that in the case of specular reflection the power of the ray is reduced by 3dB (that is R=0.5) and the ray continues as it would if there was no interaction with the object (off course direction would be changed with angle of reflection being equal to angle of incidence).

In the case of scattered ray we assume that the point of interaction between the ray and the object is a point source from where the rays are regenerated. Therefore there is a rapid drop in the E-field strength as the ray propagates away from the point of interaction.

clear all
close all

Eo=1;
r1=1:100;
r2=101:200;
r=[r1 r2];
R=0.5;

Er1=Eo./r1;
Er2=R*Eo./r2;
Er=[Er1 Er2];

Es1=Eo./r1;
Es2=Es1(end)./(r2-r1(end));
Es=[Es1 Es2];

plot(r,10*log10(Er),'b')
hold on
plot(r,10*log10(Es),'r:')
hold off

legend('Reflection','Scattering')
xlabel('Distance (m)')
ylabel('E-field (V/m)')

Reflected vs Scattered Ray

The simulation code above considers that a ray travels unobstructed for 100 m and at this point comes in contact with an object and is reflected (reflection coefficient of 0.5 is assumed) or scattered. The ray then again travels for another 100 m without coming in contact with an object.

It is observed that there is a rapid decrease in E-field strength of the scattered ray beyond 100 m. This is because the E-field strength at the point of interaction is much lower than that at the source and when this is considered to be a point source, re-generating the rays, the resulting E-field decays quite rapidly.

We assume that the reflection is specular in most of our simulations as this is closer to reality than assuming a fully scattered ray.

Eclipse

Eclipse 1.0 – A Paradigm Shift in RF Planning

NEW: Simulation of a Moving Transmitter (such as a car)

NEW: Simulation of a Moving Transmitter (such as a pedestrian)

Radio frequency planning is an essential component of network planning, roll-out, up-gradation, expansion etc. Several methods can be adopted for this from something as simple as free space models, empirical path loss models to the significantly more complicated, time consuming and expensive drive testing. Drive testing gives very accurate results but these results can be rendered useless by changing the position of an antenna or the tilt or transmit power of an antenna requiring another run in the field. One solution to this problem is ray-tracing which is very accurate but is usually considered to be very computationally expensive and of little practical value. But recent advances in computational power of machines coupled with efficient techniques have given a new lease of life to this method.

Eclipse is a near real-time simulation software for prediction of signal strength in urban areas. The software uses shooting and bouncing ray (SBR) method of ray tracing with 1 degree ray separation, 1 m step size and 9 interactions per ray path. The simulation parameters can be varied according to the resolution required. The code is highly optimized to give results in shortest possible time. It is especially useful for network planning of ultra-dense wireless networks where a dense network of antennas is placed on lamp posts instead of telecom towers. Various frequency bands can be simulated, along with different antenna radiation patterns and MIMO configurations.


Helsinki 3D Building Data

 

Path Followed by a Single Ray

 

Paths Followed by Multiple Rays

 

Received Signal Strength Over Area of Interest

Note: If you would like to run a test simulation send us a request at info@raymaps.com

Sum of Sinusoids Fading Simulator

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.

Correlation of Real and Imaginary Parts
Correlation of Real and Imaginary Parts

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.