Sphere Decoding with QR decomposition

How can I embedd shere decoding with QR decomposition in the code below to detect golden codes

clear
clc
format shorte
M=16                                                     %Modulation Order or Number of Symbols
m=log2(M)                                                %Number of Bits per Symbol
QAM_16=qammod([0:M-1],M)                                 %QAM16 Modulation or QAM16 Signal Set
Epsilon=(sum(abs(QAM_16).^2))/length(QAM_16)             %Average Sum of Energy or Energy per Symbol
Nt=2
Nr=4
Ntr=Nt*Nr
SNR_db=[0:2:20]                                          %SNR in dB
SNR_dec=(10.^(SNR_db./10))/Nt                            %SNR in Decimal or Linear SNR

Num_Err= 500*ones(1,length(SNR_db))         %Termination of Each SNR

Teta=(1+sqrt(5))/2;
Teta_=1-Teta;
Alpha=1+(j*Teta_);
Alpha_=1+j*(1-Teta_);
Gamma=j;

for kk=1:length(SNR_db)
    total_error=0;                                       %Total Errors Initialization
    t_error=0;
    num_sim=0;                                           %Initialization of Number of Simulations for Each SNR 
    while(total_error<Num_Err(kk))        
        bin1= round(rand(1,m));                          %m Bits Random Binary Input 1
        bin2= round(rand(1,m));                          %m Bits Random Binary Input 2
        bin3= round(rand(1,m));                          %m Bits Random Binary Input 3
        bin4= round(rand(1,m));                          %m Bits Random Binary Input 4
        
        
        bin1_dec1=1;
        for k=1:m
            bin1_dec1=bin1_dec1+(2^(m-k))*bin1(k);       %Binary conversion to Decimal 1
        end
        
        bin2_dec2=1;
        for k=1:m
            bin2_dec2=bin2_dec2+(2^(m-k))*bin2(k);       %Binary conversion to Decimal 2
        end
        
        bin3_dec3=1;
        for k=1:m
            bin3_dec3=bin3_dec3+(2^(m-k))*bin3(k);       %Binary conversion to Decimal 1
        end
        
        bin4_dec4=1;
        for k=1:m
            bin4_dec4=bin4_dec4+(2^(m-k))*bin4(k);       %Binary conversion to Decimal 2
        end
        
        x1=QAM_16(bin1_dec1);                          %Modulated Signal or Transmitted Signal x1
        x2=QAM_16(bin2_dec2);                          %Modulated Signal or Transmitted Signal x2
        x3=QAM_16(bin3_dec3);                         %Modulated Signal or Transmitted Signal x3
        x4=QAM_16(bin4_dec4);                         %Modulated Signal or Transmitted Signal x4
        
        x11=(1/sqrt(5))*Alpha*(x1+x2*Teta);
        x22=(1/sqrt(5))*Alpha_*(x1+x2*Teta_);
        x12=(1/sqrt(5))*Alpha*(x3+x4*Teta);
        x21=(1/sqrt(5))*Gamma*Alpha_*(x3+x4*Teta_);
        
        X1=[x11;x12];
        X2=[x21;x22];
        X=[X1,X2];
        
        n1_I=sqrt(Epsilon/(2*SNR_dec(kk)))*randn(Nr,1);  
        n1_Q=sqrt(Epsilon/(2*SNR_dec(kk)))*randn(Nr,1);  
        n1=n1_I+1j*n1_Q;
        
        n2_I=sqrt(Epsilon/(2*SNR_dec(kk)))*randn(Nr,1);  
        n2_Q=sqrt(Epsilon/(2*SNR_dec(kk)))*randn(Nr,1);  
        n2=n2_I+1j*n2_Q;
        
        h1_I=sqrt(1/2)*randn(Nr,1);                      
        h1_Q=sqrt(1/2)*randn(Nr,1);                      
        h1=h1_I+1j*h1_Q;
        
        h2_I=sqrt(1/2)*randn(Nr,1);                      
        h2_Q=sqrt(1/2)*randn(Nr,1);                      
        h2=h2_I+1j*h2_Q;
        
        y1=h1*x11+h2*x12+n1;                             %Received Signal y1
        y2=h2*x21+h1*x22+n2;                             %Received Signal y2
        
        y=[y1 y2];
        n=[n1 n2];
        
        H=[h1 h2];
        [Q,R]=qr(H);
        z=Q'*y;
        
        n_hat=Q'*n;
        n1_hat=[n_hat(1,1);n_hat(2,1);n_hat(3,1);n_hat(4,1)];
        n2_hat=[n_hat(1,2);n_hat(2,2);n_hat(3,2);n_hat(4,2)];
        
        z1=R(1,1)*x11+R(1,2)*x12+n1_hat;
        z2=R(2,1)*x21+R(2,2)*x22+n2_hat;
        
        z11= R(1,1)*x11+R(2,1)*x12+n1_hat;
        z12= R(1,2)*x11+R(2,2)*x12+n1_hat;
        z21= R(3,1)*x11+R(3,2)*x12+n2_hat;
        z22= R(4,1)*x21+R(4,2)*x22+n2_hat;
        
        
        ED_1=zeros(1,M^4);
        Store_1=zeros(M^4,4);
        Count_1=0;
        for jj1=1:M                                       %Maximum Likelihood Detection Loop
            for jj2=1:M
                for jj3=1:M
                    for jj4=1:M
                        Count_1=Count_1+1;
                        Store_1(Count_1,:)=[jj1 jj2 jj3 jj4];
                        
                        x_1=QAM_16(jj1);
                        x_2=QAM_16(jj2);                        
                        x_3=QAM_16(jj3);
                        x_4=QAM_16(jj4);                        
                        
                        x_11=(1/sqrt(5))*Alpha*(x_1+x_2*Teta);
                        x_22=(1/sqrt(5))*Alpha_*(x_1+x_2*Teta_);
                        x_12=(1/sqrt(5))*Alpha*(x_3+x_4*Teta);
                        x_21=(1/sqrt(5))*Gamma*Alpha_*(x_3+x_4*Teta_);
                        
                        X_1=[x_11;x_12];
                        X_2=[x_21;x_22];
                        
                        ED_1(Count_1) =norm(y1-(h1*x_11+h2*x_12),'fro').^2+norm(y2-(h2*x_21+h1*x_22),'fro').^2;    %Euclidean Distance for each of the 16 Symbols of 16QAM
                        %                         ED_2(Count_1) =norm(y1-(h1*x_11+h2*x_12)).^2+norm(y2-(h2*x_21+h1*x_22)).^2;    %Euclidean Distance for each of the 16 Symbols of 16QAM
                    end
                end
            end
        end                                                                          
        [Ed1,idx1]=min(ED_1);                            %Minimum Euclidean and Index of any of the 16 Symbols of 16QAM Transmitted
        %        Bin_1_Hat=Store_1(idx1,:);                      %Detected Signal X_Hat
        
        x_1_Hat=QAM_16(Store_1(idx1,1));
        x_2_Hat=QAM_16(Store_1(idx1,2));
        
        x_3_Hat=QAM_16(Store_1(idx1,3));
        x_4_Hat=QAM_16(Store_1(idx1,4));
        
        if x1~=x_1_Hat                                                       
            total_error= total_error+ 1;                     %Total Error Per Simulation
        end
        
        if x2~=x_2_Hat                                                       
            total_error= total_error+ 1;                     %Total Error Per Simulation
        end
        % % %         
        % % %         if x3~=x_3_Hat                                                       
        % % %         total_error= total_error+ 1;                     %Total Error Per Simulation
        % % %         end
        % % %         
        % % %         if x4~=x_4_Hat                                                       
        % % %         total_error= total_error+ 1;                     %Total Error Per Simulation
        % % %         end
        
        bin_1_dec1=Store_1(idx1,1);
        bin_2_dec2=Store_1(idx1,2);
        bin_3_dec3=Store_1(idx1,3);
        bin_4_dec4=Store_1(idx1,4);
        
        bin_1_hat=de2bi((bin_1_dec1-1),m,'left-msb');
        bin_2_hat=de2bi((bin_2_dec2-1),m,'left-msb');
        bin_3_hat=de2bi((bin_3_dec3-1),m,'left-msb');
        bin_4_hat=de2bi((bin_4_dec4-1),m,'left-msb');
        
        
        error=size(find(bin1-bin_1_hat),2);                                                      
        t_error= t_error+ error;                     %Total Error Per Simulation
        
        error=size(find(bin2-bin_2_hat),2);                                                      
        t_error= t_error+ error;                     %Total Error Per Simulation
        
        % % %         error=size(find(bin3-bin_3_hat),2);                                                      
        % % %         t_error= t_error+ error;                     %Total Error Per Simulation
        % % %                 
        % % %         error=size(find(bin4-bin_4_hat),2);                                                      
        % % %         t_error= t_error+ error;                     %Total Error Per Simulation
        
        num_sim=num_sim+1;
    end
    SER_s(kk)=total_error/(num_sim*2)                    %Simulated SER 
    BER_s(kk)=t_error/(num_sim*2*m)
end
semilogy(SNR_db,SER_s,'-b*','LineWidth',2)
hold on
semilogy(SNR_db,BER_s,'-k*','LineWidth',2)
grid on
ylabel('ERROR RATES')
xlabel('SNR (dB)')
title('GOLDEN CODE - ERROR RATES VS SNR (16QAM_2X4)')
legend('SER-SIM','BER-SIM')