Having incorrect results for NOMA 16QAM BER.

I've been practising and exploring NOMA QAM. I've been working on code for the sum BER of two user 16QAM NOMA, but have met with results that don't seem right. The BER straightens out at about 0.3 and this seems like a too high of BER. I've done a BPSK and QPSK version which seem correct with them going to small values but my QAM (16 QAM and up) version has too high of BER.

The code:

clc; clear variables; 

%Parameters
N = 5*10^5;                          %Number of monte carlo simulations/samples

Pt = -114:5:-24;			        %Transmit power (dBm)
pt = (10^-3)*db2pow(Pt);	        %Transmit power (linear scale)

BW = 10^6;			                %Channel bandwidth = 1MHz
No = -174 + 10*log10(BW);	        %Noise power (dBm)		
no = (10^-3)*db2pow(No);	        %Noise power (linear scale) 

SNR = Pt - No;   %Signal to noise ratio (DB) = Pt - No when not linear scale

eta = 4;         %Path loss exponent (for typical propagation environments 2-4)

d1 = 10; d2 = 10;         %Distance of users in meters
a1 = 0.8; a2 = 0.2;       %Power allocation coefficients 

%%% 16-QAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 16;          %Modulation order alphabet
k = log2(M);     %number of bit per symbols

%Rayleigh fading coefficients of both users
h1 = (sqrt(d1^-eta))*(randn(N/k,1) + 1i*randn(N/k,1))/sqrt(2);
h2 = (sqrt(d2^-eta))*(randn(N/k,1) + 1i*randn(N/k,1))/sqrt(2);

%Generate noise samples for the users
n1 = sqrt(no)*(randn(N/k,1) + 1i*randn(N/k,1))/sqrt(2); 
n2 = sqrt(no)*(randn(N/k,1) + 1i*randn(N/k,1))/sqrt(2);

%Generate random binary message data for the users
x1 = randi([0 1],N,1);
x2 = randi([0 1],N,1);

%Perform QAM modulation of data
xmod1 = qammod(x1,M,'UnitAveragePower',true,InputType="bit");
xmod2 = qammod(x2,M,'UnitAveragePower',true,InputType="bit");

%Do super position coding
x = sqrt(a1)*xmod1 + sqrt(a2)*xmod2;

for u = 1:length(Pt)
%Received signals
   y1 = sqrt(pt(u))*x.*h1 + n1;     %At user 1
   y2 = sqrt(pt(u))*x.*h2 + n2;     %At user 2

%Perform equalization
   eq1 = y1./h1;
   eq2 = y2./h2;

%Performing SIC process
%Decode at user 1 (Direct decoding)
   dec1 = qamdemod(eq1,M,'UnitAveragePower',true,OutputType="bit");         %Direct demodulation of signal 
   
%Decode at user 2
   dec12 = qamdemod(eq2,M,'UnitAveragePower',true,OutputType="bit");		%Direct demodulation to get U1's data
   dec12_remod = qammod(dec12,M,'UnitAveragePower',true,InputType="bit");	%Remodulation of U1's data
   rem2 = eq2 - sqrt(a1*pt(u))*dec12_remod;	                                %SIC to remove U1's data
   dec2 = qamdemod(rem2,M,'UnitAveragePower',true,OutputType="bit");		%Direct demodulation of remaining signal


%BER calculation
   ber1(u) = biterr(dec1, x1)/N;
   ber2(u) = biterr(dec2, x2)/N;

   ber_sum16QAM(u) = (ber1(u) + ber2(u))/2;      
end

%{
cd = comm.ConstellationDiagram(ShowReferenceConstellation=false);       %Constellation
cd(eq1)
%}

%Plot
figure()
semilogy(SNR, ber_sum16QAM, '-o', 'linewidth', 2); hold on; grid on;

xlabel('SNR (dB)');
ylabel('BER');

[Shashi Kiran]

Hi Louis,

I came across this IEEE paper which demonstrates the same observation: higher-order modulation schemes like 16-QAM and 64-QAM for NOMA result in worse BER (around 0.3-0.4).

https://ieeexplore.ieee.org/document/8399043

This might be helpful.