How to implement the sphere decoding algorithm for MIMO systems?

Question

I am fairly new to wireless communications, and I am trying to simulate an M-QAM MIMO system. I want to decode the received signal using QR decomposition and sphere decoding using Matlab.

However, I am struggling to understand how to implement the simple tree sphere decoding algorithm. Here's my understanding so far:

Given a receive signal as follows:

1. Y = HX + N;

where X is a K x 1 vector consisting of K MQAM data symbols, H is a K x K fading matrix, N is a K x 1 AWGN vector, and Y is the K x 1 received vector.

We perform QR detection as follows:

2. H = QR;

where Q is a K x K Unitary matrix, and R is a upper-triangular matrix with entry R(i, j), with i and j being the row and column indices respectively.

Then we equalize the receive signal in 1. as:

3. Z = (Q^H)*Y = RX + (Q^H)*N ,

where Q^H is the Hermitian transpose of Q.

I now want to create a function that performs sphere decoding on these K MQAM symbols based on 3. This is where I am struggling.

I have tried to implement for 16QAM symbols, using the sphere decoding algorithm for real operations.

To do so, the matrices Y, H, and X are written in terms of their real and imaginary parts as:

Y = [real(Y); imag(Y)]

H = [real(H), -imag(H); imag(H), real(H)]

X = [real(X); imag(X)]

The set of integers for 16QAM symbols is [-3, -1, 1, 3].

I then tried to design a function for the sphere decoding algorithm for real operations given an arbitrary radius d, and the 16QAM set as follows:

function output = Real_Sphere_Decoder(r, z, d)

n = length(z);
s = zeros(n, 1);
step = 2; % MQAM integer sets always spaced out by 2

%% Vectors
d_n = zeros(1, n); % radii at each level
UB = zeros(1, n); % upper bound at each level
LB = zeros(1, n); % lower bound at each level
interSum = zeros(1, n); % partial metric at each level
output = zeros(n, 1); % output variable

%% Initialise Parameters
d_n(n) = d;
UB(n) = floor((d_n(n) + z(n))/r(n, n)); % initialise upper bound
LB(n) = ceil((-d_n(n) + z(n))/r(n, n)); % initialise lower bound
minR = d_n(n); % initialise minimum euclidean distance
interSum(n) = z(n);

%% Initial 'if' statement - UB must be greater than LB
if LB(n) > UB(n)
    error('initial radius too small');
else
    s(n) = LB(n);
end

%% While loop
k = n;
while (s(n) <= UB(n))
    if s(k) > UB(k)
        k = k + 1;
        s(k) = s(k) + step;
    else % s(k) <= UB(k)
        if k > 1
            k = k - 1;
            d_n(k) = sqrt(d_n(k + 1)^2 - (interSum(k + 1) - r(k + 1, k + 1)*s(k + 1))^2);
            interSum(k) = z(k) - r(k, (k + 1):n)*s((k + 1):n);
            UB(k) = floor((d_n(k) + interSum(k))/r(k, k));
            LB(k) = ceil((-d_n(k) + interSum(k))/r(k, k));
            if LB(k)> UB(k)
                k = k + 1;
                s(k) = s(k) + step;
            else
                s(k) = LB(k);
            end
        else % k == 1
            while (s(k) <= UB(k))
                if minR > norm(z - r*s, 'fro')
                    for ii = 1:n
                    output(ii) = s(ii);
                    end
                    minR = norm(z - r*s, 'fro');
                end
                s(k) = s(k) + step;
            end
            k = k + 1;
            s(k) = s(k) + step;
        end
    end
end
    
end

The function takes in R, Z, and d, and outputs the sphere decoded signal given by the 'output' variable.

I get one of 3 errors/problems:

1. the LB is always greater than the UB when R(n, n) is negative, hence the algorithm always terminates.

2. the output variable 'output' is often empty.

3. I do not know how to round the upper and lower bounds towards the nearest integers given by the 16QAM integer set. I used the floor and ceil function, but I am not sure if that is how you do it.