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()

We have previously discussed the Bit Error Rate of M-QAM in Rayleigh Fading using Monte Carlo Simulation. We now turn our attention to calculation of Bit Error Rate (BER) of M-QAM in Rayleigh fading using analytical techniques. In particular we look at the method used in MATLAB function berfading.m. In this function the BER of 4-QAM, 16-QAM and 64-QAM is calculated from series expressions having 1, 3 and 5 terms respectively. These are given below (M is the constellation size and must be a power of 2).

Although using these expressions we get very accurate BER but it is not that simple to calculate (the expressions become even more complicated for higher constellation sizes such as 256-QAM). Therefore we try to simplify these expressions by using only the first term in each expression. To our surprise the results match quite well with the results using the exact formulae. There is very minor difference at low signal to noise ratios but that can be easily bargained for the ease of calculation.

So here is our program for calculating the BER using the approximate method.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION TO CALCULATE THE BER OF M-QAM IN RAYLEIGH FADING
% M: Input, Constellation Size
% EbNo: Input, Energy Per Bit to Noise Power Spectral Density
% ber: Output, Bit Error Rate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ber]= BER_QAM_fading (M, EbNo)
k=log2(M);
EbNoLin=10.^(EbNo/10);
gamma_c=EbNoLin*k;
if M==4
%4-QAM
ber = 1/2 * ( 1 - sqrt(gamma_c/k./(1+gamma_c/k)) );
elseif M==16
%16-QAM
ber = 3/8 * ( 1 - sqrt(2/5*gamma_c/k./(1+2/5*gamma_c/k)) );
elseif M==64
%64-QAM
ber = 7/24 * ( 1 - sqrt(1/7*gamma_c/k./(1+1/7*gamma_c/k)) );
else
%Warning
warning('M=4,16,64')
ber=zeros(1,length(EbNo));
end
semilogy(EbNo,ber,'o-')
xlabel('EbNo(dB)')
ylabel('BER')
axis([0 24 0.001 1])
grid on
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

So we see that the results match quite well with the results previously obtained through simulation. We will next tackle the problem of simplifying the expression for higher order modulations such as 256-QAM in both Rayleigh and Ricean channels.

In the previous post we had considered a static frequency-selective channel. We now consider a time-varying frequency selective channel with 7 taps. Each tap of the time domain filter has a Gaussian distributed real component with variance 1/(2*n_tap) and a Gaussian distributed imaginary component with variance (1/2*n_tap). The amplitude of each tap is thus Rayleigh distributed and the phase is Uniformly distributed. Since the power in each component is normalized by the filter length (n_tap) the BER performance would remain the same even if the filter length is changed (this has been verified experimentally).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION TO SIMULATE PERFORMANCE OF 64-OFDM IN TIME VARYING FREQUENCY SELECTIVE CHANNEL
% n_bits: Input, length of binary sequence
% n_fft: Input, length of FFT (Fast Fourier Transform)
% EbNodB: Input, energy per bit to noise power spectral density ratio
% ber: Output, bit error rate
% Copyright RAYmaps (www.raymaps.com)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[ber]= M_QAM_OFDM_fading(n_bits,n_fft,EbNodB)
Eb=7;
M=64;
k=log2(M);
n_cyc=32;
EbNo=10^(EbNodB/10);
x=transpose(round(rand(1,n_bits)));
h1=modem.qammod(M);
h1.inputtype='bit';
h1.symbolorder='gray';
y=modulate(h1,x);
n_sym=length(y)/n_fft;
n_tap=7;
for n=1:n_sym;
s_ofdm=sqrt(n_fft)*ifft(y((n-1)*n_fft+1:n*n_fft),n_fft);
s_ofdm_cyc=[s_ofdm(n_fft-n_cyc+1:n_fft); s_ofdm];
ht=(1/sqrt(2))*(1/sqrt(n_tap))*(randn(1,n_tap)+j*randn(1,n_tap));
Hf=fft(ht,n_fft);
r_ofdm_cyc=conv(s_ofdm_cyc,ht);
r_ofdm_cyc=(r_ofdm_cyc(1:n_fft+n_cyc));
wn=sqrt((n_fft+n_cyc)/n_fft)*(randn(1,n_fft+n_cyc)+j*randn(1,n_fft+n_cyc));
r_ofdm_cyc=r_ofdm_cyc+sqrt(Eb/(2*EbNo))*wn.';
r_ofdm=r_ofdm_cyc(n_cyc+1:n_fft+n_cyc);
s_est((n-1)*n_fft+1:n*n_fft)=(fft(r_ofdm,n_fft)/sqrt(n_fft))./Hf.';
end
h2=modem.qamdemod(M);
h2.outputtype='bit';
h2.symbolorder='gray';
h2.decisiontype='hard decision';
z=demodulate(h2,s_est.');
ber=(n_bits-sum(x==z))/n_bits
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

As before we have used an FFT size of 128 and cyclic prefix of 32 samples. The FFT and IIFT operations are normalized to maintain the signal to noise ratio (SNR). The extra energy transmitted in the cyclic prefix is also accounted for in the SNR calibration.

It is observed that the BER performance of 64-QAM OFDM in the time-varying frequency-selective channel is quite similar to that in the static frequency-selective channel with complex filter taps. It must be noted that with 64-QAM the goal is to achieve higher bit rate, error rates can be improved using antenna diversity and channel coding schemes.

Given below is the wrapper that should be used along with the above code. The wrapper basically calls the above routine for each value of EbNodB. The length of the binary sequence and the FFT size are other inputs to the function. The bit error rate at the specific EbNodB is the output of the function.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
k=6;
n_fft=128;
l=k*n_fft*1e3;
EbNodB=0:2:20;
for n=1:length(EbNodB);n
ber(n)=M_QAM_OFDM_fading(l,n_fft,EbNodB(n));
end;
semilogy(EbNodB,ber,'O-');
grid on
xlabel('EbNo')
ylabel('BER')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In future we would use the standard LTE channel models, namely EPA, EVA and ETU in our simulation.