In [6]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2, ncx2
from scipy.special import erfc
from scipy.integrate import quad
np.seterr(divide='ignore');

In [2]:
N_fft = 256
fs = 2000
Rb = 250/32

In [25]:
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)[0]
    return 1-p-(1-p)*(N_fft/2)/(N_fft-1)

In [13]:
quad(myfun, args = (5,), a = -np.inf, b = np.inf)


Out[13]:
(0.4401048633943861, 6.037367817273872e-09)

In [15]:
lrtc_ber_formula2(5)


Out[15]:
0.38894042946897556

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

In [26]:
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]);