Engee documentation
Notebook

2x2 MIMO spatial multiplexing: comparison of ZF and MMSE detectors

The example demonstrates the operation of a Spatial Multiplexing system in a 2x2 MIMO channel with Rayleigh fading. Two independent QPSK streams are transmitted over two antennas, and linear stream separation is applied on the receiving side. Two classical linear detectors, Zero—Forcing (ZF) and Minimum Mean Square Error (MMSE), are compared according to the MSE and BER criteria, with constellation visualization.

The reception model looks like:

where — vector of transmitted characters, — channel matrix (), — complex white Gaussian noise.

Zero-Forcing (ZF) completely suppresses inter-stream distortion by reversing the channel:

Disadvantage: with a small singular number the noise is amplified by a factor of , which dramatically worsens the MSE.

MMSE minimizes the RMS error considering the noise dispersion :

The regularizer prevents explosive growth of noise in the direction of weak eigenvalues at the cost of a small bias. That is why, with low and moderate SNRs and poor channel conditioning, MMSE significantly outperforms ZF.

Realization

The script consists of the following steps:

  1. Data generation — a QPSK character matrix (2 streams × N characters) and a random channel matrix are generated . The channel is artificially conditioned (through multiplication by fixed matrices) in order to demonstrate the advantage of MMSE.
  2. Channel and noise — the signal passes through , after which AWGN is added to it using the system object EngeeComms.AWGN (parameters EbNo, BitsPerSymbol, SignalPower).
  3. Detection — matrices are calculated and through stable LU decomposition (operator \), and are applied to the received signal.
  4. Quality assessment — the MSE is calculated for both detectors and the conditioning number is output .
  5. Visualization — 6 constellations are constructed: the received signal on two antennas (a mixture of streams), the results of ZF and MMSE for each stream.
  6. BER-swipe — in the Monte Carlo SNR cycle, BER is evaluated, and error counting is performed through the system object EngeeComms.ErrorRateCalculation.

Using Engee system Objects (EngeeComms.AWGN, EngeeComms.ErrorRateCalculation) ensures the correct calculation of energy parameters and error statistics in accordance with library standards, eliminating routine operations and making the code suitable for transfer to full-fledged communication models.

In [ ]:
# Pkg.add(["LinearAlgebra", "Random", "Statistics"])
using LinearAlgebra, Random, Statistics, EngeeComms

First, let's perform the initialization of the system:

  • Random.seed!(2024) — fixing the RNG for reproducibility.
  • Parameters: N=4000 characters/stream, SNR=8 дБ, signal strength 1. Recalculation for the AWGN block and for MMSE.
  • QPSK: the constellation is normalized to a single power, a matrix of symbols is formed .
  • Channel : the Rayleigh matrix it is artificially conditioned by multiplication by fixed matrices in order to — this will clearly demonstrate the advantage of MMSE over ZF.
  • y_ideal = H * s — the signal passes through the channel (without noise): The streams are mixed, and the task of the detectors is to separate them.
In [ ]:
Random.seed!(2024)

N = 4000
SNR_dB = 8.0
bits_per_symbol = 2
SignalPower = 1.0

EbNo_dB = SNR_dB - 10 * log10(bits_per_symbol)
σ² = 1.0 / (10^(SNR_dB / 10))

qpsk = ComplexF64[1+1im, -1+1im, -1-1im, 1-1im] ./ sqrt(2)
s = qpsk[rand(1:4, 2, N)]

H = (randn(2,2) .+ 1im .* randn(2,2)) ./ sqrt(2)
H = H * [1.0 0.3; 0.0 0.15] * [1.0 0.0; 0.4 1.0]

y_ideal = H * s
Out[0]:
2×4000 Matrix{ComplexF64}:
 -0.286017+1.34245im   0.833343-1.40624im   …  -1.88977-0.222226im
  0.179459-0.60078im  -0.512083+0.458318im     0.933405+0.321921im

Next, add the noise:

  • EngeeComms.AWGN(...) — creation of a system object that simulates additive white Gaussian noise. Parameters: EbNo (ratio of bit energy to noise), BitsPerSymbol=2 (for QPSK), SignalPower=1.0 (signal strength).
  • awgn_channel(vec(y_ideal)) — passing the signal through the channel with noise. Function vec() transforms the matrix into a column vector, since the Engee API accepts one-dimensional arrays.
  • reshape(y_vec, size(y_ideal)) — restoring the original shape of the matrix for further processing by detectors.
In [ ]:
awgn_channel = EngeeComms.AWGN(
    EbNo = EbNo_dB,
    BitsPerSymbol = bits_per_symbol,
    SignalPower = SignalPower
)

y_vec = awgn_channel(vec(y_ideal))
y = reshape(y_vec, size(y_ideal))
Out[0]:
2×4000 Matrix{ComplexF64}:
 -0.233481+1.91278im    0.304577-1.46275im  …  -1.90713-0.316705im
  0.363028-0.816211im  -0.842276+1.09409im     0.553998+0.858639im

The next stage is detection and quality assessment. What's going on here:

  1. Calculation of detector arrays W_zf and W_mmse according to the pseudo-circulation and Wiener-Hopf formulas.
  2. Applying detectors to the received signal y — separation of mixed flows.
  3. The MSE calculation is a quantitative assessment of the accuracy of character reconstruction for both methods.
In [ ]:
W_zf = (H' * H) \ H'
ŝ_zf = W_zf * y

W_mmse = (H' * H + σ² * I) \ H'
ŝ_mmse = W_mmse * y

mse_zf = mean(abs2, ŝ_zf - s)
mse_mmse = mean(abs2, ŝ_mmse - s)

println("MIMO 2×2: ZF vs MMSE")
println("SNR (Es/N0) = $SNR_dB дБ  |  Eb/No = $(round(EbNo_dB, digits=2)) дБ")
println("cond(H)     = $(round(cond(H), digits=2))")
println("MSE  ZF     = $(round(mse_zf, digits=5))")
println("MSE  MMSE   = $(round(mse_mmse, digits=5))")
MIMO 2×2: ZF vs MMSE
SNR (Es/N0) = 8.0 дБ  |  Eb/No = 4.99 дБ
cond(H)     = 23.17
MSE  ZF     = 11.04402
MSE  MMSE   = 0.50019

After that, the constellations are visualized.

  • draw(x; ...) — auxiliary function that builds a constellation in two layers:
  • translucent cloud of accepted points (I/Q);
  • red crosses — ideal QPSK positions for comparison.
  • 6 graphs: the received signal on two antennas (gray, a mixture of streams), the results of ZF (blue) and MMSE (green) for each of the two streams.
In [ ]:
function draw(x; title="", color=:blue)
    p = scatter(real.(x), imag.(x);
                markersize=3, alpha=0.35, color=color, title=title,
                xlabel="I", ylabel="Q",
                xlims=(-2.5, 2.5), ylims=(-2.5, 2.5),
                aspect_ratio=:equal, legend=:topright, grid=true)
    scatter!(p, real.(qpsk), imag.(qpsk);
             markersize=8, color=:red, marker=:xcross)
    return p
end

p_y1 = draw(y[1,:]; title="Rx antenna 1", color=:gray)
p_y2 = draw(y[2,:]; title="Rx antenna 2", color=:gray)
p_z1 = draw(ŝ_zf[1,:]; title="Stream 1 after ZF", color=:blue)
p_z2 = draw(ŝ_zf[2,:]; title="Stream 2 after ZF", color=:blue)
p_m1 = draw(ŝ_mmse[1,:]; title="Stream 1 after MMSE", color=:green)
p_m2 = draw(ŝ_mmse[2,:]; title="Stream 2 after MMSE", color=:green)

plot(p_y1, p_y2, p_z1, p_z2, p_m1, p_m2;
     layout=(3,2), size=(900, 1100),
     plot_title="MIMO 2×2: ZF vs MMSE (SNR = $SNR_dB дБ)",
     left_margin=8Plots.mm, bottom_margin=8Plots.mm)
Out[0]:

Chart output:

Top row — 4 clusters are visible on the antennas (a mixture of streams), 2 are more blurred on the antenna.

The average range (ZF) is a catastrophic spread:

  • Stream 1: a huge blue cloud, dots fill the entire field
  • Stream 2: even worse — the dots are evenly smeared, the QPSK structure is completely lost

The bottom row (MMSE) is dramatically better:

  • Stream 1: 4 clear green clusters, perfect recovery
  • Stream 2: one cluster merged in the center.

Initialization ErrorRateCalculation

  • snrs_dB = 0:2:20 — the range of signal-to-noise ratios (from 0 to 20 dB in increments of 2) for constructing BER curves.
  • ber_zf, ber_mmse — empty arrays for accumulating BER results for each SNR for the ZF and MMSE detectors.
  • Nsym = 20_000 — the number of characters in one run (iteration) to ensure statistical reliability.
  • n_iterations = 50 is the number of Monte Carlo iterations (of various random channel implementations) for each SNR point.
  • err_calc = EngeeComms.ErrorRateCalculation() — creation of the Engee system object, which will be used to calculate the number of errors in the subsequent steps of the cycle.
In [ ]:
snrs_dB = 0:2:20
ber_zf, ber_mmse = Float64[], Float64[]
Nsym = 20_000
n_iterations = 50

err_calc = EngeeComms.ErrorRateCalculation()
Out[0]:
ErrorRateCalculation:
    ReceiveDelay=0
    ComputationDelay=0
    ComputationMode=Entire frame
    OutputData=Port
    ResetPort=false
    StopSimulation=false

There are two nested loops in the code below. The external iterates through the SNR values from 0 to 20 dB in increments of 2, initializes error counters, and creates an object. EngeeComms.AWGN with the recalculated value for the current step . The internal Monte Carlo cycle performs 50 independent implementations: each generates a new random Rayleigh channel matrix. and a block of transmitted characters , the signal passes through a noisy channel, after which the ZF and MMSE detectors are used. The hard solution assigns each recovered character to the nearest QPSK point, and the system object ErrorRateCalculation calculates the number of errors that accumulate over all iterations.

At the end of the inner loop, SER is calculated as the ratio of the total number of errors to the total number of characters transmitted, and then converted to BER, taking into account the fact that for QPSK with Gray encoding . The resulting values are added to arrays ber_zf and ber_mmse.

The final graph is plotted on a logarithmic scale along the ordinate axis with two curves — ZF and MMSE.

In [ ]:
for snr_db in snrs_dB
    err_zf_total = 0
    err_mmse_total = 0
    
    ebno_db_sweep = snr_db - 10 * log10(bits_per_symbol)
    awgn_sweep = EngeeComms.AWGN(EbNo = ebno_db_sweep, BitsPerSymbol = bits_per_symbol, SignalPower = SignalPower)
    
    for _ in 1:n_iterations
        H_ = (randn(2,2) .+ 1im.*randn(2,2)) ./ sqrt(2)
        s_ = qpsk[rand(1:4, 2, Nsym)]
        
        y_ideal_sweep = H_ * s_
        y_sweep_vec = awgn_sweep(vec(y_ideal_sweep))
        y_sweep = reshape(y_sweep_vec, size(y_ideal_sweep))

        σ2_sweep = 1.0 / (10^(snr_db / 10))
        ŝz = ((H_'*H_) \ H_') * y_sweep
        ŝm = ((H_'*H_ + σ2_sweep*I) \ H_') * y_sweep

        decide(x) = sign(real(x))/sqrt(2) + 1im*sign(imag(x))/sqrt(2)
        bits_zf = decide.(ŝz)
        bits_mmse = decide.(ŝm)
        
        stats_zf = err_calc(vec(bits_zf), vec(s_))
        stats_mmse = err_calc(vec(bits_mmse), vec(s_))
        
        err_zf_total += stats_zf[2]
        err_mmse_total += stats_mmse[2]
    end
    
    ser_zf = err_zf_total / (n_iterations * Nsym)
    ser_mmse = err_mmse_total / (n_iterations * Nsym)
    
    push!(ber_zf, ser_zf / 2.0)
    push!(ber_mmse, ser_mmse / 2.0)
end

p_ber = plot(snrs_dB, ber_zf; yscale=:log10, marker=:circle, lw=2, label="ZF",
             xlabel="SNR (Es/N0), дБ", ylabel="BER", title="BER: ZF vs MMSE (2×2 MIMO)",
             grid=true, legend=:bottomleft)
plot!(p_ber, snrs_dB, ber_mmse; marker=:square, lw=2, label="MMSE", color=:green)
Out[0]:

Why the graphs matched. Unlike the demonstration of constellations, where the channel was artificially conditioned (), BER uses a "pure" Rayleigh channel. For a random matrix with complex Gaussian elements, the probability of getting a poorly conditioned implementation is low — on average not much. With good conditioning and moderate SNR, the regularizer In MMSE, it becomes negligible, and both detectors are actually reduced to pseudo-circulation. . Therefore, the ZF and MMSE curves practically merge on the graph: the theoretical gain of MMSE here is "smeared" over a variety of favorable channel implementations and is not visible against the background of statistical averaging. To identify the difference on the BER curve, it would be necessary either to fix a poorly conditioned channel (as in the first part of the example), or to shift the SNR range to very low values.

Conclusion

The example clearly demonstrates the fundamental compromise between two classical linear detectors in MIMO systems.

Zero-Forcing completely eliminates inter-thread interference, but at the cost of catastrophic noise amplification with poor channel conditioning. In our experiment, when and the SNR=8 dB MSE of the ZF detector was 11.04 — the constellation turned into a uniform cloud, and the QPSK structure was completely lost.

MMSE sacrifices full MSI suppression for the sake of noise control. The regularizer prevents explosive growth of dispersion in the directions of weak channel eigenvalues. The result: MSE = 0.50, which is 22 times better than ZF. In constellations, this appears as clear compact clusters around ideal QPSK positions.

Key insight: the advantage of MMSE over ZF is maximized with low SNRs and poor channel conditioning — precisely in those conditions where noise dominates the signal. With high SNR or good conditioning both detectors converge to a pseudo-circulation , and their performance is almost identical (as confirmed by the coincidence of BER curves in the range of 0-20 dB for a random Rayleigh channel).

Practical recommendation: in real communication systems, where the channel may temporarily become poorly conditioned (fading, antenna correlation), MMSE is preferable to ZF as a more robust detector. The computational complexity of both methods is comparable (both require matrix inversion ), therefore, the choice in favor of MMSE does not incur additional costs.