import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc
N=100000
EbN0dB=np.arange(-4,11,2)
L=16
Fc=800
Fs=L*Fc
BER=np.zeros(len(EbN0dB))
ak=np.random.randint(2,size=N)
(s_bb,t)=bpsk_mod(ak,L)
s=s_bb*np.cos(2*np.pi*Fc*t/Fs)
fig1, axs = plt.subplots(2, 2)
axs[0, 0].plot(t,s_bb)
axs[0, 0].set_xlabel('t(s)');axs[0, 1].set_ylabel(r'$s_{bb}(t)$-baseband')
axs[0, 1].plot(t,s)
axs[0, 1].set_xlabel('t(s)');axs[0, 1].set_ylabel('s(t)-with carrier')
axs[0, 0].set_xlim(0,10*L);axs[0, 1].set_xlim(0,10*L)
axs[1, 0].plot(np.real(s_bb),np.imag(s_bb),'o')
axs[1, 0].set_xlim(-1.5,1.5);axs[1, 0].set_ylim(-1.5,1.5)
for i,EbN0 in enumerate(EbN0dB):
r = awgn(s,EbN0,L)
r_bb = r*np.cos(2*np.pi*Fc*t/Fs)
ak_hat = bpsk_dmod(r_bb,L)
BER[i] = np.sum(ak !=ak_hat)/N
axs[1, 1].plot(t,r)
axs[1, 1].set_xlabel('t(s)');axs[1, 1].set_ylabel('r(t)')
axs[1, 1].set_xlim(0,10*L)
theoreticalBER = 0.5*erfc(np.sqrt(10**(EbN0dB/10)))
fig2, ax1 = plt.subplots(nrows=1,ncols = 1)
ax1.semilogy(EbN0dB,BER,'k*',label='Simulated')
ax1.semilogy(EbN0dB,theoreticalBER,'r-',label='Theoretical')
ax1.set_xlabel(r'$E_b/N_0$ (dB)')
ax1.set_ylabel(r'Probability of Bit Error - $P_b$')
ax1.set_title(['Probability of Bit Error for BPSK modulation'])
ax1.legend();