1
$\begingroup$

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;
Rayleigh_fading_scaling_factor = sqrt(0.5);
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');

Please check my results below enter image description here

I also put everything together as Dan suggested in the comments enter image description here

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

$\endgroup$
11
  • 2
    $\begingroup$ 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. $\endgroup$ Mar 20, 2022 at 10:58
  • $\begingroup$ @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. $\endgroup$
    – Pawel
    Mar 21, 2022 at 13:27
  • $\begingroup$ Thanks- can you add what you get prior to equalization and what you get in both cases without fading (just AWGN)? $\endgroup$ Mar 21, 2022 at 13:33
  • $\begingroup$ 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? $\endgroup$
    – Pawel
    Mar 21, 2022 at 13:55
  • 1
    $\begingroup$ 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. $\endgroup$
    – Pawel
    Mar 21, 2022 at 19:35

1 Answer 1

1
$\begingroup$

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)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.