Digital Communications Utilities

Digital modulation utilities (optic.comm.modulation)

grayCode(n)	Gray code generator.
grayMapping(M, constType)	Gray Mapping for digital modulations.
pamConst(M)	Generate a Pulse Amplitude Modulation (PAM) constellation.
qamConst(M)	Generate a Quadrature Amplitude Modulation (QAM) constellation.
pskConst(M)	Generate a Phase Shift Keying (PSK) constellation.
apskConst(M[, m1, phaseOffset])	Generate an Amplitude-Phase Shift Keying (APSK) constellation.
minEuclid(symb, const)	Find minimum Euclidean distance.
demap(indSymb, bitMap)	Contellation symbol index to bit sequence demapping.
modulateGray(bits, M, constType)	Modulate bit sequences to constellation symbol sequences (w/ Gray mapping).
demodulateGray(symb, M, constType)	Demodulate symbol sequences to bit sequences (w/ Gray mapping).
detector(r, σ2, constSymb[, px, rule])	Perform symbol detection using either the MAP (Maximum A Posteriori) or ML (Maximum Likelihood) rule.
softMapper(llr, M, constType[, prec])	Soft mapper for Gray-mapped modulation formats.
softEstimator(llr, bitMap, constSymb)	Estimates the mean and variance of the received symbols based on LLRs and the bit mapping.
mlse(y, h, constSymb)	Performs Maximum Likelihood Sequence Estimation (MLSE) using the Viterbi algorithm

apskConst(M, m1=None, phaseOffset=None)

Generate an Amplitude-Phase Shift Keying (APSK) constellation.

Parameters:

• M (int) – Constellation order.

• m1 (int) – Number of bits used to index the radii of the constellation.

• phaseOffset (float, optional) – Phase offset applied to each ring of the constellation [rad]. The default is pi / symbolsPerRing.

Returns:

const – APSK constellation

Return type:

np.array

Notes

symbols per ring is determined by \(m_2 = log_2(M) - m_1\), where \(m_1\) is the number of bits used to index the rings. If \(m_1\) is not provided, it is set based on M according to references [1] and [2].

References

[1] Z. Liu, et al “APSK Constellation with Gray Mapping,” IEEE Communications Letters, vol. 15, no. 12, pp. 1271-1273, 2011.

[2] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

demap(indSymb, bitMap)

Contellation symbol index to bit sequence demapping.

Parameters:

• indSymb (np.array of ints) – Indexes of received symbol sequence.

• bitMap ((M, log2(M)) np.array) – bit-to-symbol mapping.

Returns:

decBits – Sequence of demapped bits.

Return type:

np.array

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

demodulateGray(symb, M, constType)

Demodulate symbol sequences to bit sequences (w/ Gray mapping).

Hard demodulation is based on minimum Euclidean distance.

Parameters:

• symb (array of complex constellation symbols) – sequence of constellation symbols to be demodulated.

• M (int) – order of the modulation format.

• constType (string) – ‘qam’, ‘psk’, ‘apsk’, ‘pam’ or ‘ook’.

Returns:

sequence of demodulated bits.

Return type:

array of ints

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

detector(r, σ2, constSymb, px=None, rule='MAP')

Perform symbol detection using either the MAP (Maximum A Posteriori) or ML (Maximum Likelihood) rule.

Parameters:

• r (np.array) – The received signal.

• σ2 (float) – The noise variance.

• constSymb (np.array) – The constellation symbols.

• px (np.array, optional) – The prior probabilities of each symbol. If None, uniform priors are assumed.

• rule (str, optional) – The detection rule to use. Either ‘MAP’ (default) or ‘ML’.

Returns:

A tuple containing:

• np.array: The detected symbols.

• np.array: The indices of the detected symbols in the constellation.

Return type:

tuple

Notes

If px is None or rule is ‘ML’, uniform priors are assumed.

Given the received sample \(y = x + n\), where \(n\) is circular Gaussian noise with variance \(\sigma^2\), the maximum a posteriori (MAP) rule chooses the symbol with the largest posterior probability,

\[\hat{x}_{MAP} = \arg\max_{x \in \mathcal{X}}\, p(x \mid y) = \arg\max_{x \in \mathcal{X}}\left[-\frac{|y - x|^2}{\sigma^2} + \ln p(x)\right], \tag{1}\]

which minimizes the symbol error probability. For equiprobable symbols, the prior term \(\ln p(x)\) is constant and the MAP rule reduces to the maximum likelihood (ML) rule, i.e. the minimum Euclidean distance decision,

\[\hat{x}_{ML} = \arg\min_{x \in \mathcal{X}}\, |y - x|^2. \tag{2}\]

For non-uniform priors, as in probabilistically shaped constellations, the MAP decision regions shrink around the less likely symbols.

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

grayCode(n)

Gray code generator.

Parameters:

n (int) – length of the codeword in bits.

Returns:

code – list of binary strings of the gray code.

Return type:

list

Notes

In a Gray code, the binary words assigned to consecutive integers differ in a single bit. The \(i\)-th codeword is obtained from the binary representation of \(i\) as

\[g_i = i \oplus \left\lfloor i/2 \right\rfloor, \tag{1}\]

where \(\oplus\) denotes the bitwise exclusive OR. When the codewords are assigned to neighboring constellation points (Gray mapping), the most likely symbol errors, i.e. those between nearest neighbors, cause a single bit error.

grayMapping(M, constType)

Gray Mapping for digital modulations.

Parameters:

• M (int) – modulation order

• constType ('qam', 'psk', 'pam' or 'ook'.) – type of constellation.

Returns:

const – constellation symbols (sorted according their corresponding Gray bit sequence as integer decimal).

Return type:

np.array

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

minEuclid(symb, const)

Find minimum Euclidean distance.

Find closest constellation symbol w.r.t the Euclidean distance in the complex plane.

Parameters:

• symb (np.array) – Received constellation symbols.

• const (np.array) – Reference constellation.

Returns:

indexes of the closest constellation symbols.

Return type:

np.array of int

Notes

Each received symbol \(y\) is associated with the closest point of the constellation \(\mathcal{X} = \{x_0, \ldots, x_{M-1}\}\) in the Euclidean sense,

\[\hat{m} = \arg\min_{m}\, |y - x_m|^2. \tag{1}\]

For equiprobable symbols and additive white Gaussian noise, this is the maximum likelihood (ML) decision rule.

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

mlse(y, h, constSymb)

Performs Maximum Likelihood Sequence Estimation (MLSE) using the Viterbi algorithm

Parameters:

• y (array-like) – Received signal sequence

• h (array-like) – Channel impulse response

• constellation (array-like) – The available constellation symbols

Returns:

yMLSE – The MLSE decided output sequence

Return type:

ndarray

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

modulateGray(bits, M, constType)

Modulate bit sequences to constellation symbol sequences (w/ Gray mapping).

Parameters:

• bits (array of ints) – sequence of data bits.

• M (int) – order of the modulation format.

• constType (string) – ‘qam’, ‘psk’, ‘apsk’, ‘pam’ or ‘ook’.

Returns:

bits modulated to complex constellation symbols.

Return type:

array of complex constellation symbols

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

pamConst(M)

Generate a Pulse Amplitude Modulation (PAM) constellation.

Parameters:

M (int) – Number of symbols in the constellation. It must be an integer.

Returns:

1D PAM constellation.

Return type:

np.array

Notes

The \(M\) points of the pulse amplitude modulation (PAM) constellation are equally spaced and symmetric around zero,

\[\mathcal{X} = \left\{-(M-1), \ldots, -3, -1, 1, 3, \ldots, M-1\right\}, \tag{1}\]

with minimum distance \(d_{min} = 2\) and, for equiprobable symbols, average energy \(E_s = (M^2 - 1)/3\).

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

pskConst(M)

Generate a Phase Shift Keying (PSK) constellation.

Parameters:

M (int) – Number of symbols in the constellation. It must be a power of 2 positive integer.

Returns:

Complex M-PSK constellation.

Return type:

np.array

Notes

The \(M\) points of the phase shift keying (PSK) constellation lie on the unit circle, equally spaced in phase,

\[x_m = e^{j2\pi m/M}, \qquad m = 0, 1, \ldots, M-1, \tag{1}\]

so that all symbols have unit energy and the minimum distance is \(d_{min} = 2\sin(\pi/M)\).

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

qamConst(M)

Generate a Quadrature Amplitude Modulation (QAM) constellation.

Parameters:

M (int) – Number of symbols in the constellation. It must be a perfect square.

Returns:

const – Complex square M-QAM constellation.

Return type:

np.array

Notes

A square \(M\)-QAM constellation is the Cartesian product of two \(L\)-PAM constellations, \(L = \sqrt{M}\), one in the in-phase and one in the quadrature component,

\[\mathcal{X} = \left\{x_I + jx_Q \;:\; x_I, x_Q \in \{-(L-1), \ldots, -1, 1, \ldots, L-1\}\right\}, \tag{1}\]

with minimum distance \(d_{min} = 2\) and, for equiprobable symbols, average energy \(E_s = 2(M-1)/3\).

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

softEstimator(llr, bitMap, constSymb)

Estimates the mean and variance of the received symbols based on LLRs and the bit mapping.

Parameters:

• llr (ndarray of shape (numSymb, numBits)) – Log-likelihood ratios for each bit in the symbol.

• bitMap (ndarray of shape (M, numBits)) – Bit mapping of constellation points (binary matrix).

• constSymb (ndarray of shape (M,)) – Complex-valued constellation symbols.

Returns:

• softMean (ndarray of shape (numSymb,))

• softVar (ndarray of shape (numSymb,))

softMapper(llr, M, constType, prec=<class 'numpy.float32'>)

Soft mapper for Gray-mapped modulation formats.

Parameters:

• llr (1D numpy array) – Log-likelihood ratios (LLRs) of bits.

• M (int) – Modulation order.

• constType (str) – Type of constellation (‘qam’, ‘psk’, ‘pam’, ‘apsk’, or ‘ook’).

• prec (data type, optional) – Precision of the output (default is np.float32).

Returns:

• - softMean (1D numpy array) – Soft mean of the constellation symbols.

• - softVar (1D numpy array) – Soft variance of the constellation symbols.

Metrics for signal and performance characterization (optic.comm.metrics)

bert(Irx[, bitsTx, seed])	Calculate Bit Error Rate (BER) and Q-factor for optical communication using On-Off Keying (OOK).
fastBERcalc(rx, tx, M, constType[, px])	Monte Carlo BER/SER/SNR calculation.
calcLLR(rxSymb, σ2, constSymb, bitMap, px[, ...])	LLR calculation assuming a circular AGWN channel model.
calcExtrLLR(bitLLR, x, xMu, xNu, M, ...[, ...])	Calculate the extrinsic LLRs assuming an auxiliary Gaussian channel model.
monteCarloGMI(rx, tx, M, constType[, px, bitMap])	Monte Carlo based generalized mutual information (GMI) estimation.
monteCarloMI(rx, tx, M, constType[, px])	Monte Carlo based mutual information (MI) estimation.
calcMI(rx, tx, σ2, constSymb, pX)	Mutual information (MI) calculation for AWGN channels.
Qfunc(x)	Calculate function Q(x).
calcEVM(symb, M, constType[, symbTx])	Calculate error vector magnitude (EVM) metrics.
theoryBER(M, EbN0, constType)	Theoretical (approx.) bit error probability for PAM/QAM/PSK in AWGN channel.
theoryMI(M, constType, SNR[, pX, symmetry, ...])	Calculate mutual information for discrete input continuous output the memoryless AWGN channel (DCMC).
theoryGMI(M, constType, SNR[, bitMap, pX, ...])	Calculate generalized mutual information (GMI) for discrete input continuous output the memoryless AWGN channel (DCMC).
calcLinOSNR(Ns, Pin, α, Ls, OSNRin[, NF, ...])	Calculate the OSNR evolution in a multi-span fiber transmission system.

Qfunc(x)

Calculate function Q(x).

Parameters:

x (scalar) – function input.

Returns:

value of Q(x).

Return type:

scalar

Notes

The Gaussian Q-function is the tail probability of the standard normal distribution,

\[Q(x) = \frac{1}{\sqrt{2\pi}}\int_{x}^{\infty} e^{-u^2/2}\,du = \frac{1}{2}\,\mathrm{erfc}\left(\frac{x}{\sqrt{2}}\right). \tag{1}\]

It gives the probability that a zero-mean Gaussian random variable with unit variance exceeds \(x\), and appears in most error probability expressions for AWGN channels.

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

bert(Irx, bitsTx=None, seed=123)

Calculate Bit Error Rate (BER) and Q-factor for optical communication using On-Off Keying (OOK).

Parameters:

• Irx (np.array) – Received signal intensity values.

• bitsTx (np.array, optional) – Transmitted bit sequence. If not provided, a random bit sequence is generated.

• seed (int, optional) – Random seed for bit sequence generation when bitsTx is not provided.

Returns:

• BER (float) – Bit Error Rate, a measure of the number of erroneous bits relative to the total bits.

• Q (float) – Q-factor, a measure of signal quality in the communication system.

Notes

This function calculates the BER and Q-factor for an optical communication system using On-Off Keying (OOK) modulation. The received signal intensity Irx and an optional transmitted bit sequence bitsTx are required. If bitsTx is not provided, a random bit sequence is generated using the specified seed.

The following statistics are calculated for the received signal: - \(I_1\): The average value of the signal when the transmitted bit is 1. - \(I_0\): The average value of the signal when the transmitted bit is 0. - \(\sigma_1\): The standard deviation of the signal when the transmitted bit is 1. - \(\sigma_0\): The standard deviation of the signal when the transmitted bit is 0.

The optimal decision threshold Id and the Q-factor are calculated based on the signal statistics.

The function then applies the optimal decision rule to estimate the received bit sequence bitsRx. The Bit Error Rate (BER) is calculated by comparing bitsRx to bitsTx.

Assuming that the received signal, conditioned on each transmitted bit, is Gaussian distributed with mean \(I_b\) and standard deviation \(\sigma_b\), \(b \in \{0, 1\}\), the decision threshold that minimizes the bit error rate is well approximated by

\[I_D = \frac{\sigma_0 I_1 + \sigma_1 I_0}{\sigma_0 + \sigma_1}, \tag{1}\]

for which the error probabilities of the two symbols are equal. The quality of the signal is summarized by the Q-factor,

\[Q = \frac{I_1 - I_0}{\sigma_1 + \sigma_0}, \tag{2}\]

which is related to the bit error probability by

\[P_b = \frac{1}{2}\,\mathrm{erfc}\left(\frac{Q}{\sqrt{2}}\right) \approx \frac{e^{-Q^2/2}}{Q\sqrt{2\pi}}. \tag{3}\]

For example, \(Q = 6\) corresponds to \(P_b \approx 10^{-9}\).

References

[1] Agrawal, Govind P. Fiber-optic communication systems. John Wiley & Sons, 2012.

calcEVM(symb, M, constType, symbTx=None)

Calculate error vector magnitude (EVM) metrics.

Parameters:

• symb (np.array) – Sequence of noisy symbols.

• M (int) – Constellation order.

• constType (string) – Modulation type: ‘pam’, ‘qam’ or ‘psk’

• symbTx (np.array, optional) – Sequence of transmitted symbols (noiseless). The default is [].

Returns:

EVM – Squared error vector magnitude per signal dimension, i.e. the ratio between the power of the error vector and the power of the reference symbols (linear scale). The rms EVM in percent is 100*np.sqrt(EVM), and 10*np.log10(EVM) gives the EVM in dB.

Return type:

np.array

Notes

The error vector magnitude (EVM) measures the deviation of the received symbols \(y[k]\) from their references \(x[k]\). Both are normalized to unit average power, and the value returned for each signal dimension is the ratio of the power of the error vector to the power of the reference symbols,

\[\mathrm{EVM}^2 = \frac{\sum_k |y[k] - x[k]|^2}{\sum_k |x[k]|^2}. \tag{1}\]

If the transmitted symbols are provided, they are used as references (after a constant complex gain alignment for QAM and PSK); otherwise, the references are the closest constellation points (data-aided or decision-directed EVM). The rms EVM in percent is \(100\sqrt{\mathrm{EVM}^2}\), and, for an AWGN channel with data-aided references, \(\mathrm{EVM}^2 \approx 1/\mathrm{SNR}\).

References

[1] R. A. Shafik, et al, “On the error vector magnitude as a performance metric and comparative analysis”, em 2006 International Conference on Emerging Technologies, 2006, p. 27–31. doi: 10.1109/ICET.2006.335992.

[2] H. A. Mahmoud e H. Arslan, “Error vector magnitude to SNR conversion for nondata-aided receivers”, IEEE Transactions on Wireless Communications, vol. 8, nº 5, p. 2694–2704, 2009, doi: 10.1109/TWC.2009.080862.

calcExtrLLR(bitLLR, x, xMu, xNu, M, constSymb, bitMap, px=None, prec=<class 'numpy.float32'>)

Calculate the extrinsic LLRs assuming an auxiliary Gaussian channel model.

Parameters:

• bitLLR (np.array of shape (q*numSymb,)) – received bit LLRs

• x (np.array of shape (numSymb,)) – received symbols

• xMu (np.array of shape (numSymb,)) – mean of the received symbols

• xNu (np.array of shape (numSymb,)) – variance of the received symbols

• M (int) – modulation order

• constSymb (np.array of shape (M,)) – constellation symbols

• bitMap (np.array of shape (M, q)) – bit mapping of the constellation symbols

• px (np.array of shape (M,), optional) – prior probabilities of the constellation symbols, if None, uniform distribution is used

Returns:

LLRe – extrinsic LLRs for each bit

Return type:

np.array of shape (q*numSymb,)

calcLLR(rxSymb, σ2, constSymb, bitMap, px, maxLog=False)

LLR calculation assuming a circular AGWN channel model.

Parameters:

• rxSymb (np.array) – Received symbol sequence.

• σ2 (scalar) – Noise variance.

• constSymb ((M, 1) np.array) – Constellation symbols.

• bitMap ((M, log2(M)) np.array) – Bit-to-symbol mapping.

• px ((M, 1) np.array) – Prior symbol probabilities.

• maxLog (bool, optional) – If True, use the Max-Log approximation for LLR calculation.

Returns:

LLRs – sequence of calculated LLRs.

Return type:

np.array

Notes

For a transmitted symbol \(x\) drawn from the constellation \(\mathcal{X}\) with probability \(p(x)\), and an additive circular Gaussian noise with variance \(\sigma^2\), the log-likelihood ratio (LLR) of the \(i\)-th bit of the symbol, given the received sample \(y\), is

\[\Lambda_i(y) = \ln\frac{P(b_i = 0 \mid y)}{P(b_i = 1 \mid y)} = \ln\frac{\sum_{x \in \mathcal{X}_i^0} p(x)\, \exp\left(-|y-x|^2/\sigma^2\right)} {\sum_{x \in \mathcal{X}_i^1} p(x)\,\exp\left(-|y-x|^2/\sigma^2\right)}, \tag{1}\]

where \(\mathcal{X}_i^b\) is the subset of constellation points whose \(i\)-th bit is \(b\). The max-log approximation replaces each sum by its largest term,

\[\Lambda_i(y) \approx \max_{x \in \mathcal{X}_i^0}\left[-\frac{|y-x|^2}{\sigma^2} + \ln p(x)\right] - \max_{x \in \mathcal{X}_i^1}\left[-\frac{|y-x|^2}{\sigma^2} + \ln p(x)\right], \tag{2}\]

which avoids the evaluation of exponentials and logarithms at the cost of a small loss of accuracy at low SNR. Positive LLRs favor the bit value 0.

References

[1] A. Alvarado, T. Fehenberger, B. Chen, e F. M. J. Willems, “Achievable Information Rates for Fiber Optics: Applications and Computations”, Journal of Lightwave Technology, vol. 36, nº 2, p. 424–439, jan. 2018, doi: 10.1109/JLT.2017.2786351.

calcLinOSNR(Ns, Pin, α, Ls, OSNRin, NF=4.5, Fc=193100000000000.0, Bref=12500000000.0)

Calculate the OSNR evolution in a multi-span fiber transmission system.

Parameters:

• Ns (int) – Number of spans of fiber + EDFA.

• Pin (scalar) – Fiber launch power.

• α (scalar) – Fiber attenuation coefficient in dB/km.

• Ls (scalar) – Length of fiber spans in km.

• OSNRin (scalar) – OSNR at the input of the first span.

• NF (scalar, optional) – Noise figure of the EDFA amplifiers. The default is 4.5.

• Fc (scalar, optional) – Optical central frequency. The default is 193.1e12.

• Bref (scalar, optional) – Reference bandwidth for OSNR measurement. The default is 12.5e9.

Returns:

OSNR – OSNR values in dB at the output of each fiber span.

Return type:

np.array

Notes

In a chain of \(N_s\) fiber spans of length \(L_s\), each one followed by an EDFA whose gain \(G = \alpha L_s\) (in dB) compensates the span loss, the signal power at the output of each amplifier is equal to the launch power \(P_{in}\), whereas the ASE noise accumulates. Each amplifier adds the ASE noise power

\[P_{ASE} = 2N_{ASE}B_{ref}, \qquad N_{ASE} = (G-1)\,n_{sp}\,h\nu, \tag{1}\]

in both polarizations, within the reference bandwidth \(B_{ref}\) (usually 12.5 GHz, i.e. 0.1 nm at 1550 nm), where \(n_{sp} = (G\cdot NF - 1)/[2(G-1)]\) is obtained from the noise figure \(NF\) (see optic.models.devices.edfa). Since the gain of each amplifier compensates exactly the loss of the preceding span, the noise power at the output of the \(k\)-th amplifier is

\[P_{n,k} = G\,\frac{P_{n,k-1}}{G} + P_{ASE} = P_{n,k-1} + P_{ASE}, \qquad P_{n,0} = \frac{P_{in}}{\mathrm{OSNR}_{in}}, \tag{2}\]

where \(P_{n,0}\) is the noise at the input of the link, given by the input OSNR. Hence, \(P_{n,k} = P_{n,0} + kP_{ASE}\) and the OSNR after the \(k\)-th span is

\[\mathrm{OSNR}_k = \frac{P_{in}}{P_{n,0} + kP_{ASE}}. \tag{3}\]

For a negligible input noise, the OSNR decreases by 3 dB each time the number of spans doubles.

References

[1] J. G. Proakis; M. Salehi, Communication Systems Engineering, 2nd Edition. Pearson, 2002.

[2] R. -J. Essiambre, et al, “Capacity Limits of Optical Fiber Networks,” Journal of Lightwave Technology, vol. 28, no. 4, p. 662-701, 2010, doi: 10.1109/JLT.2009.2039464.

[3] R. Schober, P. Bayvel, e F. D. Pasquale, “Analytical model for the calculation of the optical signal-to-noise ratio (SNR) of WDM EDFA chains”, Optical and Quantum Electronics, vol. 31, no 3, p. 237–241. 1999, doi: 10.1023/A:1006948826091.

calcMI(rx, tx, σ2, constSymb, pX)

Mutual information (MI) calculation for AWGN channels.

Parameters:

• rx (np.array) – Received symbol sequence.

• tx (np.array) – Transmitted symbol sequence.

• σ2 (scalar) – Noise variance.

• constSymb ((M,) np.array) – Constellation symbols.

• pX ((M,) np.array) – Prob. mass function (p.m.f.) of the constellation symbols.

Returns:

Estimated mutual information.

Return type:

scalar

Notes

The mutual information is computed as

\[I(X;Y) = H(X) - H(X \mid Y), \qquad H(X) = -\sum_{x \in \mathcal{X}} p(x)\log_2 p(x), \tag{1}\]

where the conditional entropy is estimated by Monte Carlo averaging over the \(N\) pairs of transmitted and received symbols \((x_k, y_k)\),

\[H(X \mid Y) \approx -\frac{1}{N}\sum_{k=1}^{N} \log_2\frac{p(y_k \mid x_k)\,p(x_k)} {\sum_{x \in \mathcal{X}} p(y_k \mid x)\,p(x)}. \tag{2}\]

The channel transition probability is that of the AWGN channel with noise variance \(\sigma^2\), i.e. \(p(y \mid x) \propto \exp\left(-|y-x|^2/\sigma^2\right)\) for complex symbols and \(p(y \mid x) \propto \exp\left[-(y-x)^2/(2\sigma^2)\right]\) for real symbols (the normalization constants cancel out in Eq. (2)). When the actual channel is not AWGN, the result is a lower bound on the mutual information of the channel, achievable by a receiver designed for the AWGN channel.

fastBERcalc(rx, tx, M, constType, px=None)

Monte Carlo BER/SER/SNR calculation.

Parameters:

• rx (np.array) – Received symbol sequence.

• tx (np.array) – Transmitted symbol sequence.

• M (int) – Modulation order.

• constType (string) – Modulation type: ‘qam’, ‘psk’, ‘pam’ or ‘ook’.

• px ((M, 1) np.array) – Prior symbol probabilities.

Returns:

• BER (np.array) – Bit-error-rate.

• SER (np.array) – Symbol-error-rate.

• SNR (np.array) – Estimated SNR from the received constellation.

Notes

Before the error counting, the received symbols \(y[k]\) are aligned with the transmitted ones \(x[k]\): for QAM and PSK, a constant complex gain \(\hat{h} = \frac{1}{N}\sum_k x[k]/y[k]\) compensates a residual amplitude scaling and phase rotation, and both sequences are normalized to unit power. The signal-to-noise ratio is then estimated as

\[\widehat{\mathrm{SNR}} = \frac{\sum_k |x[k]|^2}{\sum_k |y[k] - x[k]|^2}. \tag{1}\]

The received symbols are demodulated by minimum Euclidean distance hard decisions and Gray demapping, and the error rates are estimated by counting,

\[\mathrm{BER} = \frac{N_{b,err}}{N\log_2 M}, \qquad \mathrm{SER} = \frac{N_{s,err}}{N}, \tag{2}\]

where \(N_{b,err}\) is the number of wrong bits and \(N_{s,err}\) is the number of symbols with at least one wrong bit.

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

monteCarloGMI(rx, tx, M, constType, px=None, bitMap=None)

Monte Carlo based generalized mutual information (GMI) estimation.

Parameters:

• rx (np.array) – Received symbol sequence.

• tx (np.array) – Transmitted symbol sequence.

• M (int) – Modulation order.

• constType (string) – Modulation type: ‘qam’ or ‘psk’

• px ((M, 1) np.array) – Prior symbol probabilities. The default is [].

• bitMap ((M, b) np.array) – Bit mapping matrix. The default is None.

Returns:

• GMI (np.array) – Generalized mutual information values.

• NGMI (np.array) – Normalized mutual information.

Notes

The generalized mutual information (GMI) is an achievable rate, in bits per symbol, for systems with bit-wise decoding, where a binary FEC decoder processes the bit LLRs independently. For a constellation with \(m = \log_2 M\) bits per symbol, it is estimated from \(N\) transmitted symbols as

\[\mathrm{GMI} \approx H(X) - \sum_{i=1}^{m}\frac{1}{N}\sum_{k=1}^{N} \log_2\left(1 + e^{(2b_{k,i} - 1)\Lambda_{k,i}}\right), \tag{1}\]

where \(H(X) = -\sum_x p(x)\log_2 p(x)\) is the entropy of the transmitted symbols, \(b_{k,i}\) is the \(i\)-th bit of the \(k\)-th symbol and \(\Lambda_{k,i}\) is its LLR (see calcLLR), computed assuming an AWGN channel whose noise variance is estimated from the data. Each term of the inner sum is equal to \(-\log_2 P(b_{k,i} \mid y_k)\), the information that is still missing about the transmitted bit after the observation of \(y_k\). The normalized GMI,

\[\mathrm{NGMI} = \frac{\mathrm{GMI}}{H(X)}, \tag{2}\]

is the fraction of the transmitted information that is recovered, and is commonly used as a threshold for soft-decision FEC decoders.

References

[1] A. Alvarado, T. Fehenberger, B. Chen, e F. M. J. Willems, “Achievable Information Rates for Fiber Optics: Applications and Computations”, Journal of Lightwave Technology, vol. 36, nº 2, p. 424–439, jan. 2018, doi: 10.1109/JLT.2017.2786351.

monteCarloMI(rx, tx, M, constType, px=None)

Monte Carlo based mutual information (MI) estimation.

Parameters:

• rx (np.array) – Received symbol sequence.

• tx (np.array) – Transmitted symbol sequence.

• M (int) – Modulation order.

• constType (string) – Modulation type: ‘qam’ or ‘psk’

• px ((M, 1) np.array) – p.m.f. of the constellation symbols. The default is [].

Returns:

MI – Estimated MI values.

Return type:

np.array

Notes

The mutual information \(I(X;Y)\) between the transmitted and the received symbols is the maximum achievable rate, in bits per symbol, for a given input distribution \(p(x)\) with symbol-wise decoding. After aligning and normalizing the received symbols with respect to the transmitted ones, the noise variance \(\sigma^2\) is estimated from the data, and \(I(X;Y)\) is estimated by assuming an AWGN channel with this noise variance (see calcMI).

References

[1] A. Alvarado, T. Fehenberger, B. Chen, e F. M. J. Willems, “Achievable Information Rates for Fiber Optics: Applications and Computations”, Journal of Lightwave Technology, vol. 36, nº 2, p. 424–439, jan. 2018, doi: 10.1109/JLT.2017.2786351.

theoryBER(M, EbN0, constType)

Theoretical (approx.) bit error probability for PAM/QAM/PSK in AWGN channel.

Parameters:

• M (int) – Modulation order.

• EbN0 (scalar) – Signal-to-noise ratio (SNR) per bit in dB.

• constType (string) – Modulation type: ‘pam’, ‘qam’ or ‘psk’

Returns:

Pb – Theoretical probability of bit error.

Return type:

scalar

Notes

The values of error probability obtained with this function are good approximations for moderate to high SNR regime (see [1]). All cases assume Gray mapped constellations. For low SNR values and high constellation cardinalities (\(P_b\)>1e-1), the results should underestimate the real error probability.

For an average energy per bit to noise power spectral density ratio \(E_b/N_0\) and \(k = \log_2 M\) bits per symbol, the approximations used are:

• Square M-QAM, with \(L = \sqrt{M}\):

  \[P_b \approx \frac{2\left(1 - 1/L\right)}{\log_2 L}\, Q\left(\sqrt{\frac{3\log_2 L}{L^2 - 1}\,\frac{2E_b}{N_0}}\right). \tag{1}\]

• M-PSK:

  \[P_b \approx \frac{2}{k}\, Q\left(\sqrt{\frac{2kE_b}{N_0}}\,\sin\frac{\pi}{M}\right). \tag{2}\]

• M-PAM:

  \[P_b \approx \frac{2(M-1)}{kM}\, Q\left(\sqrt{\frac{6\log_2 M}{M^2 - 1}\,\frac{E_b}{N_0}}\right). \tag{3}\]

Here \(Q(\cdot)\) is the Gaussian Q-function (see Qfunc). The PSK and PAM expressions approximate the symbol error probability and assume that, with Gray mapping, each symbol error causes a single bit error, so that \(P_b \approx P_s/k\).

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

theoryGMI(M, constType, SNR, bitMap=None, pX=None, lim=None, tol=0.001)

Calculate generalized mutual information (GMI) for discrete input continuous output the memoryless AWGN channel (DCMC).

Parameters:

• M (int) – Number of symbols in the constellation.

• constType (str) – Type of constellation (‘qam’, ‘psk’).

• SNR (float) – Signal-to-noise ratio in dB.

• bitMap (array_like) – Bit mapping of the constellation symbols (shape: M x log2(M)).

• pX (array_like, optional) – Probability of each transmitted symbol (default is None).

• lim (float, optional) – Limit for numerical integration (default is np.inf).

• tol (float, optional) – Tolerance for numerical integration error (default is 1e-3).

Returns:

Theoretical GMI for the given parameters.

Return type:

float

Notes

The generalized mutual information of a bit-interleaved coded modulation (BICM) system with \(m = \log_2 M\) bits per symbol is the sum of the mutual informations between each bit \(B_i\) of the symbol and the channel output,

\[\mathrm{GMI} = \sum_{i=1}^{m} I(B_i; Y) = \sum_{i=1}^{m}\left[H(B_i) - H(B_i \mid Y)\right], \tag{1}\]

where the channel law of each bit is the mixture of the AWGN transition probabilities of the constellation points that carry it,

\[p(y \mid B_i = b) = \frac{1}{P(B_i = b)} \sum_{x \in \mathcal{X}_i^b} p(x)\,p(y \mid x). \tag{2}\]

The conditional entropies \(H(B_i \mid Y)\) are evaluated by numerical integration. The GMI depends on the bit mapping and is upper bounded by the mutual information of the constellation (see theoryMI), being close to it for Gray mapping.

References

[1] A. Alvarado, T. Fehenberger, B. Chen, e F. M. J. Willems, “Achievable Information Rates for Fiber Optics: Applications and Computations”, Journal of Lightwave Technology, vol. 36, nº 2, p. 424–439, jan. 2018, doi: 10.1109/JLT.2017.2786351.

theoryMI(M, constType, SNR, pX=None, symmetry=True, lim=inf, tol=0.001)

Calculate mutual information for discrete input continuous output the memoryless AWGN channel (DCMC).

Parameters:

• M (int) – Number of symbols in the constellation.

• constType (str) – Type of constellation (‘qam’, ‘psk’).

• SNR (float) – Signal-to-noise ratio in dB.

• pX (array_like, optional) – Probability of each transmitted symbol (default is None).

• symmetry (bool, optional) – Flag to exploit rotational symmetry of the constellation (default is True).

• lim (int, optional) – Limit for numerical integration (default is np.inf).

• tol (float, optional) – Tolerance for numerical integration error (default is 1e-3).

Returns:

Mutual information for the given parameters.

Return type:

float

Notes

For a discrete-input continuous-output memoryless channel (DCMC) with additive white Gaussian noise, the mutual information is

\[I(X;Y) = H(X) - \sum_{x \in \mathcal{X}} p(x)\int p(y \mid x) \log_2\frac{p(y)}{p(y \mid x)\,p(x)}\,dy, \tag{1}\]

with \(p(y) = \sum_{x} p(y \mid x)\,p(x)\) and, for complex constellations normalized to unit average energy,

\[p(y \mid x) = \frac{1}{2\pi\sigma^2} \exp\left(-\frac{|y-x|^2}{2\sigma^2}\right), \qquad \sigma^2 = \frac{1}{2\,\mathrm{SNR}}, \tag{2}\]

where \(\sigma^2\) is the noise variance per real dimension (for PAM, a real-valued Gaussian with \(\sigma^2 = 1/\mathrm{SNR}\) is used). The integral in Eq. (1) is evaluated numerically. Since the integrand of constellation points with the same magnitude is identical for rotationally symmetric constellations, the integral may be computed once for each group of such points (symmetry=True).

References

[1] A. Alvarado, T. Fehenberger, B. Chen, e F. M. J. Willems, “Achievable Information Rates for Fiber Optics: Applications and Computations”, Journal of Lightwave Technology, vol. 36, nº 2, p. 424–439, jan. 2018, doi: 10.1109/JLT.2017.2786351.

Sources of discrete sequences (optic.comm.sources)

bitSource(param)	Generate a sequence of bits of length nBits either as a random bit sequence or a pseudo-random binary sequence (PRBS).
prbsGenerator([order, length, seed])	Generate a Pseudo-Random Binary Sequence (PRBS) of the given order.
symbolSource(param)	Generate a random symbol sequence based on the specified modulation scheme, order, and pmf.
cazacSequence(N[, M])	Generate a CAZAC (Zadoff-Chu) sequence of length N.

bitSource(param)

Generate a sequence of bits of length nBits either as a random bit sequence or a pseudo-random binary sequence (PRBS).

Parameters:

param (optic.utils.parameters) –

Parameters for the bit source:

• nBits : int, optional. The number of bits in the sequence. [default: 1000]

• mode : str, optional. The mode of the bit generation. If ‘random’, a sequence of random bits is generated. If ‘prbs’, a pseudo-random binary sequence (PRBS) is generated [default: ‘random’].

• order : int, optional. The order of the PRBS generator. Only used if mode is ‘prbs’ [default: 23].

• seed : int, optional. The seed for the random number generator. Only applicable when mode is ‘random’ [default: None].

Returns:

bits – An array of bits of length nBits, either randomly generated or from a PRBS.

Return type:

np.array

References

[1] Wikipedia, “Pseudorandom binary sequence,” https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence

[2] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

cazacSequence(N, M=1)

Generate a CAZAC (Zadoff-Chu) sequence of length N.

Parameters:

• N (int) – The length of the CAZAC sequence.

• M (int, optional) – The root of the CAZAC sequence. Default is 1.

Returns:

sequence – A NumPy array containing the generated CAZAC sequence.

Return type:

np.array

Notes

Constant amplitude zero autocorrelation (CAZAC) sequences are used for synchronization and channel estimation. This function generates the Zadoff-Chu sequence of length \(N\) and root \(M\),

\[x[n] = \exp\left[-j\pi M\frac{n(n+c_f)}{N}\right], \qquad n = 0, 1, \ldots, N-1, \tag{1}\]

where \(M\) and \(N\) are coprime, and \(c_f = N \bmod 2\), i.e. the exponent is \(n(n+1)\) for odd \(N\) and \(n^2\) for even \(N\). The sequence has constant amplitude, \(|x[n]| = 1\), and an ideal periodic autocorrelation,

\[\sum_{n=0}^{N-1} x[n]\, x^*[(n + k) \bmod N] = N\,\delta[k], \tag{2}\]

so that its cyclically shifted versions are orthogonal.

References

[1] D. Chu, “Polyphase codes with good periodic correlation properties (Corresp.),” IEEE Transactions on Information Theory, 18 (4), pp. 531-532, 1972.

prbsGenerator(order=23, length=None, seed=1)

Generate a Pseudo-Random Binary Sequence (PRBS) of the given order.

Parameters:

• order (int) – The order of the PRBS sequence. Supported orders are 7, 9, 11, 13, 15, 23, 31.

• length (int, optional) – The length of the PRBS sequence to generate. If not specified, the length is \(2^{order} - 1\).

• seed (int, optional) – The seed for the linear feedback shift register (LFSR). Default is 1.

Returns:

bits – A NumPy array of bits representing the PRBS sequence.

Return type:

np.array

Notes

A pseudo-random binary sequence (PRBS) of order \(n\) is generated by an \(n\)-stage linear feedback shift register (LFSR). At each clock cycle, the register outputs its most significant bit and is shifted by one position, and the feedback bit, computed as the XOR of the register taps given by the generator polynomial, is inserted. The bit sequence satisfies the linear recursion over GF(2) defined by the generator polynomial, for instance

\[b_k = b_{k-7} \oplus b_{k-6} \quad \text{for PRBS-7}, \qquad g(x) = x^7 + x^6 + 1. \tag{1}\]

The generator polynomials used are \(x^7 + x^6 + 1\), \(x^9 + x^5 + 1\), \(x^{11} + x^9 + 1\), \(x^{13} + x^{12} + x^2 + x + 1\), \(x^{15} + x^{14} + 1\), \(x^{23} + x^{18} + 1\) and \(x^{31} + x^{28} + 1\). Since these polynomials are primitive, the output is a maximal-length sequence (m-sequence): it repeats with period \(2^n - 1\), contains \(2^{n-1}\) ones and \(2^{n-1} - 1\) zeros per period, and its periodic autocorrelation (in bipolar form) is equal to \(2^n - 1\) at zero lag and to \(-1\) at any other lag, which makes it resemble a random sequence.

References

[1] Wikipedia, “Pseudorandom binary sequence,” https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence

symbolSource(param)

Generate a random symbol sequence based on the specified modulation scheme, order, and pmf.

Parameters:

param (optic.utils.parameters object) –

Parameters of the symbol source:

• nSymbols : int, optional. The number of symbols to generate. [default: 1000]

• M : int, optional. The modulation order, defining the size of the constellation. [default: 4]

• constType : str, optional. The type of modulation scheme. Supported types are ‘qam’, ‘pam’, ‘psk’, and ‘apsk’. [default: ‘qam’].

• dist : str, optional. The probability distribution for generating symbols. Options are ‘uniform’ or ‘maxwell-boltzmann’ [default: ‘uniform’].

• shapingFactor : float, optional. The shaping factor applied when dist is ‘maxwell-boltzmann’. Controls the shaping of the constellation points. [default: 0.0].

• px : array-like, optional. Custom probability distribution for the constellation points. If None, the distribution is determined by dist. [default: None].

• seed : int, optional. Seed for the random number generator to ensure reproducibility [default: None].

Returns:

symbols – A NumPy array containing the generated symbols from the specified constellation, based on the given probability distribution.

Return type:

np.array

Notes

The function generates symbols using a specified modulation scheme following a uniform or Maxwell-Boltzmann distribution to the constellation points. The Maxwell-Boltzmann distribution is shaped by the shapingFactor. If a custom probability distribution px is provided, it will override the default distribution.

If the constType is set to ‘qam’, ‘pam’, ‘psk’, or ‘apsk’, the corresponding constellation is used. Custom modulation schemes are not supported.

The seed parameter ensures the same sequence of symbols is generated across different runs when set.

With the Maxwell-Boltzmann distribution, the probability of each constellation point decreases exponentially with its energy,

\[p(x) = \frac{e^{-\lambda|x|^2}}{\sum_{x' \in \mathcal{X}} e^{-\lambda|x'|^2}}, \tag{1}\]

where \(\lambda \geq 0\) is the shaping factor (\(\lambda = 0\) gives the uniform distribution). This probabilistic constellation shaping reduces the average energy for a given entropy \(H(X)\), approaching the Shannon capacity of the AWGN channel more closely than uniform signaling. In all cases, the constellation is normalized to unit average energy, \(\sum_x p(x)|x|^2 = 1\).

References

[1] Junho Cho and Peter J. Winzer, “Probabilistic Constellation Shaping for Optical Fiber Communications,” J. Lightwave Technol. 37, 1590-1607 (2019).

OFDM utilities (optic.comm.ofdm)

hermit(V)	Hermitian simmetry block.
zeroPad(x, L)	Pad an input array with zeros on both sides.
calcSymbolRate(M, Rb, Nfft, Np, G, hermitSym)	Calculate the symbol rate of a given OFDM configuration.
modulateOFDM(symb[, param])	Modulate OFDM signal.
demodulateOFDM(sig[, param])	Demodulate OFDM signal.

calcSymbolRate(M, Rb, Nfft, Np, G, hermitSym)

Calculate the symbol rate of a given OFDM configuration.

Parameters:

• M (int) – Constellation order

• Rb (float) – Bit rate

• Nfft (int) – Size of FFT

• Np (int) – number of pilot subcarriers

• G (scalar) – cyclic prefix length

• hermitSym (boolean) – True: Real OFDM symbols / False: Complex OFDM symbols

Returns:

Rs – OFDM symbol rate

Return type:

scalar

Notes

Each OFDM frame of \(N_{FFT} + G\) samples, where \(G\) is the length of the cyclic prefix, carries \(N_d\) data symbols of \(\log_2 M\) bits each. To transport the bit rate \(R_b\), the sampling (symbol) rate of the OFDM signal must be

\[R_s = \frac{R_b}{\dfrac{N_d}{N_{FFT} + G}\,\log_2 M}, \tag{1}\]

where \(N_d = N_{FFT} - N_p\) for complex OFDM and \(N_d = N_{FFT}/2 - 1 - N_p\) with Hermitian symmetry, \(N_p\) being the number of pilot subcarriers.

demodulateOFDM(sig, param=None)

Demodulate OFDM signal.

Parameters:

• sig (np.np.array) – Complex-valued array representing the OFDM signal sequence received at one sample per symbol.

• param (optic.utils.parameters object, optional) –

  Parameters for OFDM demodulation.

  • param.Nfft : scalar, optional. Size of the FFT [default: 512].

  • param.G : scalar, optional. Cyclic prefix length [default: 4].

  • param.hermitSymmetry : bool, optional. If True, indicates real OFDM symbols; if False, indicates complex OFDM symbols [default: False].

  • param.pilot : complex-valued scalar, optional. Pilot symbol [default: 1 + 1j].

  • param.pilotCarriers : np.array, optional. Indexes of pilot subcarriers [default: an empty array].

  • param.nullCarriers : np.array, optional. Indexes of null subcarriers [default: an empty array].

  • param.returnChannel : bool, optional. If True, return the estimated channel [default: False].

Returns:

If returnChannel is False, returns a complex-valued array representing the demodulated symbols sequence received. If returnChannel is True, returns a tuple containing the demodulated symbols sequence received and the estimated channel.

Return type:

np.array or tuple

Notes

• The input signal must be sampled at one sample per symbol.

• This function performs demodulation of the OFDM signal according to the provided parameters, including channel estimation and single tap equalization.

For each received frame, the cyclic prefix is removed and the subcarrier symbols are obtained with a DFT,

\[Y[k] = \frac{1}{\sqrt{N_{FFT}}}\sum_{n=0}^{N_{FFT}-1} y[n]\, e^{-j2\pi kn/N_{FFT}} = H[k]X[k] + N[k], \tag{1}\]

where \(H[k]\) is the frequency response of the channel at the \(k\)-th subcarrier. If pilot subcarriers are used, the channel is estimated at the pilot positions \(k_p\), as \(\hat{H}[k_p] = Y[k_p]/X_p\), where \(X_p\) is the pilot symbol. The magnitude and the phase of \(\hat{H}\) are then linearly interpolated over all subcarriers and averaged over the frames, and the data symbols are recovered by one-tap (zero-forcing) equalization,

\[\hat{X}[k] = \frac{Y[k]}{\hat{H}[k]}. \tag{2}\]

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

hermit(V)

Hermitian simmetry block.

Parameters:

V (complex-valued np.array) – input array

Returns:

Vh – vector with hermitian simmetry

Return type:

complex-valued np.array

Notes

The discrete Fourier transform of a real-valued sequence of length \(N\) satisfies the Hermitian symmetry \(X[N-k] = X^*[k]\). To generate a real-valued OFDM signal (e.g. for intensity modulation), the vector \(V\) of \(N_s\) subcarrier symbols is arranged as

\[V_h = \left[0,\; V[0], \ldots, V[N_s-1],\; 0,\; V^*[N_s-1], \ldots, V^*[0]\right], \tag{1}\]

with length \(2N_s + 2\), where the DC and the Nyquist subcarriers are set to zero. The inverse DFT of \(V_h\) is then real-valued.

modulateOFDM(symb, param=None)

Modulate OFDM signal.

Parameters:

• symb (np.np.array) – Complex-valued array of modulation symbols representing the symbols sequence to be transmitted.

• param (optic.utils.parameters object, optional) –

  Parameters for OFDM modulation.

  • param.Nfft : scalar, optional. Size of the FFT. [default: 512].

  • param.G : scalar, optional. Cyclic prefix length. [default: 4].

  • param.hermitSymmetry : bool, optional. If True, indicates real OFDM symbols; if False, indicates complex OFDM symbols. [default: False].

  • param.pilot : complex-valued scalar, optional. Pilot symbol. [default: 1 + 1j].

  • param.pilotCarriers : np.array, optional. Indexes of pilot subcarriers. [default: empty array].

  • param.nullCarriers : np.array, optional. Indexes of null subcarriers. [default: empty array].

  • param.SpS : int, optional. Oversampling factor. [default: 2].

Returns:

Complex-valued array representing the OFDM symbols sequence transmitted.

Return type:

np.array

Notes

In orthogonal frequency division multiplexing (OFDM), the data symbols are transmitted in parallel over \(N_{FFT}\) orthogonal subcarriers. Each OFDM frame is built by placing the data symbols on the data subcarriers, the pilot symbol on the pilot subcarriers, and zeros on the null subcarriers, which yields the vector \(X[k]\), \(k = 0, \ldots, N_{FFT}-1\) (with Hermitian symmetry, see hermit, the time-domain signal is real-valued). The time-domain samples are obtained with an inverse DFT,

\[x[n] = \frac{1}{\sqrt{N_{FFT}}}\sum_{k=0}^{N_{FFT}-1} X[k]\, e^{j2\pi kn/N_{FFT}}, \qquad n = 0, \ldots, N_{FFT}-1. \tag{1}\]

Eq. (1) corresponds to SpS = 1; oversampling by SpS is implemented by zero padding the spectrum before the inverse DFT, with the normalization adjusted to preserve the signal power. Finally, a cyclic prefix is inserted by copying the last \(G\) samples of the frame to its beginning,

\[\tilde{x}[n] = x[(n - G) \bmod N_{FFT}], \qquad n = 0, \ldots, N_{FFT} + G - 1. \tag{2}\]

As long as the channel memory is shorter than the cyclic prefix, the linear convolution with the channel impulse response becomes a circular convolution within each frame, so that each subcarrier experiences a single complex gain \(H[k]\).

References

[1] Proakis, J. G., & Salehi, M. Digital Communications (5th Edition). McGraw-Hill Education, 2008.

zeroPad(x, L)

Pad an input array with zeros on both sides.

Parameters:

• x (array_like) – Input array to be padded.

• L (int) – Number of zeros to pad on each side of the array.

Returns:

padded_array – Padded array with zeros added at the beginning and end.

Return type:

np.array

Notes

This function pads the input array x with L zeros on both sides, effectively increasing its length by 2*L.

Forward error correction (FEC) utilities (optic.comm.fec)

par2gen(H)	Convert a binary parity-check matrix H into a systematic generator matrix G over GF(2).
gaussElim(M)	Perform Gaussian elimination over GF(2) to reduce a binary matrix to row echelon form.
encoder(G, bits[, systematic])	Encode binary sequences using a generator matrix over GF(2).
encodeDVBS2(bits, A)	Encode multiple binary sequences using a DVB-S2 LDPC parity-check matrix.
encodeTriang(bits, P1, P2)	Encode binary sequences using two parity matrices for LDPC encoding.
sumProductAlgorithm(llrs, checkNodes, ...[, ...])	Performs belief propagation decoding using the sum-product algorithm (SPA) for multiple codewords.
minSumAlgorithm(llrs, checkNodes, varNodes, ...)	Performs belief propagation decoding using the Min-Sum Algorithm (MSA) for multiple codewords.
encodeLDPC(bits, param)	Encode binary sequences using a parity-check matrix of a LDPC code.
decodeLDPC(llrs, param)	Decode multiple LDPC codewords using the belief propagation algorithms.
writeAlist(H, filename)	Save a binary parity-check matrix H (numpy array) to ALIST format.
readAlist(filename)	Read an ALIST file and reconstruct the binary parity-check matrix H.
triangularize(H)	Convert binary matrix H into lower-triangular form using only row and column permutations.
triangP1P2(H)	Convert a binary parity-check matrix H into a lower-triangular form and extract matrices P1 and P2.
inverseMatrixGF2(A)	Invert a square binary matrix over GF(2) using Gauss-Jordan elimination.
plotBinaryMatrix(H)	Plot the binary matrix H with dots at positions where H[i,j] = 1.
parseAlist(path)	Parse an LDPC ALIST file and return basic parameters.
summarizeAlistFolder(folderPath)	Scan a folder for .alist files and print summary table.
hammingParityCheckMatrix(m[, extended])	Generate the parity-check matrix H for a Hamming code with m check bits .
encodeHamming(bits, param)	Encode binary sequences using a Hamming code.

decodeLDPC(llrs, param)

Decode multiple LDPC codewords using the belief propagation algorithms.

Parameters:

• llrs (np.array of shape (n, numCodewords)) – Array of log-likelihood ratios (LLRs) for each bit of the received codewords. Codewords are assumed to be disposed in columns.

• param (optic.utils.parameters object) –

  Object containing the following attributes:

  • H : Sparse binary parity-check matrix of shape (m, n) [default: None].

  • maxIter : Maximum number of iterations for belief propagation [default: 25].

  • alg : Decoding algorithm to use (‘SPA’ for sum-product or ‘MSA’ for min-sum) [default: ‘SPA’].

  • prec : Numerical precision to use in computations (default is np.float32) [default: np.float32].

Returns:

• decodedBits (np.array of shape (n, numCodewords)) – Array of decoded bits for each codeword.

• outputLLRs (np.array of shape (n, numCodewords)) – Array of updated log-likelihood ratios (LLRs) after decoding.

encodeDVBS2(bits, A)

Encode multiple binary sequences using a DVB-S2 LDPC parity-check matrix.

Parameters:

• bits (np.array of shape (k, N)) – Binary input sequences to be encoded. Each column represents a bit sequence of length \(k\).

• A (np.array of shape (m, k)) – Matrix corresponding to the first \(k\) columns of the parity-check matrix \(H\).

Returns:

codewords – Binary encoded codewords. Each column is a codeword of length \(n\) corresponding to the respective input bit sequence.

Return type:

np.array of shape (n, N)

encodeHamming(bits, param)

Encode binary sequences using a Hamming code.

Parameters:

• bits (np.array of shape (k, N)) – Binary input sequences to encode. Each column is a bit sequence of length \(k\).

• param (optic.utils.parameters object) –

  Object containing the following attributes:

  • m : Number of check bits for the Hamming code [default: 3].

  • extended : boolean indication of whether to use the extended Hamming code [default: False].

Returns:

codewords – Encoded binary codewords. Each column corresponds to a codeword of length \(n\) resulting from applying the Hamming encoding to the respective input bit sequence.

Return type:

np.array of shape (n, N)

References

[1] R. W. Hamming, “Error detecting and error correcting codes,” Bell System Technical Journal, vol. 29, no. 2, pp. 147-160, April 1950.

[2] S. Lin and D. J. Costello, “Error Control Coding: Fundamentals and Applications,” 2nd Edition, Pearson, 2004.

encodeLDPC(bits, param)

Encode binary sequences using a parity-check matrix of a LDPC code.

Parameters:

• bits (np.array of shape (k, N)) – Binary input sequences to be encoded. Each column is a bit sequence of length \(k\) bits.

• param (optic.utils.parameters object) –

  Object containing the following attributes:

  • mode : Mode of operation (‘DVBS2’, ‘IEEE_802.11nD2’, or ‘AR4JA’)[default: ‘DVBS2’].

  • n : Codeword length \(n\) [default: 64800].

  • R : Code rate \(R\) [default: ‘4/5’].

  • H : Binary parity-check matrix \(H\) of shape \((n - k, n)\) [default: None].

  • G : Binary generator matrix \(G\) of shape \((k, n)\) [default: None].

  • systematic : boolean indicator if the generator matrix is assumed to be in systematic form. [default: True]

  • P1 : Matrix of shape (m, k) used for encoding in triangular mode [default: None].

  • P2 : Matrix of shape (m, m) used for encoding in triangular mode [default: None].

  • path : Path to the folder containing ALIST files for the specified mode [default: None].

Returns:

codewords – Binary encoded codewords. Each column is a codeword of length \(n\) corresponding to the respective input bit sequence.

Return type:

np.array of shape (n, N)

References

[1] T. J. Richardson and R. L. Urbanke, “Efficient encoding of low-density parity-check codes,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 638-656, Feb 2001.

encodeTriang(bits, P1, P2)

Encode binary sequences using two parity matrices for LDPC encoding.

Parameters:

• bits (np.array of shape (k, N)) – Binary input sequences. Each column is a bit sequence to be encoded.

• P1 (np.array of shape (m1, k)) – First parity matrix.

• P2 (np.array of shape (m2, k)) – Second parity matrix.

Returns:

codewords – Encoded codewords, one per column.

Return type:

np.array of shape (k + m1 + m2, N)

References

[1] T. J. Richardson and R. L. Urbanke, “Efficient encoding of low-density parity-check codes,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 638-656, Feb 2001.

encoder(G, bits, systematic=True)

Encode binary sequences using a generator matrix over GF(2).

Parameters:

• G (np.array of shape (k, n)) – Binary generator matrix.

• bits (np.array of shape (k, N)) – Binary input sequences to encode. Each column is a bit sequence of length \(k\).

• systematic (bool, optional) – If True, the generator matrix is assumed to be in systematic form. If False, the generator matrix is treated as a general linear transformation (default is True).

Returns:

codewords – Encoded binary codewords. Each column corresponds to a codeword of length \(n\) resulting from applying the generator matrix to the respective input bit sequence.

Return type:

np.array of shape (n, N)

gaussElim(M)

Perform Gaussian elimination over GF(2) to reduce a binary matrix to row echelon form.

Parameters:

M (np.array of uint8) – Input binary matrix (dense, 2D). All operations are over GF(2).

Returns:

matrix – Matrix in row echelon form (mod 2).

Return type:

np.array of uint8

hammingParityCheckMatrix(m, extended=False)

Generate the parity-check matrix H for a Hamming code with m check bits .

Parameters:

• m (int) – Number of check bits.

• extended (bool, optional) – If True, generate the parity-check matrix for the extended Hamming code (default is False).

Returns:

H – Parity-check matrix for the Hamming code, where m = n - k and k is the number of message bits.

Return type:

np.array of shape (m, n)

Notes

• For standard Hamming codes, the code parameters are \(n = 2^m - 1\) and \(k = n - m\).

• For extended Hamming codes, the code parameters are \(n = 2^m\) and \(k = n - m - 1\).

References

[1] R. W. Hamming, “Error detecting and error correcting codes,” Bell System Technical Journal, vol. 29, no. 2, pp. 147-160, April 1950.

[2] S. Lin and D. J. Costello, “Error Control Coding: Fundamentals and Applications,” 2nd Edition, Pearson, 2004.

inverseMatrixGF2(A)

Invert a square binary matrix over GF(2) using Gauss-Jordan elimination.

Parameters:

A (np.array of shape (n, n), dtype=np.uint8) – Binary square matrix.

Returns:

• Ainv (np.array of shape (n, n), dtype=np.uint8) – Inverse of A over GF(2), if invertible. If not invertible, returns the identity matrix.

• success (bool) – True if A is invertible, False otherwise.

minSumAlgorithm(llrs, checkNodes, varNodes, maxIter, prec=<class 'numpy.float32'>)

Performs belief propagation decoding using the Min-Sum Algorithm (MSA) for multiple codewords. Optimized to use sparse edge-based message passing and O(d_c) check node updates.

Parameters:

• llrs (np.array of shape (n, numCodewords)) – Log-likelihood ratios (LLRs) of the received codeword bits.

• checkNodes (list of np.array) – List of length \(m\), where each entry contains the indices of variable nodes connected to the corresponding check node.

• varNodes (list of np.array) – List of length \(n\), where each entry contains the indices of check nodes connected to the corresponding variable node.

• maxIter (int) – Maximum number of iterations for belief propagation.

• prec (data-type, optional) – Numerical precision to use in computations (default is np.float32).

Returns:

• finalLLR (np.array of shape (n,)) – Updated LLR values for the decoded codeword after the final iteration.

• numIter (int) – Number of iterations performed before successful decoding or reaching maxIter.

• frameDecodingFail (np.array of shape (numCodewords,)) – Array indicating whether decoding was successful (0) or failed (1) for each codeword. A value of 0 indicates successful decoding, while 1 indicates failure.

References

[1] M. P. C. Fossorier, M. Mihaljevic and H. Imai, “Reduced complexity iterative decoding of low-density parity check codes based on belief propagation,” IEEE Transactions on Communications, vol. 47, no. 5, pp. 673-680, May 1999

par2gen(H)

Convert a binary parity-check matrix H into a systematic generator matrix G over GF(2).

Parameters:

H (np.array of shape (n - k, n)) – Parity-check matrix with entries in {0, 1}.

Returns:

• G (np.array of shape (k, n)) – Systematic generator matrix corresponding to the input parity-check matrix H. The form of G is \([I_k | P]\), where \(I_k\) is the identity matrix and \(P\) is a binary matrix.

• colSwaps (np.array of shape (n,)) – Indices representing the column permutations applied to H to obtain Hm. These permutations are required to match the systematic form used in G.

• Hm (np.array of shape (n - k, n)) – The modified parity-check matrix after Gaussian elimination and column reordering. This matrix has the identity portion on the right-hand side, corresponding to a standard systematic form \([P^T | I_{n-k}]\).

Notes

This method assumes that the parity-check matrix H is full rank and that the identity portion can be isolated on the right via column permutations.

parseAlist(path)

Parse an LDPC ALIST file and return basic parameters.

Parameters:

path (str) – Path to the folder with the ALIST files.

plotBinaryMatrix(H)

Plot the binary matrix H with dots at positions where H[i,j] = 1.

Parameters:

H (np.array of shape (m, n)) – Binary matrix.

readAlist(filename)

Read an ALIST file and reconstruct the binary parity-check matrix H.

Parameters:

filename (str) – Path to the ALIST file.

Returns:

H – Reconstructed binary parity-check matrix.

Return type:

np.array of shape (m, n)

sumProductAlgorithm(llrs, checkNodes, varNodes, maxIter, prec=<class 'numpy.float32'>)

Performs belief propagation decoding using the sum-product algorithm (SPA) for multiple codewords.

Parameters:

• llrs (np.array of shape (n, numCodewords)) – Array of log-likelihood ratios (LLRs) for each bit of the received codeword.

• checkNodes (list of np.array) – List of length \(m\), where each element is a 1D array containing the indices of variable nodes (bits) involved in the corresponding check node (parity-check equation).

• varNodes (list of np.array) – List of length \(n\), where each element is a 1D array containing the indices of check nodes that the corresponding variable node participates in.

• maxIter (int) – Maximum number of belief propagation iterations.

• prec (data-type, optional) – Data type for the computations (default is np.float32).

Returns:

• finalLLR (np.array of shape (n,)) – Updated log-likelihood ratios after message passing.

• numIter (int) – Number of iterations executed until decoding converged or reached maxIter.

• frameDecodingFail (np.array of shape (numCodewords,)) – Array indicating whether decoding was successful (0) or failed (1) for each codeword. A value of 0 indicates successful decoding, while 1 indicates failure.

References

[1] F. R. Kschischang, B. J. Frey and H. . -A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 498-519, Feb 2001.

[2] T. J. Richardson and R. L. Urbanke, “The capacity of low-density parity-check codes under message-passing decoding,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 599-618, Feb 2001.

summarizeAlistFolder(folderPath)

Scan a folder for .alist files and print summary table.

Parameters:

folderPath (str) – Path to the folder containing ALIST files.

triangP1P2(H)

Convert a binary parity-check matrix H into a lower-triangular form and extract matrices P1 and P2.

Parameters:

H (np.array of shape (m, n)) – Binary parity-check matrix. It is used to derive the matrices P1 and P2.

Returns:

• P1 (np.array of shape (m1, k)) – First parity matrix.

• P2 (np.array of shape (m2, k)) – Second parity matrix.

• triangH (np.array of shape (m, n)) – Triangularized H matrix.

References

[1] T. J. Richardson and R. L. Urbanke, “Efficient encoding of low-density parity-check codes,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 638-656, Feb 2001.

triangularize(H)

Convert binary matrix H into lower-triangular form using only row and column permutations.

Parameters:

H (np.array of shape (m, n), dtype=np.uint8) – Binary parity-check matrix \(H\).

Returns:

• triangH (np.array of shape (m, n), dtype=np.uint8) – Triangularized matrix.

• rowPerm (np.array of shape (m,), dtype=np.int32) – Row permutation indices.

• colPerm (np.array of shape (n,), dtype=np.int32) – Column permutation indices.

References

[1] T. J. Richardson and R. L. Urbanke, “Efficient encoding of low-density parity-check codes,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 638-656, Feb 2001.

writeAlist(H, filename)

Save a binary parity-check matrix H (numpy array) to ALIST format.

Parameters:

• H (np.array of shape (m, n)) – Binary parity-check matrix.

• filename (str) – Name of the ALIST file to be written.