In [1]:
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib
# showing figures inline
%matplotlib inline
In [2]:
# plotting options
font = {'size' : 20}
plt.rc('font', **font)
plt.rc('text', usetex=True)
matplotlib.rc('figure', figsize=(18, 6) )
In [3]:
# max. numbers of errors and/or symbols
# used as loop condition for the simulation
max_errors = 1e2
max_syms = 1e4
# Eb/N0
EbN0_db_min = 0
EbN0_db_max = 12
EbN0_db_step = 2
# initialize Eb/N0 array
EbN0_db_range = np.arange( EbN0_db_min, EbN0_db_max, EbN0_db_step )
EbN0_range = 10**( EbN0_db_range / 10 )
# constellation points
mod_points_bpsk = [-1, 1]
In [4]:
# initialize BER array
# NOTE: dtype is required since zeros_like(arg) is of same type as arg
ber_bpsk = np.zeros_like( EbN0_db_range, dtype=float)
# theoretical ber for bpsk as on slides
# NOTE: using sigma^2 = N0 / 2 which will be justified later
# NOTE: Q(x) = 1 - Phi( x )
ber_bpsk_theo = 1 - stats.norm.cdf( np.sqrt( 2 * EbN0_range ) )
# loop for snr
for ind_snr, val_snr in enumerate( EbN0_range ):
# initialize error counter
num_errors_bpsk = 0
num_syms = 0
# get noise variance
sigma2 = 1. / ( val_snr )
# loop for errors
while ( num_errors_bpsk < max_errors and num_syms < max_syms ):
# generate data and modulate by look-up
d = np.random.randint( 0, 2)
s_bpsk = mod_points_bpsk[ d ]
# add noise
noise = np.sqrt( sigma2 / 2 ) * ( np.random.randn() + 1j * np.random.randn() )
r_bpsk = s_bpsk + noise
# demod by comparing real part with 0
d_est_bpsk = int( np.real( r_bpsk ) > 0 )
# count errors
num_errors_bpsk += int( d_est_bpsk != d )
# increase symbol counter
num_syms += 1
# ber as relative amount of errors
ber_bpsk[ ind_snr ] = num_errors_bpsk / num_syms
print('Es/N0 planned (dB) = {:2.1f}\n'.format( 10*np.log10(val_snr) ) )
In [5]:
# plot theoretical values
ax_bpsk = plt.plot( EbN0_db_range, ber_bpsk_theo, linewidth=2.0, label = "BPSK, theo.")
# get color and use same color for simulation
col_bpsk = ax_bpsk[0].get_color()
plt.plot( EbN0_db_range, ber_bpsk, '^', color=col_bpsk, linewidth=2.0, markersize = 12, label = "BPSK, sim." )
plt.yscale('log')
plt.grid(True)
plt.legend(loc='upper right')
plt.xlabel('$E_b/N_0$ (dB)')
plt.ylabel('BER')
Out[5]:
In [ ]: