Utilities

class Constellation(family: str, order: int, *, labelling: str = 'gray', norm: bool = True)

A modulation alphabet, and everything that is a property of it.

Signal Model

The alphabet get_alphabet() returns, plus the facts that always travel with it. A memoryless modulation is a set \(\mathcal{A} = \{a_0, \ldots, a_{M-1}\}\), and everything a study asks of it – how many bits it carries, what its average energy is, how close its nearest neighbours are, what error rate it reaches at a given SNR, what rate it can carry – is determined by the family and the order \(M\). Passing those two around separately is how a page ends up drawing the closed form of a 16-QAM under the measurement of a 64-QAM, which nothing checks.

The object is a drop-in for the array it wraps: np.asarray returns the alphabet, so every block and every function that takes one keeps working.

param family:

Modulation family: "PAM", "PSK" or "QAM".

type family:

str

param order:

Modulation order \(M\).

type order:

int

param labelling:

Bit-to-symbol mapping, "gray" (default) or "binary".

type labelling:

str, optional, keyword-only

param norm:

Rescale to unit average symbol energy. Default True, which is the convention the rest of the library assumes.

type norm:

bool, optional, keyword-only

alphabet

The symbols \(a_m\), in labelling order.

Type:

np.ndarray

bits_per_symbol

\(k = \log_2 M\).

Type:

int

energy

Average symbol energy \(E_s\), over an equiprobable input.

Type:

float

min_distance

Smallest distance between two symbols, which sets the high-SNR error rate.

Type:

float

raises ValueError:

If the family, the order or the labelling is unknown. The check is get_alphabet()’s, made once here instead of at every call site.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 4.

Examples

>>> qam = Constellation("QAM", 16)
>>> qam.bits_per_symbol, round(qam.energy, 6), round(qam.min_distance, 4)
(4, 1.0, 0.6325)
>>> np.asarray(qam).shape            # a drop-in for the array
(16,)
>>> print(f"{qam.metrics(10.0)['ser']:.4e}")      # at Eb/N0 = 10 dB
7.0043e-03

info() → dict[str, Any]

What the constellation is, as a dictionary ready to print.

Returns:

dict – family, order, labelling, bits_per_symbol, energy, min_distance, and papr_dB – the last being the peak-to-average power ratio of the constellation, a lower bound on what the shaped waveform will show.

Examples

>>> print(Constellation("PSK", 8).info()["papr_dB"])
0.0

plot(ax: Axes | None = None, **kwargs: Any) → Axes

Draw the constellation diagram (decision D25).

Parameters:

• ax (matplotlib.axes.Axes or None, optional) – Axis to draw on. If None, a new figure and axis are created.

• **kwargs – Forwarded to plot_alphabet().

Returns:

matplotlib.axes.Axes – The axis containing the plot.

Examples

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> Constellation("QAM", 16).plot().get_title()
'16-QAM'

metrics(snr_dB: ArrayLike, *, per: str = 'bit', channel: str = 'awgn', diversity: int = 1, metrics: Tuple[str, ...] = ('ser', 'ber'), px: ndarray | None = None) → dict[str, ndarray | float]

What this constellation is worth at a given SNR.

Signal Model

A front-end to the four quantities a page draws against a measurement, all of them functions of the same family, order and SNR:

• "ser" and "ber", the closed-form error rates of compute_metric_awgn_theo() or compute_metric_rayleigh_theo();

• "mi" and "gmi", the mutual information and the bit-interleaved rate of constellation_capacity() and bicm_capacity().

The two families are parameterized differently by convention: the error rates against \(E_b/N_0\), the rates against the SNR per symbol, which is \(k\) times larger. The object knows \(k\), so per= says which of the two is being passed and the conversion happens here rather than at the call site – get it wrong and the curve moves by \(10\log_{10} k\) dB, with nothing to signal it.

"mi" and "gmi" are quadratures over \(M^2\) points, not closed forms, so they are not in the default and cost seconds rather than microseconds on a large constellation.

param snr_dB:

Signal-to-noise ratio in dB – the quantity a chain is swept over.

type snr_dB:

float or np.ndarray

param per:

"bit" (default) if snr_dB is \(E_b/N_0\), "symbol" if it is \(E_s/N_0\).

type per:

str, optional, keyword-only

param channel:

"awgn" (default) or "rayleigh", for the error rates. The rates below are AWGN only.

type channel:

str, optional, keyword-only

param diversity:

Number of combined branches, for channel="rayleigh".

type diversity:

int, optional, keyword-only

param metrics:

Which quantities to return, among "ser", "ber", "mi" and "gmi". Default ("ser", "ber").

type metrics:

tuple of str, optional, keyword-only

param px:

Symbol distribution, for a shaped input. The constellation is rescaled to unit energy under this law before the rates are computed – a shaped input compared as it stands is simply a quieter one. Not available for the error rates, whose closed forms assume an equiprobable input.

type px:

np.ndarray, optional, keyword-only

returns:

dict of str to float or np.ndarray – One entry per requested metric, in the order requested, each with the shape of snr_dB.

raises ValueError:

If per, channel or a metric name is unknown, if a rate is asked of a fading channel, or if px is given together with an error rate.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Section 4.3; G. Ungerboeck, “Channel coding with multilevel/phase signals”, IEEE Trans. Inf. Theory 28(1), 1982 (the rate a constellation carries).

Examples

>>> qpsk = Constellation("PSK", 4)
>>> theory = qpsk.metrics(np.array([4.0, 8.0]))
>>> print(np.array2string(theory["ber"], precision=6))
[0.012423 0.000191]
>>> rate = qpsk.metrics(20.0, per="symbol", metrics=("mi",))
>>> print(round(float(rate["mi"]), 3))
2.0

get_alphabet(modulation: str, order: int, type: str = 'gray', norm: bool = True) → ndarray

Return the symbol alphabet of a memoryless modulation of order \(M\).

Signal Model

A memoryless modulation is described by its alphabet \(\mathcal{A} = \{a_0, \ldots, a_{M-1}\}\): the mapper turns the symbol index \(s[n] \in \{0, \ldots, M-1\}\) into the transmitted sample

\[x[n] = a_{s[n]}\]

The raw constellations are tabulated on the usual grids – odd integers \(\pm 1, \pm 3, \ldots\) for PAM and QAM, the unit circle for PSK – and the row order of the file is the bit-to-symbol mapping (Gray or natural binary).

With norm=True (the default) the alphabet is rescaled to unit average symbol energy, symbols being assumed equiprobable:

\[a_m \leftarrow \frac{a_m}{\sqrt{E_s}}, \qquad E_s = \frac{1}{M}\sum_{m=0}^{M-1} \left|a_m\right|^2\]

so that \(\mathbb{E}\left[|x[n]|^2\right] = 1\). This is the convention the rest of the library assumes: with \(E_s = 1\) an absolute noise variance \(\sigma^2\) is \(N_0\), and the MMSE regularizer \(\sigma^2 \mathbf{I}\) of mmse_estimator() needs no rescaling. With norm=False the raw grid is returned instead (\(E_s = 10\) for 16-QAM, for instance) and those expressions must be rescaled by \(E_s\).

param modulation:

Modulation family: "PAM", "PSK" or "QAM".

type modulation:

str

param order:

Modulation order \(M\), i.e. \(M = 2^k\) symbols carrying \(k = \log_2(M)\) bits each. Tabulated values are 4, 16, 32, 64, 128 and 256.

type order:

int

param type:

Bit-to-symbol mapping: "gray" or "bin" (natural binary). Default is "gray".

type type:

str, optional

param norm:

If True, normalize the alphabet to unit average symbol energy \(E_s = 1\). Default is True.

type norm:

bool, optional

returns:

np.ndarray – Complex alphabet \(\mathcal{A}\) of length \(M\), indexed by the symbol index \(s\).

Notes

The tables are CSV files of the data subdirectory of this module, named <modulation>_<order>_<type>.csv, with one row per symbol (index, real part, imaginary part) sorted by increasing index.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 3 (memoryless modulation methods and Gray mapping).

Examples

>>> alphabet = get_alphabet("QAM", 4)
>>> print(np.round(alphabet, 4))
[-0.7071+0.7071j -0.7071-0.7071j  0.7071+0.7071j  0.7071-0.7071j]
>>> print(round(float(np.mean(np.abs(alphabet)**2)), 6))
1.0
>>> print(np.round(get_alphabet("QAM", 16, norm=False)[:4], 1))
[-3.+3.j -3.+1.j -3.-3.j -3.-1.j]

plot_alphabet(alphabet: ndarray, ax: Axes | None = None, label: str = 'alphabet', title: str = 'Constellation', **kwargs: Any) → Axes

Plot the constellation diagram of a symbol alphabet.

Each symbol \(a_m\) of the alphabet \(\mathcal{A} = \{a_0, \ldots, a_{M-1}\}\) is drawn as a point of the complex plane, its real part on the horizontal axis and its imaginary part on the vertical one. This is a plotting helper, not a processing block: it transforms no signal.

Parameters:

• alphabet (np.ndarray) – Complex-valued symbol alphabet \(\mathcal{A}\) to plot, typically returned by get_alphabet().

• ax (matplotlib.axes.Axes or None, optional) – Axis to draw on. If None, a new figure and axis are created.

• label (str, optional) – Label for the scatter plot. Default is "alphabet".

• title (str, optional) – Title of the plot. Default is "Constellation".

• **kwargs – Additional keyword arguments forwarded to ax.plot.

Returns:

matplotlib.axes.Axes – The axis containing the plot (decision D25).

Examples

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> ax = plot_alphabet(get_alphabet("PSK", 4))
>>> ax.get_xlabel()
'real part'

sym_2_bin(sym: ndarray, width: int = 4) → ndarray

Expand symbol indices into their binary representation, MSB first.

Signal Model

Each symbol index \(s[n] \in \{0, \ldots, M-1\}\) is written in base 2 on \(k\) bits (\(k =\) width, with \(M \leq 2^k\)):

\[s[n] = \sum_{i=0}^{k-1} b[nk + i] \, 2^{\,k-1-i}, \qquad b[\cdot] \in \{0, 1\}\]

The bits of a symbol are emitted most-significant first and the symbols are concatenated, so a sequence of \(N\) indices produces \(kN\) bits. This is the same ordering as the bit table of SymbolDemapper, which is what makes the BER of comnumpy.core.metrics comparable across blocks.

Axes: element-wise – each index is expanded independently; the input is read as a flat sequence and the output is 1-D of length \(kN\).

param sym:

Sequence of \(N\) integer symbol indices \(s[n]\).

type sym:

array-like

param width:

Number of bits per symbol \(k\). Default is 4.

type width:

int, optional

returns:

np.ndarray – 1-D integer array of \(kN\) bits (0s and 1s).

Notes

Each symbol is represented by exactly width bits: a shorter binary representation is left-padded with zeros. An index larger than \(2^k - 1\) is not rejected – numpy.binary_repr then emits more than width bits and the output loses its block structure.

Examples

>>> print(sym_2_bin([0, 5, 9], width=4))
[0 0 0 0 0 1 0 1 1 0 0 1]
>>> print(sym_2_bin([0, 3], width=2))
[0 0 1 1]

hard_projector(z: ndarray, alphabet: ndarray) → Tuple[ndarray, ndarray]

Project symbols onto the nearest constellation point (hard decision).

Signal Model

Given an equalized observation \(z[n]\) and the alphabet \(\mathcal{A} = \{a_0, \ldots, a_{M-1}\}\), the decision is the minimum-distance (nearest-neighbour) rule

\[\hat{s}[n] = \arg\min_{m \in \{0, \ldots, M-1\}} \left| z[n] - a_m \right|^2, \qquad \hat{x}[n] = a_{\hat{s}[n]}\]

For the AWGN model \(z[n] = x[n] + b[n]\) with \(b[n] \sim \mathcal{CN}(0, \sigma^2)\) and equiprobable symbols, this rule is the maximum-likelihood (hence MAP) detector: the likelihood decreases monotonically with the Euclidean distance, so \(\sigma^2\) does not enter the decision.

Axes: element-wise – the search is carried out over the alphabet on an axis appended internally, so both outputs keep the shape of z. Ties are broken by numpy.argmin (smallest index wins).

param z:

Observation \(z[n]\), of any shape.

type z:

np.ndarray

param alphabet:

1-D constellation alphabet \(\mathcal{A}\) of length \(M\).

type alphabet:

np.ndarray

returns:

• s (np.ndarray) – Integer indices \(\hat{s}[n]\) of the nearest constellation points, with the shape of z.

• x (np.ndarray) – Nearest constellation symbols \(\hat{x}[n] = a_{\hat{s}[n]}\), with the shape of z.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 4 (optimum receivers for the AWGN channel).

Examples

>>> alphabet = get_alphabet("QAM", 4)
>>> z = np.array([0.9+1.1j, -0.2-0.8j])
>>> s, x = hard_projector(z, alphabet)
>>> print(s)
[2 1]
>>> print(np.round(x, 4))
[ 0.7071+0.7071j -0.7071-0.7071j]

soft_projector(z: ndarray, alphabet: ndarray, sigma2: float, kernel: ndarray | None = None) → ndarray

Compute the soft (MMSE) symbol estimate under a Gaussian noise model.

Signal Model

For the scalar observation \(z[n] = x[n] + b[n]\) with \(b[n] \sim \mathcal{CN}(0, \sigma^2)\) and equiprobable symbols drawn from \(\mathcal{A} = \{a_0, \ldots, a_{M-1}\}\), the MMSE estimate is the posterior mean

\[\hat{x}[n] = \mathbb{E}\left[x[n] \mid z[n]\right] = \frac{\displaystyle\sum_{m=0}^{M-1} a_m \, e^{-|z[n]-a_m|^2/\sigma^2}} {\displaystyle\sum_{m=0}^{M-1} e^{-|z[n]-a_m|^2/\sigma^2}}\]

The exponent carries \(\sigma^2\) and not \(2\sigma^2\) because the noise is circularly symmetric complex: its density is \((\pi\sigma^2)^{-1} e^{-|b|^2/\sigma^2}\), \(\sigma^2\) being the total variance of the two quadratures.

Passing a kernel \(g_{m}[n]\) replaces \(a_m\) in the numerator, which turns the same weights into any posterior average

\[\frac{\sum_m g_{m}[n] \, e^{-|z[n]-a_m|^2/\sigma^2}} {\sum_m e^{-|z[n]-a_m|^2/\sigma^2}} = \mathbb{E}\left[g \mid z[n]\right]\]

The AMP detectors of comnumpy.mimo.detectors use \(g_{m}[n] = |a_m - \hat{x}[n]|^2\) to obtain the posterior variance that drives their state evolution.

Axes: element-wise – each observation is processed independently; the input is flattened, so the output is 1-D of size z.size.

param z:

Observation \(z[n]\), flattened internally.

type z:

np.ndarray

param alphabet:

1-D constellation alphabet \(\mathcal{A}\) of length \(M\).

type alphabet:

np.ndarray

param sigma2:

Noise variance \(\sigma^2\) (total, both quadratures). Only its real part is used. Large values flatten the weights towards the alphabet mean, small values recover the hard decision of hard_projector().

type sigma2:

float

param kernel:

Values \(g_{m}[n]\) averaged in place of \(a_m\), of shape (z.size, M) (or broadcastable to it). If None, the alphabet itself is used and the posterior mean is returned.

type kernel:

np.ndarray or None, optional

returns:

np.ndarray – Soft estimates \(\hat{x}[n]\), 1-D of size z.size.

References

• S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, Prentice Hall, 1993 (the MMSE estimator is the posterior mean).

• C. Jeon, R. Ghods, A. Maleki and C. Studer, “Optimality of large MIMO detection via approximate message passing,” Proc. IEEE Int. Symp. Information Theory (ISIT), 2015, pp. 1227-1231 (the F function and its variance kernel).

Examples

>>> alphabet = get_alphabet("QAM", 4)
>>> z = np.array([0.9+1.1j, -0.2-0.8j])
>>> print(np.round(soft_projector(z, alphabet, 0.1), 4))
[ 0.7071+0.7071j -0.7022-0.7071j]
>>> print(np.round(soft_projector(z, alphabet, 100.0), 4))
[ 0.009+0.011j -0.002-0.008j]

esn0_to_snr_dB(esn0_dB: ArrayLike, oversampling: int = 1) → NDArray[float64]

Convert a symbol-energy-to-noise ratio \(E_s/N_0\) (dB) into the SNR (dB) expected by AWGN.

Signal Model

Three quantities must not be confused. \(E_s\) is the energy of one constellation symbol, \(E_b\) the energy of one information bit and \(N_0\) the noise power spectral density; the SNR of AWGN is a per-sample power ratio

\[\mathrm{SNR} = \frac{P_x}{\sigma^2}, \qquad P_x = \mathbb{E}\left[|x[n]|^2\right]\]

At one sample per symbol \(P_x = E_s\) and \(\sigma^2 = N_0\), so the two quantities coincide. Oversampling by \(L\) samples per symbol multiplies the noise bandwidth – hence the per-sample noise variance \(\sigma^2 = L N_0\) – while an energy-preserving pulse shaping leaves the measured signal power at \(E_s\), so

\[\mathrm{SNR_{dB}} = E_s/N_0\big|_{dB} - 10\log_{10}(L)\]

Axes: element-wise – a pure scalar conversion, applied pointwise to an array of operating points when sweeping.

param esn0_dB:

Symbol-energy-to-noise-spectral-density ratio \(E_s/N_0\big|_{dB}\), in dB.

type esn0_dB:

float or np.ndarray

param oversampling:

Number of samples per symbol \(L\). Default is 1.

type oversampling:

int, optional

returns:

float – The SNR in dB, relative to the measured signal power.

See also

ebn0_to_snr_dB

same conversion from the energy per information bit.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 4 (energy-per-symbol and energy-per-bit normalizations).

Examples

>>> print(float(esn0_to_snr_dB(10)))
10.0
>>> print(round(float(esn0_to_snr_dB(10, oversampling=4)), 4))
3.9794

ebn0_to_snr_dB(ebn0_dB: ArrayLike, bits_per_symbol: int, code_rate: float = 1.0, oversampling: int = 1) → NDArray[float64]

Convert a bit-energy-to-noise ratio \(E_b/N_0\) (dB) into the SNR (dB) expected by AWGN.

Signal Model

One constellation symbol carries \(k = \log_2(M)\) channel bits, of which \(kR\) are information bits when a code of rate \(R\) is used. The energy budget of a symbol is therefore shared by \(kR\) information bits:

\[E_s = k R \, E_b\]

Combined with the per-sample SNR of AWGN (see esn0_to_snr_dB(), where \(\sigma^2 = L N_0\) for \(L\) samples per symbol), this gives

\[\mathrm{SNR_{dB}} = E_b/N_0\big|_{dB} + 10\log_{10}(k R) - 10\log_{10}(L)\]

The conversion requires chain-level knowledge (\(k\), \(R\), \(L\)), which is why it lives here and not inside the AWGN block (decision D41). Comparing modulations at equal \(E_b/N_0\) is the fair comparison – at equal \(E_s/N_0\), a dense constellation is unduly favoured because it spends the same energy on more bits.

Axes: element-wise – a pure scalar conversion, applied pointwise to an array of operating points when sweeping.

param ebn0_dB:

Information-bit-energy-to-noise-spectral-density ratio \(E_b/N_0\big|_{dB}\), in dB.

type ebn0_dB:

float or np.ndarray

param bits_per_symbol:

Number of bits carried by one constellation symbol \(k = \log_2(M)\) (e.g. 4 for 16-QAM).

type bits_per_symbol:

int

param code_rate:

FEC code rate \(R \in (0, 1]\). Default is 1.0 (uncoded, \(E_s = k E_b\)).

type code_rate:

float, optional

param oversampling:

Number of samples per symbol \(L\). Default is 1.

type oversampling:

int, optional

returns:

float – The SNR in dB, relative to the measured signal power.

See also

esn0_to_snr_dB

same conversion from the energy per symbol.

References

J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 4 (energy-per-symbol and energy-per-bit normalizations).

Examples

>>> print(round(float(ebn0_to_snr_dB(10, bits_per_symbol=4)), 4))
16.0206
>>> print(round(float(ebn0_to_snr_dB(10, bits_per_symbol=2, code_rate=0.5)), 4))
10.0

zf_estimator(Y: ndarray, H: ndarray) → ndarray

Perform Zero Forcing (ZF) linear equalization using the channel matrix pseudoinverse.

Signal Model

For the flat MIMO observation model

\[\mathbf{y}[n] = \mathbf{H}\mathbf{x}[n] + \mathbf{b}[n], \qquad \mathbf{b}[n] \sim \mathcal{CN}\left(\mathbf{0}, \sigma^2 \mathbf{I}_{N_r}\right)\]

with \(\mathbf{H}\) of size \(N_r \times N_t\), the ZF equalizer applies the Moore-Penrose pseudoinverse:

\[\mathbf{z}[n] = \mathbf{H}^{\dagger}\mathbf{y}[n] = \left(\mathbf{H}^H\mathbf{H}\right)^{-1}\mathbf{H}^H \mathbf{y}[n] = \mathbf{x}[n] + \mathbf{H}^{\dagger}\mathbf{b}[n]\]

the closed form on the right holding when \(\mathbf{H}\) has full column rank (\(N_r \geq N_t\)). Interference between streams is cancelled exactly, at the price of noise enhancement: the residual noise has covariance \(\sigma^2\left(\mathbf{H}^H\mathbf{H}\right)^{-1}\), which blows up when \(\mathbf{H}\) is ill-conditioned. mmse_estimator() trades a residual bias for a bounded noise gain.

Axes: declared axis – expects the MIMO layout (..., ant, N): \(N_r\) receive antennas on axis -2 of Y, time on axis -1, with H a 2-D (N_r, N_t) matrix. The output carries \(N_t\) streams on axis -2.

param Y:

Received samples \(\mathbf{y}[n]\) stacked column-wise, of shape (..., N_r, N).

type Y:

np.ndarray

param H:

Channel matrix \(\mathbf{H}\) of shape (N_r, N_t).

type H:

np.ndarray

returns:

np.ndarray – Equalized signal \(\mathbf{z}[n]\) of shape (..., N_t, N), still to be mapped onto the alphabet by hard_projector().

References

• D. Tse and P. Viswanath, Fundamentals of Wireless Communication, Cambridge University Press, 2005, Chapter 8.

• J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 9 (linear equalization).

Examples

>>> H = np.array([[1.0, 0.5], [0.0, 2.0]])
>>> X = np.array([[1.0, -1.0], [1.0, 1.0]])
>>> Y = np.matmul(H, X)
>>> print(np.round(zf_estimator(Y, H), 4))
[[ 1. -1.]
 [ 1.  1.]]

mmse_estimator(Y: ndarray, H: ndarray, sigma2: float) → ndarray

Perform Minimum Mean Square Error (MMSE) linear equalization.

Signal Model

For the flat MIMO observation model

\[\mathbf{y}[n] = \mathbf{H}\mathbf{x}[n] + \mathbf{b}[n], \qquad \mathbf{b}[n] \sim \mathcal{CN}\left(\mathbf{0}, \sigma^2 \mathbf{I}_{N_r}\right)\]

the linear estimator minimizing \(\mathbb{E}\left[\|\mathbf{z}[n]-\mathbf{x}[n]\|^2\right]\) is

\[\mathbf{z}[n] = \left(\mathbf{H}^H\mathbf{H} + \sigma^2 \mathbf{I}_{N_t}\right)^{-1} \mathbf{H}^H\mathbf{y}[n]\]

The regularization term is the noise-to-signal ratio \((\sigma^2/E_s)\,\mathbf{I}_{N_t}\); it reduces to \(\sigma^2 \mathbf{I}_{N_t}\), as implemented, because the library normalizes the constellation to unit symbol energy \(E_s = 1\) (see get_alphabet()) and the streams are assumed uncorrelated, \(\mathbb{E}[\mathbf{x}[n]\mathbf{x}^H[n]] = \mathbf{I}_{N_t}\). Two limits are worth remembering: \(\sigma^2 \to 0\) gives back the zero-forcing solution of zf_estimator(), while a large \(\sigma^2\) shrinks the estimate towards zero rather than amplifying the noise. The MMSE output is biased, which is why it is followed by a decision on the alphabet, not used as an unbiased estimate.

Axes: declared axis – expects the MIMO layout (..., ant, N): \(N_r\) receive antennas on axis -2 of Y, time on axis -1, with H a 2-D (N_r, N_t) matrix. The output carries \(N_t\) streams on axis -2.

param Y:

Received samples \(\mathbf{y}[n]\) stacked column-wise, of shape (..., N_r, N).

type Y:

np.ndarray

param H:

Channel matrix \(\mathbf{H}\) of shape (N_r, N_t).

type H:

np.ndarray

param sigma2:

Noise variance \(\sigma^2\) per receive antenna, in the same power units as the unit-energy constellation.

type sigma2:

float

returns:

np.ndarray – Equalized signal \(\mathbf{z}[n]\) of shape (..., N_t, N), still to be mapped onto the alphabet by hard_projector().

References

• D. Tse and P. Viswanath, Fundamentals of Wireless Communication, Cambridge University Press, 2005, Chapter 8.

• J. G. Proakis, M. Salehi, Digital Communications, 5th ed., McGraw-Hill, 2008, Chapter 9 (MMSE linear equalization).

Examples

>>> H = np.array([[1.0, 0.5], [0.0, 2.0]])
>>> X = np.array([[1.0, -1.0], [1.0, 1.0]])
>>> Y = np.matmul(H, X)
>>> print(np.round(mmse_estimator(Y, H, 0.1), 4))
[[ 0.9151 -0.8931]
 [ 0.9868  0.9647]]