Capacity

Channel capacity: the reference axis every other metric is read against.

A code is not judged in the abstract but by its distance to capacity, and a modulation not by its SER alone but by the rate it can carry. This module provides the closed forms where they exist (Gaussian input over AWGN, Rayleigh ergodic capacity) and the numerical ones where they do not (constellation-constrained and BICM capacity, MIMO), so that a simulated curve can be placed against the limit it is chasing.

Every function returns bits per channel use.

awgn_capacity(snr: ndarray | float) → ndarray

Shannon capacity of the complex AWGN channel.

Signal Model

For \(y = x + b\) with \(b \sim \mathcal{CN}(0, \sigma^2)\) and an average power constraint \(\mathbb{E}\left[|x|^2\right] \leq P\), the capacity is reached by a Gaussian input and equals

\[C = \log_2\left(1 + \rho\right), \qquad \rho = P / \sigma^2\]

param snr:

Signal-to-noise ratio \(\rho\), in linear scale.

type snr:

float or np.ndarray

returns:

np.ndarray – Capacity in bits per channel use.

References

C. E. Shannon, “A mathematical theory of communication”, Bell System Technical Journal 27, 1948; T. M. Cover, J. A. Thomas, Elements of Information Theory, 2nd ed., Wiley, 2006, Chapter 9.

Examples

>>> print(np.round(awgn_capacity(np.array([1.0, 10.0, 100.0])), 4))
[1.     3.4594 6.6582]

constellation_capacity(alphabet: ndarray, snr: ndarray | float, *, px: ndarray | None = None, n_nodes: int = 40, method: str = 'gauss-hermite') → ndarray

Mutual information of a discrete constellation over AWGN.

Signal Model

A real system does not transmit a Gaussian: it transmits one of \(M\) constellation points, equiprobably by default. The achievable rate is then the mutual information

\[I(X; Y) = \log_2 M - \frac{1}{M} \sum_{m=1}^{M} \mathbb{E}_{b}\left[ \log_2 \sum_{m'=1}^{M} \exp\left(-\frac{|x_m - x_{m'} + b|^2 - |b|^2}{\sigma^2}\right) \right]\]

where the expectation is over the noise \(b \sim \mathcal{CN}(0, \sigma^2)\). This is the constrained capacity: it saturates at \(\log_2 M\) however high the SNR, and it is the ceiling a coded system using that constellation can approach – not the Shannon capacity, which assumes a Gaussian input.

Probabilistic shaping makes the input non-uniform, and then the \(1/M\) factors are no longer constants to pull out of the sum:

\[I(X; Y) = - \sum_{m} P_X(x_m) \; \mathbb{E}_{b}\left[ \log_2 \sum_{m'} P_X(x_{m'}) \exp\left(-\frac{|x_m - x_{m'} + b|^2 - |b|^2}{\sigma^2}\right) \right]\]

which is the same expression with px in both places and reduces to the one above when it is uniform. Passing px is how the rate of a shaped constellation is read exactly rather than estimated from samples with compute_mi().

The expectation is a two-dimensional integral against a Gaussian weight, evaluated by quadrature – no random draws. Which rule does it is left open by method; see _noise_quadrature() for what separates them.

param alphabet:

Constellation points, of unit average energy (as returned by get_alphabet()).

type alphabet:

np.ndarray

param snr:

Signal-to-noise ratio \(\rho = 1/\sigma^2\), linear.

type snr:

float or np.ndarray

param px:

Input distribution \(P_X\), of length \(M\). Default (None) is uniform. Note that the SNR is still read from the alphabet’s own scale, so a shaped input of lower energy is compared at lower power unless the alphabet is rescaled – see maxwell_boltzmann().

type px:

np.ndarray, optional, keyword-only

param n_nodes:

Quadrature nodes per noise dimension, so the cost grows as \(n^2\). Default 40.

type n_nodes:

int, optional, keyword-only

param method:

"gauss-hermite" (default), the rule matched to the Gaussian weight, or "simpson", the classical composite rule on a truncated grid. Same quantity, different numerics: the second is there to be able to check the first without leaving the library.

type method:

str, optional, keyword-only

returns:

np.ndarray – Mutual information in bits per channel use, in \([0, \log_2 M]\) – at most \(H(P_X)\) when px is given.

raises ValueError:

If method is not one of the two rules above, or if px does not have one non-negative probability per constellation point.

References

G. Ungerboeck, “Channel coding with multilevel/phase signals”, IEEE Trans. Inf. Theory 28(1), 1982; Cover & Thomas, 2nd ed., Chapter 9.

Examples

>>> from comnumpy.core.utils import get_alphabet
>>> qpsk = get_alphabet("PSK", 4)
>>> print(np.round(constellation_capacity(qpsk, np.array([1.0, 1e4])), 4))
[0.9719 2.    ]

The classical rule lands on the same number, which is the point of keeping it – it checks the default rule without leaving the library:

>>> qam16 = get_alphabet("QAM", 16)
>>> gauss = constellation_capacity(qam16, 10.0)
>>> simpson = constellation_capacity(qam16, 10.0, method="simpson")
>>> print(np.round([float(gauss), float(simpson)], 3))
[3.164 3.164]

A non-uniform input carries less than \(\log_2 M\) however clean the channel, because it is its entropy that saturates:

>>> law = np.array([0.4, 0.1, 0.1, 0.4])
>>> print(round(float(constellation_capacity(
...     get_alphabet("PSK", 4), 1e6, px=law)), 4))
1.7219

bicm_capacity(alphabet: ndarray, snr: ndarray | float, *, px: ndarray | None = None, n_nodes: int = 40, method: str = 'gauss-hermite') → ndarray

BICM capacity: the bound a bit-interleaved coded system chases.

Signal Model

A chain that maps bits to symbols, demaps to per-bit LLRs and then decodes does not achieve the constellation capacity: the parallel binary channels seen by the decoder ignore the dependence between the bits of a symbol. The achievable rate is the sum of the per-bit mutual informations

\[C_{\mathrm{BICM}} = \sum_{i=1}^{m} I(B_i; Y), \qquad m = \log_2 M\]

which is the right reference for the SymbolDemapper(soft=True) + FEC chain of this library. The gap to constellation_capacity() is the price of the bit-wise interface, and it depends on the labelling – Gray mapping minimizes it, which is why it is the default.

param alphabet:

Constellation points of unit average energy, in the labelling order used by the mapper (index = symbol value).

type alphabet:

np.ndarray

param snr:

Signal-to-noise ratio, linear.

type snr:

float or np.ndarray

param n_nodes:

Quadrature nodes per noise dimension. Default 40.

type n_nodes:

int, optional, keyword-only

param method:

"gauss-hermite" (default) or "simpson", as in constellation_capacity().

type method:

str, optional, keyword-only

returns:

np.ndarray – BICM capacity in bits per channel use.

raises ValueError:

If the constellation size is not a power of two, or method is not one of the two rules.

References

G. Caire, G. Taricco, E. Biglieri, “Bit-interleaved coded modulation”, IEEE Trans. Inf. Theory 44(3), 1998.

Examples

>>> from comnumpy.core.utils import get_alphabet
>>> qam16 = get_alphabet("QAM", 16)
>>> snr = np.array([1.0, 10.0])
>>> print(np.round(bicm_capacity(qam16, snr), 4))       # bit-wise interface
[0.8993 3.1636]
>>> print(np.round(constellation_capacity(qam16, snr), 4))  # the ceiling above it
[0.9897 3.1639]

A shaped input is passed the same way as to constellation_capacity(), and the ceiling moves with it: the sum of the per-bit entropies, not \(\log_2 M\).

>>> from comnumpy.core.shaping import maxwell_boltzmann
>>> law = maxwell_boltzmann(qam16, lam=0.5)
>>> print(np.round(bicm_capacity(qam16, snr, px=law), 4))
[0.7928 3.0653]

rayleigh_ergodic_capacity(snr: ndarray | float) → ndarray

Ergodic capacity of a flat Rayleigh fading channel.

Signal Model

With \(y = h x + b\), \(h \sim \mathcal{CN}(0,1)\) known at the receiver only, the ergodic capacity averages the instantaneous one over the fading:

\[C = \mathbb{E}_{h}\left[\log_2\left(1 + \rho |h|^2\right)\right] = \frac{1}{\ln 2} \, e^{1/\rho} \, E_1\!\left(1/\rho\right)\]

since \(|h|^2 \sim \mathrm{Exp}(1)\), with \(E_1\) the exponential integral. Fading always costs: this is strictly below awgn_capacity() at the same \(\rho\), by Jensen’s inequality.

param snr:

Average signal-to-noise ratio \(\rho\), linear.

type snr:

float or np.ndarray

returns:

np.ndarray – Ergodic capacity in bits per channel use.

References

D. Tse, P. Viswanath, Fundamentals of Wireless Communication, Cambridge University Press, 2005, Section 5.4.

Examples

>>> capacity = rayleigh_ergodic_capacity(np.array([1.0, 10.0, 100.0]))
>>> print(np.round(capacity, 4))
[0.8603 2.9065 5.884 ]
>>> bool(np.all(capacity < awgn_capacity(np.array([1.0, 10.0, 100.0]))))
True

mimo_ergodic_capacity(n_tx: int, n_rx: int, snr: ndarray | float, *, n_realizations: int = 2000, rng: Generator | None = None) → ndarray

Ergodic capacity of an i.i.d. Rayleigh MIMO channel.

Signal Model

With equal power split across the transmit antennas and channel knowledge at the receiver only,

\[C = \mathbb{E}_{\mathbf{H}}\left[ \log_2 \det\left( \mathbf{I}_{N_r} + \frac{\rho}{N_t} \mathbf{H} \mathbf{H}^H \right) \right]\]

with \(\mathbf{H}\) of i.i.d. \(\mathcal{CN}(0,1)\) entries. The famous consequence is the multiplexing gain: at high SNR the capacity grows as \(\min(N_t, N_r) \log_2 \rho\), i.e. spatial dimensions buy rate, not just reliability.

No closed form exists for finite dimensions, so the expectation is estimated by Monte Carlo.

param n_tx:

Numbers of transmit and receive antennas \(N_t\), \(N_r\).

type n_tx:

int

param n_rx:

Numbers of transmit and receive antennas \(N_t\), \(N_r\).

type n_rx:

int

param snr:

Total signal-to-noise ratio \(\rho\), linear.

type snr:

float or np.ndarray

param n_realizations:

Channel draws averaged. Default 2000.

type n_realizations:

int, optional, keyword-only

param rng:

Random generator.

type rng:

numpy.random.Generator, optional, keyword-only

returns:

np.ndarray – Ergodic capacity in bits per channel use.

References

E. Telatar, “Capacity of multi-antenna Gaussian channels”, European Trans. Telecommunications 10(6), 1999; Tse & Viswanath, Chapter 8.

Examples

>>> rng = np.random.default_rng(0)
>>> siso = mimo_ergodic_capacity(1, 1, 10.0, n_realizations=20000, rng=rng)
>>> bool(abs(siso - rayleigh_ergodic_capacity(10.0)) < 0.05)   # 1x1 = SISO
True

outage_capacity(snr: float, outage: float = 0.01) → float

Outage capacity of a flat Rayleigh channel.

Signal Model

A slowly fading channel has no ergodic behaviour over one codeword: the rate either is or is not supported. The \(\varepsilon\)-outage capacity is the largest rate supported with probability \(1 - \varepsilon\). With \(|h|^2 \sim \mathrm{Exp}(1)\), \(\mathbb{P}\left[|h|^2 < t\right] = 1 - e^{-t}\), hence

\[C_{\varepsilon} = \log_2\left(1 - \rho \ln\left(1 - \varepsilon\right)\right)\]

It is far below the ergodic capacity at small \(\varepsilon\) – the price of not being able to average over the fading.

param snr:

Average signal-to-noise ratio \(\rho\), linear.

type snr:

float

param outage:

Outage probability \(\varepsilon\) in \((0, 1)\). Default 0.01.

type outage:

float, optional

returns:

float – Outage capacity in bits per channel use.

References

Tse & Viswanath, Fundamentals of Wireless Communication, 2005, Section 5.4.1.

Examples

>>> print(round(outage_capacity(100.0, outage=0.01), 4))
1.0036
>>> print(round(rayleigh_ergodic_capacity(100.0).item(), 4))
5.884

waterfilling(gains: ndarray, snr: float) → tuple[ndarray, float]

Optimal power allocation over parallel Gaussian channels.

Signal Model

Given \(K\) parallel channels of gain \(g_k\) – the subcarriers of an OFDM system, the eigenmodes of a MIMO channel – the capacity-achieving allocation solves

\[\max_{p_k \geq 0} \sum_k \log_2\left(1 + \rho g_k p_k\right) \quad \text{subject to} \quad \sum_k p_k = K\]

whose solution pours power onto the good channels up to a common water level \(\mu\):

\[p_k = \left(\mu - \frac{1}{\rho g_k}\right)^{+}\]

Channels below the water level get nothing: it is not worth spending power where the noise dominates.

param gains:

Channel power gains \(g_k\), non-negative.

type gains:

np.ndarray

param snr:

Reference signal-to-noise ratio \(\rho\) at unit power.

type snr:

float

returns:

tuple of (np.ndarray, float) – The allocation \(p_k\) (summing to \(K\)) and the resulting capacity in bits per channel use, per channel.

raises ValueError:

If a gain is negative or all gains are zero.

References

T. M. Cover, J. A. Thomas, Elements of Information Theory, 2nd ed., Wiley, 2006, Section 9.4 (parallel Gaussian channels).

Examples

>>> power, capacity = waterfilling(np.array([1.0, 1.0, 1.0]), snr=10.0)
>>> print(np.round(power, 6), round(capacity, 4))
[1. 1. 1.] 3.4594
>>> power, capacity = waterfilling(np.array([4.0, 1.0, 0.01]), snr=10.0)
>>> print(np.round(power, 4), round(capacity, 4) + 0.0)
[1.5375 1.4625 0.    ] 3.3105