# MIMO detection using ZF and MMSE

I have a problem related to my uni assignment. The task seems rather simple - plot SER as a function of Es/N0 for 2x2 MIMO and 16-QAM modulation assuming AWGN noise and Rayleigh flat-fading channel. But the results I got are, shall we say, questionable.

First of all, instead of waterfall-like curve I got a straight line. This implies that SER decays exponentially, not linearly, with Es/N0. I can blame this on interference between MIMO layers. The second issue is much serious though - performance of MMSE and ZF is almost identical, even for low Es/N0. Furthermore, I really need high Es/N0 values to reach 10^-4 SER.

I have no idea what I did incorrectly, therefore I would really appreciate your feedback. My code is as follows

close all;
clear all;

% Declare useful constants
N_users = 2;
N_tx = 1;
N_rx = 2;
s_variance = 10;
detection_threshold = 2;

SNR = 15 : 2 : 55;
SNR_linear = 10 .^ (SNR ./ 10);

% SER
M = length(SNR); % number of SNR steps
SER_MMSE = zeros(1, M);
SER_ZF = zeros(1, M);

for idx = 1 : M
a = (2e6 - 1e4) / (M - 1);
b = 1e4 - a;
L = a * idx + b;
fprintf("Processing SNR = %d, number of realisations = %d\n", SNR(idx), L);
for l = 1 : L
s = 2 * randi(4, N_users, N_tx) - 5 + 1.0j * (2 * randi(4, N_users, N_tx) - 5);

n_scaling_factor = sqrt(N_rx * s_variance / (2 * SNR_linear(idx)));
n = n_scaling_factor * complex(randn(N_rx, N_tx), randn(N_rx, N_tx));
H = Rayleigh_fading_scaling_factor .* complex(randn(N_rx, N_users), randn(N_rx, N_users));

y = H * s + n;

y_mmse = (H' * H + (1 / SNR_linear(idx) * eye(N_users, N_users))) \ (H' * y);
s_hat_mmse = sign(real(y_mmse)) .* (1 * (abs(real(y_mmse)) < detection_threshold) + 3 * (abs(real(y_mmse)) > detection_threshold)) + 1.0j * sign(imag(y_mmse)) .* (1 * (abs(imag(y_mmse)) < detection_threshold) + 3 * (abs(imag(y_mmse)) > detection_threshold));
SER_MMSE(idx) = SER_MMSE(idx) + sum(s_hat_mmse ~= s, 'all');

y_zf = (H' * H) \ (H' * y);
s_hat_zf = sign(real(y_zf)) .* (1 * (abs(real(y_zf)) < detection_threshold) + 3 * (abs(real(y_zf)) > detection_threshold)) + 1.0j * sign(imag(y_zf)) .* (1 * (abs(imag(y_zf)) < detection_threshold) + 3 * (abs(imag(y_zf)) > detection_threshold));
SER_ZF(idx) = SER_ZF(idx) + sum(s_hat_zf ~= s, 'all');
end
SER_MMSE(idx) = SER_MMSE(idx) / (L * N_tx * N_users);
SER_ZF(idx) = SER_ZF(idx) / (L * N_tx * N_users);
end

figure;
semilogy(SNR, SER_MMSE);
hold on;
semilogy(SNR, SER_ZF);
grid on;
legend('MMSE', 'ZF', 'Location', 'NorthEast');
xlabel('E_s / N_0 [dB]');
ylabel('SER');
title('SER vs SNR');


I also put everything together as Dan suggested in the comments

If there is no fading and no equalization, or I use identity matrix as H, results I get match what I found in the literature. But if there is either fading and no equalization or no fading and but equlization, then SER is close to 100%. When I put fading and equalization, then results are "in between". But this still does not explain why ZF and MMSE gives almost the same results.

Kind regards

• If you could show your plots that may help illustrate your question and get an answer. It's less likely someone will have the time to go through your whole code line by line. Mar 20, 2022 at 10:58
• @DanBoschen Thank you for your suggestion, I updated my post. Also, I suppose that actual ZF / MMSE is correct. But what I am not sure of is if I generate noise correctly. Mar 21, 2022 at 13:27
• Thanks- can you add what you get prior to equalization and what you get in both cases without fading (just AWGN)? Mar 21, 2022 at 13:33
• I am not sure what do you want me to do. If I replace y = H * s + n with y = s + n, but still do the equalization, then I get really poor results. Also, I don't get what do you mean by "prior to equalization". Do you mean hard detection results? Mar 21, 2022 at 13:55
• I suppose I know why ZF and MMSE results are similar regardless of Es/N0. Whereas ZF amplifies the noise more significantly than MMSE, this really becomes a problem if H is not full rank. Provided that I do hard detection, even though noise amplification is smaller for MMSE, I am not able to utilize this. Mar 21, 2022 at 19:35

Indeed, ZF and MMSE are not suitable when the attend ratio is low ($$N_{r}/N_{t}$$), some may find detectors like K-Best (Sphere Decoding) more suitable. I would suggest trying M-MIMO, like $$8*128$$, $$16*128$$, and perhaps will obtain a better result.

Moreover, sometimes Matlab gets really strange when doing round-off functions, for me, Python is much better. I will provide a sample code as follows (I generally use MMSE in my research):


import time
import torch
import numpy as np
import multiprocessing
import matplotlib.pyplot as plt
from math import sqrt
from scipy.linalg import toeplitz
# Parameters Setting
TestDataLen = 30000
TxAntNum = 8   # Number of transmitting antennas
RxAntNum = 16   # Number of receiving antennas tested
DataLen  = 1    # Length of data sequence
MaxIter = 4     # Max iteration number of iterative algorithms
delta = 0.5     # Damping factor

# SNR setting
SNRdBLow = 0   # Minimal value of SNR in dB
SNRdBHigh = 12   # Maximal value of SNR in dB
SNRIterval = 2  # Interval value of SNR sequence
SNRNum = int((SNRdBHigh-SNRdBLow)/SNRIterval)+1   # Number of SNR sequence
SNRdB = np.linspace(SNRdBLow, SNRdBHigh, SNRNum)
TSratio = 0
RSratio = 0
Kron = np.sqrt(TSratio*RSratio)
# Constellation Setting
ModType = 8
if ModType==2:
Model = '4QAM'
Cons = torch.tensor([-1., 1.])
fnorm = 1/np.sqrt(2)
bitperSym = 1
normCons = fnorm*Cons
elif ModType==4:
Model = '16QAM'
Cons = torch.tensor([-3., -1., 1., 3.])
bitCons = torch.tensor([[0, 0], [0, 1], [1, 0], [1, 1]])
bitperSym = 2
fnorm = 1/np.sqrt(10)
normCons = fnorm*Cons
elif ModType==8:
Model = '64QAM'
Cons = torch.tensor([-7., -5., -3., -1., 1., 3., 5., 7.])
bitCons = torch.tensor([[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]])
bitperSym = 3
fnorm = 1/np.sqrt(42)
normCons = fnorm*Cons
elif ModType==16:
Model = '256QAM'
Cons = torch.tensor([-15., -13., -11., -9., -7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13., 15.])
bitCons = torch.tensor([[0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 1, 1], [0, 1, 0, 0], [0, 1, 0, 1], [0, 1, 1, 0], [0, 1, 1, 1],
[1, 0, 0, 0], [1, 0, 0, 1], [1, 0, 1, 0], [1, 0, 1, 1], [1, 1, 0, 0], [1, 1, 0, 1], [1, 1, 1, 0], [1, 1, 1, 1]])
bitperSym = 4
fnorm = 1/np.sqrt(170)
normCons = fnorm*Cons

# Variable Initialization
error_MMSE = np.zeros(SNRNum)

# Correlation channel
def KronChannel(TxAntNum, RxAntNum, DataLen, TSratio, RSratio):
# 1. Generate Rayleigh channel matrix
Hiid = np.sqrt(1 / 2) * (torch.randn(DataLen, RxAntNum, TxAntNum) + 1j * torch.randn(DataLen, RxAntNum, TxAntNum))

# 2. Generate the Transmit UpperMatrix
indexT = torch.tensor(toeplitz(-torch.arange(TxAntNum), torch.arange(TxAntNum))).unsqueeze(dim=-3)
randT = torch.rand([DataLen, 1, 1])
phaseT = torch.exp(1j * randT * np.pi / 2 * indexT)
ampT = TSratio ** torch.abs(indexT)
Rt = ampT * phaseT
Ct = torch.linalg.cholesky(Rt, upper=True)

# 3. Generate the Receive UpperMatrix
indexR = torch.tensor(toeplitz(-torch.arange(RxAntNum), torch.arange(RxAntNum))).unsqueeze(dim=-3)
randR = torch.rand([DataLen, 1, 1])
phaseR = torch.exp(1j * randR * np.pi / 2 * indexR)
ampR = RSratio ** torch.abs(indexR)
Rr = ampR * phaseR
Cr = torch.linalg.cholesky(Rr, upper=True)

Hkron = Cr.matmul(Hiid.matmul(Ct.conj().transpose(-1, -2)))

return Hkron

# Generating Data Sent, Recived, Noise, and Transmitting Matrix
def GenerateTestData(SampleNum, TxAntNum, RxAntNum, DataLen, SNRdBLow, SNRdBHigh):
x_ = torch.zeros([SampleNum, 2*TxAntNum, DataLen])
y_ = torch.zeros([SampleNum, 2*RxAntNum, DataLen])
H_ = torch.zeros([SampleNum, 2*RxAntNum, 2*TxAntNum])
Nv_ = torch.zeros([SampleNum, 1, 1])
H = KronChannel(TxAntNum, RxAntNum, SampleNum, TSratio, RSratio)
for itr in range(SampleNum):
RandSNRdB = np.random.uniform(low=SNRdBLow, high=SNRdBHigh)
SNR = 10**(RandSNRdB/10)
# Generate real and image part of data sequence seperately
TxDataSyms = torch.randint(0, ModType, size=(2*TxAntNum, DataLen))
# TxData_r, TxData_i = Modulation(TxDataBits, ModType, Cons)
TxData_r = Cons[TxDataSyms[:TxAntNum, :]]
TxData_i = Cons[TxDataSyms[TxAntNum:, :]]
# Transform complex Tx signals to real
TxData = torch.cat((TxData_r, TxData_i), dim=0)
x_[itr, :, :] = TxData
# TxSymbol = np.concatenate((TxPilot, TxData), 1)
TxSymbol = TxData
# Generate channel matrix (Rayleigh channel)
# Hc = np.sqrt(1/2)*(torch.randn(RxAntNum, TxAntNum) + 1j*torch.randn(RxAntNum, TxAntNum))
Hc = H[itr]
# Transform complex channle matrix to real
HMat = torch.cat((torch.cat((torch.real(Hc), -torch.imag(Hc)), 1),
torch.cat((torch.imag(Hc), torch.real(Hc)), 1)), 0).float()
# Normalize the column of real channel matrix
HMat /= np.sqrt(torch.norm(HMat)**2/(2*TxAntNum))
# Data send via the channels without AWGN
RxSymbol_noAWGN = torch.matmul(HMat, TxSymbol)
# Calculate the norm of channel matrix
Hnorm = torch.norm(HMat)**2
# Noise variance & adding AWGN
Nv = (1*Hnorm)/(2*SNR*(2*RxAntNum)*fnorm**2)
RxSymbol = RxSymbol_noAWGN + np.sqrt(Nv)*torch.randn(2*RxAntNum, DataLen)
Hhat = HMat
y_[itr, :, :] = RxSymbol
H_[itr, :, :] = Hhat
Nv_[itr, :, :] = Nv

return x_, y_, H_, Nv_

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
print(device)

def MMSEtest(x, y, H, Nv):
error_MMSE = 0

HTy = H.transpose(-2, -1).matmul(y)
HTH = H.transpose(-2, -1).matmul(H)

# MIMO detection (MMSE) using perfect CSI
Sigma = torch.inverse(HTH + Nv * fnorm ** 2 * torch.unsqueeze(torch.eye(2 * TxAntNum), 0))
xhat = torch.matmul(Sigma, HTy)

# calculation of BER
_, indices = torch.min((xhat - Cons) ** 2, dim=-1, keepdim=True)
_, indices_x = torch.min((x - Cons) ** 2, dim=-1, keepdim=True)
Rxbit = bitCons[indices]
xbit = bitCons[indices_x]
comp = torch.where(
Rxbit != xbit)  # return tuple (Tensor, Tensor, Tensor), where contains the 3D coordinates of non-zero elements
error_MMSE += len(comp[0])

return error_MMSE

if __name__ == "__main__":
print("cpu count:", multiprocessing.cpu_count(), "\n")
for nEN in range(SNRNum):
print(SNRdB[nEN])

for i in range (5):
# we do this because the process of "GenerateTestData" sometimes uses
# too much space, thus, we must need this loop to maintain the
# number of TestDataLen

x, y, H, Nv = GenerateTestData(TestDataLen, TxAntNum, RxAntNum,
DataLen, SNRdB[nEN], SNRdB[nEN])

error_MMSE[nEN] = MMSEtest(x, y, H, Nv)

ber_MMSE = error_MMSE / (2 * TestDataLen * TxAntNum * bitperSym)

plt.figure(1)
# plt.style.use("_classic_test_patch")
p1 = plt.semilogy(SNRdB, ber_MMSE, 'b-o', label='MMSE')
plt.legend()
plt.grid()
plt.xlabel('SNR')
plt.ylabel('BER')
plt.title(str(RxAntNum) + r'$$\times$$' + str(TxAntNum) + ', MIMO, ' + str(Model) + ', Kron' +str(Kron))

plt.show()


P.S. (I prefer to use BER instead of SER hhh)