This code is provided as supplementary material of the lecture Channel Coding 2 - Advanced Methods.
This code illustrates
In [2]:
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
Conditional pdf $f_{Y|X}(y|x)$ for a channel with noise variance (per dimension) $\sigma_n^2$. This is merely the Gaussian pdf with mean $x$ and variance $\sigma_n^2$
In [3]:
def f_YgivenX(y,x,sigman):
return np.exp(-((y-x)**2)/(2*sigman**2))/np.sqrt(2*np.pi)/sigman
Output pdf $f_Y(y) = \frac12[f_{Y|X}(y|X=+1)+f_{Y|X}(y|X=-1)]$
In [4]:
def f_Y(y,sigman):
return 0.5*(f_YgivenX(y,+1,sigman)+f_YgivenX(y,-1,sigman))
This is the function we like to integrate, $f_Y(y)\cdot\log_2(f_Y(y))$. We need to take special care of the case when the input is 0, as we defined $0\cdot\log_2(0)=0$, which is usually treated as "nan"
In [5]:
def integrand(y, sigman):
value = f_Y(y,sigman)
if value < 1e-20:
return_value = 0
else:
return_value = value * np.log2(value)
return return_value
Compute the capacity using numerical integration. We have \begin{equation*} C_{\text{BI-AWGN}} = -\int_{-\infty}^\infty f_Y(y)\log_2(f_Y(y))\mathrm{d}y - \frac12\log_2(2\pi e\sigma_n^2) \end{equation*}
In [6]:
def C_BIAWGN(sigman):
# numerical integration of the h(Y) part
integral = integrate.quad(integrand, -np.inf, np.inf, args=(sigman))[0]
# take into account h(Y|X)
return -integral - 0.5*np.log2(2*np.pi*np.exp(1)*sigman**2)
This is an alternative way of calculating the capacity by approximating the integral using the Gauss-Hermite Quadrature (https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature). The Gauss-Hermite quadrature states that \begin{equation*} \int_{-\infty}^\infty e^{-x^2}f(x)\mathrm{d}x \approx \sum_{i=1}^nw_if(x_i) \end{equation*} where $w_i$ and $x_i$ are the respective weights and roots that are given by the Hermite polynomials.
We have to rearrange the integral $I = \int_{-\infty}^\infty f_Y(y)\log_2(f_Y(y))\mathrm{d}y$ a little bit to put it into a form suitable for the Gauss-Hermite quadrature \begin{align*} I &= \frac{1}{2}\sum_{x\in\{\pm 1\}}\int_{-\infty}^\infty f_{Y|X}(y|X=x)\log_2(f_Y(y))\mathrm{d}y \\ &= \frac{1}{2}\sum_{x\in\{\pm 1\}}\int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}\sigma_n}e^{-\frac{(y-x)^2}{2\sigma_n^2}}\log_2(f_Y(y))\mathrm{d}y \\ &\stackrel{(a)}{=} \frac{1}{2}\sum_{x\in\{\pm 1\}}\int_{-\infty}^\infty \frac{1}{\sqrt{\pi}}e^{-z^2}\log_2(f_Y(\sqrt{2}\sigma_n z + x))\mathrm{d}z \\ &\approx \frac{1}{2\sqrt{\pi}}\sum_{x\in\{\pm 1\}} \sum_{i=1}^nw_i \log_2(f_Y(\sqrt{2}\sigma_n x_i + x)) \end{align*} where in $(a)$, we substitute $z = \frac{y-x}{\sqrt{2}\sigma}$
In [29]:
# alternative method using Gauss-Hermite Quadrature (see https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature)
# use 40 components to approximate the integral, should be sufficiently exact
x_GH, w_GH = np.polynomial.hermite.hermgauss(40)
def C_BIAWGN_GH(sigman):
integral_xplus1 = np.sum(w_GH * [np.log2(f_Y(np.sqrt(2)*sigman*xi + 1, sigman)) for xi in x_GH])
integral_xminus1 = np.sum(w_GH * [np.log2(f_Y(np.sqrt(2)*sigman*xi - 1, sigman)) for xi in x_GH])
integral = (integral_xplus1 + integral_xminus1)/2/np.sqrt(np.pi)
return -integral - 0.5*np.log2(2*np.pi*np.exp(1)*sigman**2)
Compute the capacity for a range of of $E_s/N_0$ values (given in dB)
In [33]:
esno_dB_range = np.linspace(-16,10,100)
# convert dB to linear
esno_lin_range = [10**(esno_db/10) for esno_db in esno_dB_range]
# compute sigma_n
sigman_range = [np.sqrt(1/2/esno_lin) for esno_lin in esno_lin_range]
capacity_BIAWGN = [C_BIAWGN(sigman) for sigman in sigman_range]
# capacity of the AWGN channel
capacity_AWGN = [0.5*np.log2(1+1/(sigman**2)) for sigman in sigman_range]
Plot the capacity curves as a function of $E_s/N_0$ (in dB) and $E_b/N_0$ (in dB). In order to calculate $E_b/N_0$, we recall from the lecture that \begin{equation*} \frac{E_s}{N_0} = r\cdot \frac{E_b}{N_0}\qquad\Rightarrow\qquad\frac{E_b}{N_0} = \frac{1}{r}\cdot \frac{E_s}{N_0} \end{equation*} Next, we know that the best rate that can be achieved is the capacity, i.e., $r=C$. Hence, we get $\frac{E_b}{N_0}=\frac{1}{C}\cdot\frac{E_s}{N_0}$. Converting to decibels yields \begin{align*} \frac{E_b}{N_0}\bigg|_{\textrm{dB}} &= 10\cdot\log_{10}\left(\frac{1}{C}\cdot\frac{E_s}{N_0}\right) \\ &= 10\cdot\log_{10}\left(\frac{1}{C}\right) + 10\cdot\log_{10}\left(\frac{E_s}{N_0}\right) \\ &= \frac{E_s}{N_0}\bigg|_{\textrm{dB}} - 10\cdot\log_{10}(C) \end{align*}
In [34]:
fig = plt.figure(1,figsize=(15,7))
plt.subplot(121)
plt.plot(esno_dB_range, capacity_AWGN)
plt.plot(esno_dB_range, capacity_BIAWGN)
plt.xlim((-10,10))
plt.ylim((0,2))
plt.xlabel('$E_s/N_0$ (dB)',fontsize=16)
plt.ylabel('Capacity (bit/channel use)',fontsize=16)
plt.grid(True)
plt.legend(['AWGN','BI-AWGN'],fontsize=14)
# plot Eb/N0 . Note that in this case, the rate that is used for calculating Eb/N0 is the capcity
# Eb/N0 = 1/r (Es/N0)
plt.subplot(122)
plt.plot(esno_dB_range - 10*np.log10(capacity_AWGN), capacity_AWGN)
plt.plot(esno_dB_range - 10*np.log10(capacity_BIAWGN), capacity_BIAWGN)
plt.xlim((-2,10))
plt.ylim((0,2))
plt.xlabel('$E_b/N_0$ (dB)',fontsize=16)
plt.ylabel('Capacity (bit/channel use)',fontsize=16)
plt.grid(True)
In [35]:
from scipy.stats import norm
# first compute the BSC error probability
# the Q function (1-CDF) is also often called survival function (sf)
delta_range = [norm.sf(1/sigman) for sigman in sigman_range]
capacity_BIAWGN_hard = [1+delta*np.log2(delta)+(1-delta)*np.log2(1-delta) for delta in delta_range]
fig = plt.figure(1,figsize=(15,7))
plt.subplot(121)
plt.plot(esno_dB_range, capacity_AWGN)
plt.plot(esno_dB_range, capacity_BIAWGN)
plt.plot(esno_dB_range, capacity_BIAWGN_hard)
plt.xlim((-10,10))
plt.ylim((0,2))
plt.xlabel('$E_s/N_0$ (dB)',fontsize=16)
plt.ylabel('Capacity (bit/channel use)',fontsize=16)
plt.grid(True)
plt.legend(['AWGN','BI-AWGN', 'Hard BI-AWGN'],fontsize=14)
# plot Eb/N0 . Note that in this case, the rate that is used for calculating Eb/N0 is the capcity
# Eb/N0 = 1/r (Es/N0)
plt.subplot(122)
plt.plot(esno_dB_range - 10*np.log10(capacity_AWGN), capacity_AWGN)
plt.plot(esno_dB_range - 10*np.log10(capacity_BIAWGN), capacity_BIAWGN)
plt.plot(esno_dB_range - 10*np.log10(capacity_BIAWGN_hard), capacity_BIAWGN_hard)
plt.xlim((-2,10))
plt.ylim((0,2))
plt.xlabel('$E_b/N_0$ (dB)',fontsize=16)
plt.ylabel('Capacity (bit/channel use)',fontsize=16)
plt.grid(True)