function y1 = spheredecod3(rev, C, M, codebook)

% Sphere decoding as the channel matrix is available
% only for real number matrices and vectors
% rec: received signal vector
% C:  Initial radius
% M: channel transfer function matrix (real number matrix)
% codebook (be real number)

[m,n]=size(M);
T = zeros(1,n);

[Q1,Q] = qr(M);
Q = Q(1:n, :);
ro = Q1'*rev;
for ii = 1:n
    if Q(ii,ii) == 0
        Q(ii,ii) = 0.001;
    end
    for jj = ii+1:n
        Q(ii,jj) = Q(ii,jj)/Q(ii,ii);
    end
    ro(ii) = ro(ii)/Q(ii,ii);
    Q(ii,ii) = Q(ii,ii)^2;
end

d = C;
T(n) = C;
S = ro;
cn = length(codebook);
len = zeros(1,n);
Y = zeros(n, cn);
ind = zeros(1,n);
y1 = zeros(n,1);
yt = [];

while 1
    k = n;
    while k >= 1
        U = sqrt(T(k)/Q(k,k))+S(k);
        ii = cn;
        while ii >= 1 && U < codebook(ii)-10^(-11)
            ii = ii - 1;
        end
        
        L = -sqrt(T(k)/Q(k,k))+S(k);
        jj = 1;
        while jj <= cn && L > codebook(jj)+10^(-11)
            jj = jj + 1;
        end
    
        len(k) = ii-jj+1;
        
        Y(k,:) = zeros(1, cn);
        yt = [];
        for kk = 1:len(k)
            yt = [yt codebook(jj+kk-1)];
        end
        [yt1, idex] = sort(abs(yt - S(k)));
        yt = yt(idex);
        for kk = 1:len(k)
            Y(k,kk) = yt(kk);
        end
        
        ind(k)=0;
        d1 = d+1;
        while d1 > d - 10^(-11)
            ind(k) = ind(k) + 1;
            while ind(k) > len(k)
                if k == n
                    return
                end
                k = k+1;
                ind(k) = ind(k) + 1;
            end
            
            if k > 1
                T(k-1) = T(k) - Q(k,k)*(S(k)-Y(k,ind(k)))^2;
                tmp = 0;
                for ii = k:n
                    tmp = tmp + Q(k-1, ii)*Y(ii,ind(ii)) ;
                end
                S(k-1) = ro(k-1) - tmp;
                k = k - 1;
                break
            end
            d1 = T(n) - T(1) + Q(1,1)*(S(1) - Y(1, ind(1)))^2;
        end %while d1 > d
        if d1 < d-10^(-11)
            for l=1:n
                y1(l) = Y(l, ind(l));
            end
            d = d1;
            T(n)=d1;
            break;
        end
    end %while k>=1
    
end %while 1 k=n