``````

In :

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2, ncx2
from scipy.special import erfc
np.seterr(divide='ignore');

``````
``````

In :

N_fft = 256
fs = 2000
Rb = 250/32

``````
``````

In :

def lrtc_ber_sim(ebn0_dB, n_sims = 100000):
ebn0 = 10**(0.1*ebn0_dB)
snr = ebn0 * (Rb * N_fft/ fs) # by the way the factor in () is just 1, so snr = ebn0
noise = np.random.randn(n_sims, N_fft) + 1j * np.random.randn(n_sims, N_fft)
noise[:,0] += np.sqrt(2*snr) # the factor 2 is because since the noise is complex, it has power 2
# the following is a simple way to count bit errors
# it works because the FFT bin number 0 is the one in which the signal is
# present, while the remaining bins are interchangeable
return np.average(np.argmax(np.abs(noise), axis = 1) % 2)

def lrtc_ber_formula(ebn0_dB):
ebn0 = 10**(0.1*ebn0_dB)
cnr = ebn0 * Rb / fs * N_fft+1;
return 0.5*(1-chi2.cdf(cnr, 2)**N_fft)

def myfun(x, cnr_current):
return ncx2.pdf(x,2,cnr_current*2) * chi2.cdf(x,2)**(N_fft-1)

def lrtc_ber_formula2(ebn0_dB):
ebn0 = 10**(0.1*ebn0_dB)
cnr = ebn0 * Rb / fs * N_fft
p = quad(myfun, args = (cnr,), a = -np.inf, b = +np.inf)
return 1-p-(1-p)*(N_fft/2)/(N_fft-1)

``````
``````

In :

quad(myfun, args = (5,), a = -np.inf, b = np.inf)

``````
``````

Out:

(0.4401048633943861, 6.037367817273872e-09)

``````
``````

In :

lrtc_ber_formula2(5)

``````
``````

Out:

0.38894042946897556

``````
``````

In :

ebn0s = np.linspace(0, 18, 20)
bers_sim = np.array([lrtc_ber_sim(ebn0) for ebn0 in ebn0s])

``````
``````

In :

plt.figure(figsize = (12, 8), facecolor = 'w')
ebn0s_formula = np.linspace(0, 18, 1000)
plt.plot(ebn0s, np.log10(bers_sim))
plt.plot(ebn0s_formula, np.log10(lrtc_ber_formula(ebn0s_formula)))
plt.plot(ebn0s, [np.log10(lrtc_ber_formula2(e)) for e in ebn0s])
plt.plot(ebn0s_formula, np.log10(0.5 * erfc(np.sqrt(0.5*10**(0.1*ebn0s_formula)))))
plt.legend(['Simulation', 'Formula', 'Formula 2', 'FSK coherent'])
plt.xlabel('Eb/N0 (dB)')
plt.ylabel('log10(BER)')
plt.ylim([-10,0]);

``````
``````

``````