function [Ser,Ber] = Bler_calculate(B,BS_antennas,User_antennas,P_max,K,P,R,N_ris,sigma2,H_bkp,F_rkp,G_brp,W,Theta)

num = 20000;
snr_dB = -5:20;
Ser = zeros(1,length(snr_dB));
Ber = zeros(1,length(snr_dB));
snr = 10.^(snr_dB/10);
%Rece_Ampli = 10^(-3)* 25.17935662^(- 3);
SNR = snr/(10^(-7));

a = sqrt(0.1);
QAM_16 = [a+1i*a; a+3i*a; 3*a+1i*a; 3*a+3*1i*a;
          a-1i*a; a-3i*a; 3*a-1i*a; 3*a-3*1i*a;
          -a+1i*a; -a+3i*a; -3*a+1i*a; -3*a+3i*a;
          -a-1i*a; -a-3i*a; -3*a-1i*a; -3*a-3i*a];
      
BIT_16 = [0,0,0,0 ; 0,0,0,1 ; 0,0,1,0 ; 0,0,1,1 ;...
          0,1,0,0 ; 0,1,0,1 ; 0,1,1,0 ; 0,1,1,1 ;...
          1,0,0,0 ; 1,0,0,1 ; 1,0,1,0 ; 1,0,1,1 ;...
          1,1,0,0 ; 1,1,0,1 ; 1,1,1,0 ; 1,1,1,1];
H_bkp=reshape(H_bkp,[BS_antennas,User_antennas]);
F_rkp=reshape(F_rkp,[N_ris,User_antennas]);
G_brp=reshape(G_brp,[N_ris,BS_antennas]);

onehot = eye(16);
s = repmat(QAM_16,[num,1]);
onehot_s = repmat(onehot,[num,1]);
randIndex = randperm(size(s,1));
s = s(randIndex,:);
onehot_s = onehot_s(randIndex,:);

H_replace = H_bkp' + F_rkp' * Theta' * G_brp;
% [~,m] = size(H_replace);
% 
% [U,S,V] = svd(H_replace');
% %size(U);
% W = U(:,1);
% W2 = V(:,1);

x = W*s.';

%r = H_replace * W;
%pr = pinv(r);

for i_snr = 1:length(SNR)
    SNR(i_snr);
    noise = sqrt(1/(2*SNR(i_snr)))*(randn(User_antennas,16*num)+1i*randn(User_antennas,16*num));
    
    y = H_replace * x + noise; 
    
    %% SIC
    %yi = y;
    %H_sic = H_replace;
    %x_sic = zeros(m,num*16);
    %for i = 1:m
    %    G = pinv(H_sic);
    %    nG = sum(abs(G).^2,2);
    %    [~, index]=sort(nG);
    %    xj = G(index(i),:)*yi;
    %    xi = zeros(m, num*16);
    %    xi(index(i),:) = xj;
    %    x_sic(index(i),:) = xj;
    %    yi = yi - H_sic * xi;
    %    H_sic(:,index(i)) = 0;
    %end
    %x_hat = W'*x_sic;
    %% MMSE
    %Weight_Matrix = (H_replace'*H_replace + eye(8)/SNR(i_snr))\H_replace'; % ZF¾ùºâ¾ØÕó
    %x_hat = W' * Weight_Matrix * y;
    
    %% ZF
    Weight_Matrix = pinv(H_replace);%(H_replace'*H_replace)\H_replace'; % ZF¾ùºâ¾ØÕó
    x_hat = W' * Weight_Matrix * y;
    
    %% SVD
    %x_hat = W2' * y/norm(r);
        
    %mean(sum(abs(x).^2,1))  
    %mean(sum(abs(G_brp*x).^2,1))
    %mean(sum(abs(Theta * G_brp * x).^2,1))
    %mean(sum(abs(F_rkp' * Theta' * G_brp * x).^2,1))
    %mean(sum(abs(H_bkp' * x).^2,1))
    %mean(sum(abs(y).^2,1))
    
    b = [-2 0 2]*a; 
    c = [-3 -1 1 3]*a;
    Xr = real(x_hat);  Xi = imag(x_hat);

    R(Xr<b(1)) = c(1);  I(Xi<b(1)) = c(1);
    R(b(1)<=Xr&Xr<b(2)) = c(2);  I(b(1)<=Xi&Xi<b(2)) = c(2);
    R(b(2)<=Xr&Xr<b(3)) = c(3);  I(b(2)<=Xi&Xi<b(3)) = c(3);
    R(b(3)<=Xr) = c(4);  I(b(3)<=Xi) = c(4);

    X_hat = R + 1i*I;

    dis = abs(s - X_hat.');
    ser(dis>0.01)=1;
    ser(dis<=0.01)=0;
    Ser(i_snr) = sum(ser)/(num*16);
    
    dis = abs(X_hat.' - QAM_16.');
    onehot_p(dis>0.01) = 0;
    onehot_p(dis<=0.01) = 1;
    onehot_p = reshape(onehot_p, [num*16,16]);
    Ber(i_snr) = sum(sum(abs(onehot_p*BIT_16 - onehot_s*BIT_16)))/(num*16*4);
end
Ser;
Ber;

end

