Compensators#

class BlindDualMIMOCompensator(L: int = 10, *, alphabet: ndarray, mu: float = 0.0001, oversampling: int = 1, mode: Literal['cma', 'rde', 'dd'] = 'cma', sub_block_length: int = 20, name: str = 'mimo filter')[source]#

Blind adaptive 2x2 MIMO equalizer (CMA, RDE or DD stochastic gradient).

Signal Model#

At each step \(n\), the two output polarizations are computed from the estimated equalizer matrix \(\mathbf{W}\) of size \(2 \times 2(2L+1)\) – written \(\mathbf{W}\), not \(\mathbf{H}\), because everywhere else in this package \(\mathbf{H}\) is the channel; the attribute is H_ for historical reasons, the trailing underscore marking it as estimated (D23):

\[y_i[n] = \mathbf{w}_i^H \tilde{\mathbf{x}}[n], \qquad \tilde{\mathbf{x}}[n] = \begin{bmatrix} x_0[n] & \cdots & x_0[n-2L] & x_1[n] & \cdots & x_1[n-2L] \end{bmatrix}^T\]

where \(\mathbf{w}_i^T\) is the \(i\)-th row of \(\mathbf{W}\) and \(x_0\), \(x_1\) are the two received polarizations. The equalizer is updated by stochastic gradient descent, \(\mathbf{W} \leftarrow \mathbf{W} + \mu \mathbf{g}[n]\), minimizing one of the following blind losses:

Constant Modulus Algorithm (CMA): \(\mathcal{L}_{CMA}(y[n]) = \left| R - |y[n]|^2 \right|^2\), where \(R = \mathbb{E}[|s|^4] / \mathbb{E}[|s|^2]\) is derived from the alphabet;
Radius Directed Equalization (RDE): \(\mathcal{L}_{RDE}(y[n]) = \left| \mathcal{P}_{rad}^2(|y[n]|) - |y[n]|^2 \right|^2\), where \(\mathcal{P}_{rad}\) projects onto the list of alphabet radii;
Decision Directed (DD): \(\mathcal{L}_{DD}(y[n]) = \left| \mathcal{P}_{\mathcal{M}}(y[n]) - y[n] \right|^2\), where \(\mathcal{P}_{\mathcal{M}}\) projects onto the alphabet.

The equalizer is adaptive (decision D22): every forward call updates \(\mathbf{H}\) sample by sample, and partial_fit(X) runs one adaptation pass over X from the current state.

Staged adaptation#

CMA is the only one of the three that converges from a cold start. RDE and DD both need decisions that are already roughly right: started from the identity equalizer on 16-QAM, RDE stalls – 3.5 dB where the noise floor is 24. The usual recipe is to run CMA first and switch once it has opened the eye, which is what a ['cma', 'rde'] schedule means in the literature. Two ways to do it here:

between passes, with partial_fit() and set_params() (or a plain assignment) to change mode and mu;
inside a single pass, by overriding process_after_iteration(), which is called after every output sample with the most recent outputs and exists for exactly this (see examples/mimo/one_shot_CMA.py).

Measured on a static Jones rotation with a 24 dB noise floor, the three stages give 21.0, 23.9 and 23.98 dB: CMA opens the eye, RDE takes it to within a tenth of the floor, and DD closes what is left.

Axes: polarization pair – expects (..., 2, N), produces (..., 2, N // oversampling) (one output sample per oversampling input samples: fractionally spaced equalizer). Leading axes are batch, and the butterfly is a per-pair algorithm: each pair gets its own equalizer state, adapted independently. After a batched pass H_ carries the batch axes in front, (..., 2, 2(2L+1)), so a further partial_fit resumes each pair from its own state.

Delay. The centre tap of the initial equalizer sits \(L\) input samples back, so output sample \(k\) is aligned with input sample \(k \cdot os - L\): the block has a group delay of \(L\) input samples, and comparing its output with the transmitted symbols without accounting for it reads as pure noise. The first \((2L+1)/os\) output samples are never written – the filter has no history yet – and are left at zero.

Cost, and why the loop is a loop. While it adapts, this block costs about 18 us per output sample: \(\mathbf{H}\) is updated from \(y[n]\) before \(y[n+1]\) is computed, so the samples cannot be processed together – that is what a stochastic gradient is, and no amount of vectorization removes it. The time goes into a dozen tiny numpy calls per iteration, each dominated by Python overhead rather than arithmetic, so the only real way down is a compiled kernel.

There is one exact exception, and it is the useful one. With \(\mu = 0\) the equalizer is frozen, every output is the same linear map, and the pass becomes a sequence of matrix products – measured 25 to 30 times faster on 100 000 output samples, and bit-for-bit the same output. That is the “fit, then apply” regime of decision D22: adapt on a preamble with partial_fit(), set mu = 0, and run the payload through. The fast path is skipped when process_after_iteration() is overridden, since that hook is defined per output sample.

param L:: Half-length of the filter (each row of \(\mathbf{H}\) holds \(2L+1\) taps per polarization). Default is 10.
type L:: int, optional
param alphabet:: Modulation alphabet, used to derive the CMA radius \(R\), the RDE radius list and the DD decisions.
type alphabet:: np.ndarray, keyword-only
param mu:: Step size \(\mu\) of the gradient update. Default is 1e-4.
type mu:: float, optional, keyword-only
param oversampling:: Oversampling factor of the input. When greater than one, the algorithm implements a fractionally spaced equalizer. Default is 1.
type oversampling:: int, optional, keyword-only
param mode:: Blind loss used for the update. Default is "cma". See Staged adaptation above: only "cma" converges from a cold start.
type mode:: {“cma”, “rde”, “dd”}, optional, keyword-only
param sub_block_length:: Number of recent output samples handed to process_after_iteration(). Default is 20.
type sub_block_length:: int, optional, keyword-only
param name:: Name of the processor instance. Default is "mimo filter".
type name:: str, optional, keyword-only

H_#

Estimated equalizer matrix of shape (2, 2*(2*L+1)) (decision D23: the trailing underscore distinguishes it from the configured channel H of FlatMIMOChannel). Initialized to the identity equalizer (center tap) and updated at each call.

Type:: np.ndarray

raises ShapeError:: If the input does not carry a polarization pair (..., 2, N), or if H_ carries the batch shape of a previous pass and the new input carries a different one – call initialize_H() to restart from the identity equalizer.

References

D. N. Godard, “Self-recovering equalization and carrier tracking in two-dimensional data communication systems,” IEEE Transactions on Communications, vol. 28, no. 11, pp. 1867-1875, 1980.
M. S. Faruk, S. J. Savory, “Digital signal processing for coherent transceivers employing multilevel formats,” Journal of Lightwave Technology, vol. 35, no. 5, pp. 1125-1141, 2017.

Examples

>>> alphabet = np.array([1+1j, 1-1j, -1+1j, -1-1j]) / np.sqrt(2)
>>> compensator = BlindDualMIMOCompensator(L=2, alphabet=alphabet)
>>> rng = np.random.default_rng(0)
>>> X = alphabet[rng.integers(0, 4, size=(2, 200))]
>>> Y = compensator(X)
>>> print(Y.shape, compensator.H_.shape)
(2, 200) (2, 10)

process_after_iteration(n: int, Y_sub: ndarray) → None[source]#

Hook called after each output sample; does nothing by default.

Override it to schedule the stages of the adaptation (switch mode once the eye is open) or to correct a residual phase.

Parameters:

n (int) – Index of the output sample just produced.
Y_sub (np.ndarray) – The last sub_block_length outputs, (2, K), oldest first.

partial_fit(X: ndarray)[source]#: One adaptation pass over X (adaptive regime of decision D22): the estimate H_ keeps evolving from its current state.