function [x_intf,x_sigma]= New_Turbo_IC(x_rake,Former_Rou,Former_Mask,Latter_Rou,Latter_Mask, Noise_variance,x_mean,x_variance,SubslotData_length,Path_number,Det_intl_table,nIter)


% Gray coded 16-QAM mapping vector
MapMatrix16QAM = [-1-j -1-3*j -1+j -1+3*j -3-j -3-3*j -3+j -3+3*j 1-j 1-3*j 1+j 1+3*j 3-j 3-3*j 3+j 3+3*j];
MapMatrix16QAM = MapMatrix16QAM/sqrt(10);

[M,N]=size(Former_Mask); L=N/M;
Former_Mask2=abs(Former_Mask.^2);
Latter_Mask2=abs(Latter_Mask.^2);

variance1 = sum(vec2mat(sum(Former_Mask2),M))+Former_Rou*Noise_variance;
[Rou_sorted,Order1]=sort((Former_Rou.^2)./(variance1+0.001));
Order1 = Order1(end:-1:1);

variance2 = sum(vec2mat(sum(Latter_Mask2),M))+Latter_Rou*Noise_variance;
[Rou_sorted,Order2]=sort((Latter_Rou.^2)./(variance2+0.001));
Order2 = Order2(end:-1:1);


x_mean0=x_mean;
x_variance0=x_variance;

for t=1:nIter
    
    %	Former Half Subslot
    for mm=1:M
        m = Order1(mm);
        Mask_m  = Former_Mask(:,m:M:end);
        Mask2_m = Former_Mask2(:,m:M:end);

        for n=1:SubslotData_length/2 
            x_intf(m,n) = x_rake(m,n)-sum(sum(x_mean0(:,n+(0:L-1)).*Mask_m));

            x_sigma(m,n) = abs(sum(sum(x_variance0(:,n+(0:L-1)).*Mask2_m))+Former_Rou(m)*Noise_variance+0.001);
            
            E=abs(x_intf(m,n)-Former_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n);
            
            P =exp(-E);      
            
            if max(P)<0.005
                x_mean0(m,Path_number-1+n)=0;
                x_variance0(m,Path_number-1+n)=1;      
            else
                
                [Y,I] = sort(P);
                I_reseved=I(11:16);
                P(I_reseved)=P(I_reseved)/sum(P(I_reseved));
                x_mean0(m,Path_number-1+n)=P(I_reseved)*MapMatrix16QAM(I_reseved).';
                x_variance0(m,Path_number-1+n)=abs(P(I_reseved)*abs((MapMatrix16QAM(I_reseved)).^2)'-abs(x_mean0(m,Path_number-1+n))^2);
            end
%             P =exp(-abs(x_intf(m,n)-Former_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n))/abs(sum(exp(-abs(x_intf(m,n)-Former_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n)))+0.001);
% 
%             x_mean0(m,Path_number-1+n)=P*MapMatrix16QAM.';
%             x_variance0(m,Path_number-1+n)=P*abs(MapMatrix16QAM.^2)'-abs(x_mean0(m,Path_number-1+n))^2;   
        end
    end
    
    %	Latter Half Subslot
    for mm=1:M
        m = Order2(mm);
        Mask_m  = Latter_Mask(:,m:M:end);
        Mask2_m = Latter_Mask2(:,m:M:end);

        for n=1+SubslotData_length/2:SubslotData_length 
            x_intf(m,n) = x_rake(m,n)-sum(sum(x_mean0(:,n+(0:L-1)).*Mask_m));

            x_sigma(m,n) = abs(sum(sum(x_variance0(:,n+(0:L-1)).*Mask2_m))+Latter_Rou(m)*Noise_variance+0.001);
            
            E=abs(x_intf(m,n)-Latter_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n);
            
            P =exp(-E);      
            
            if max(P)<0.005
                x_mean0(m,Path_number-1+n)=0;
                x_variance0(m,Path_number-1+n)=1;      
            else
                
                [Y,I] = sort(P);
                I_reseved=I(11:16);
                P(I_reseved)=P(I_reseved)/sum(P(I_reseved));
                x_mean0(m,Path_number-1+n)=P(I_reseved)*MapMatrix16QAM(I_reseved).';
                x_variance0(m,Path_number-1+n)=abs(P(I_reseved)*abs((MapMatrix16QAM(I_reseved)).^2)'-abs(x_mean0(m,Path_number-1+n))^2);
            end

%             P =exp(-abs(x_intf(m,n)-Latter_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n))/abs(sum(exp(-abs(x_intf(m,n)-Latter_Rou(m)*MapMatrix16QAM).^2/x_sigma(m,n)))+0.001);
% 
%             x_mean0(m,Path_number-1+n)=P*MapMatrix16QAM.';
%             x_variance0(m,Path_number-1+n)=P*abs(MapMatrix16QAM.^2)'-abs(x_mean0(m,Path_number-1+n))^2;   
        end
    end
        
end