Calculating the Bit error rate (BER) for a Visible Light Communication (VLC)

I am trying to calculate the bit error rate (BER) for a visible light communication LED source.

Based on the code of this link Number of bit errors and bit error rate (BER) I calculate the BER. Unfortunately, the results that I got were unreasonable as shown below:

but according to this link https://fr.mathworks.com/help/comm/ref/biterr.html?s_tid=srchtitle_biterr_1

their results like below:

Why is my estimated BER so far off from the theoretical one?

Is there anyone who calculated the BER for a VLC before? May I get an assistance, please?

close all;
clear variables;
clc;

%% Simulation Parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    Main Simulation Parameters                     %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------%
% NUMBER OF LIGHT SOURCES %
%-------------------------%
N_t = 1; % Number of light sources 

% Set the simulation parameters.
M = 64;                 % Modulation order
k = log2(M);            % Bits per symbol
EbNoVec = (5:15)';      % Eb/No values (dB)
numSymPerFrame = 100;   % Number of QAM symbols per frame

n = 34560; % Number of transmitted bits
N_iter = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          AP Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------%
% LIGHT SOURCES GEOMETRY %
%------------------------%
L = 20; W = 20; H = 3; % Length, width and height of the room (m)
theta_half = 30;
m = -log(2)./log(cosd(theta_half)); % Lambertian order of emission
coord_t = [0 0 0]; % Positions of the light sources

n_t_LED = [0, 0, -1]; n_t_LED = n_t_LED/norm(n_t_LED); % Normalized normal vector of each light source
n_t = repmat(n_t_LED, N_t, 1); % Normalized normal vectors of the light sources

%-------------------------------------%
% LIGHT SOURCES ELECTRICAL PARAMETERS %
%-------------------------------------%
P_LED = 2.84; % Average electrical power consummed by each light source (W)

param_t = {coord_t, n_t, P_LED, m};

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          Rx Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------%
% PHOTODETECTOR PARAMETERS %
%--------------------------%
A_det = 1e-4;        % Photoreceiver sensitive area (m²)
FOV = 60*pi/180;     % Fielf-of-view of the photoreceiver

T_s = 1;                      % Gain of the optical filter (ignore if not used)
index = 1.5;                  % Refractive index of the Rx concentrator/lens (ignore if not used)
G_Con = (index^2)/sin(FOV);   % Gain of an optical concentrator; ignore if no lens is used

n_r = [0, 0, 1];     % Normal vector of the photoreceiver
n_r = n_r/norm(n_r); % Normal vector of the photoreceiver (normalized)

%---------------------------%
% RECEIVER PLANE PARAMETERS %
%---------------------------%
step = 0.5; % Distance between each receiving point (m)
X_r = -L/2:step:L/2; % Range of Rx points along x axis
Y_r = -W/2:step:W/2; % Range of Rx points along y axis
N_rx = length(X_r); N_ry = length(Y_r); % Number of reception points simulated along the x and y axis
z_ref = 0.85; % Height of the receiver plane from the ground (m)
z = z_ref-H; % z = -1.65; % Height of the Rx points ("-" because coordinates system origin at the center of the ceiling)
if( abs(z) > H )
    fprintf('ERROR: The receiver plane is out of the room.\n');
    return
end

param_r = {A_det, n_r, FOV}; % Vector of the Rx parameters used for channel simulation

berEst = zeros(size(EbNoVec));

 for n = 1:length(EbNoVec)
     
    snrdB = EbNoVec(n) + 10*log10(k);
    numErrs = 0;
    numBits = 0;

    while numErrs < 200 && numBits < 1e7
        data = randi([0 1],numSymPerFrame,k);
        dataSym = bi2de(data);
        txSignal = qammod(dataSym,M);
        % LOS received optical power calculation
        H0_LOS = zeros(N_rx,N_ry,N_t);
        P_rec_dbm = zeros(N_rx,N_ry,N_t);
        %rxSignal = zeros(N_t,length(txSignal));
        T = param_t{1}(1,:);
        P_t = param_t{3};
        for r_x = 1:N_rx
            for r_y = 1:N_ry
                for i_t = 1:N_t
                    x = X_r(r_x); y = Y_r(r_y);
                    R = [x,y,z];
                    v_tr = (R-T)./norm(R-T);
                    d_tr = sqrt(dot(R-T,R-T));
                    phi = 0;
                    psi = 0;
                    H0_LOS(r_x,r_y,i_t)= (m+1)*A_det/(2*pi*d_tr^2)*cosd(phi)^m*cosd(psi);
                end
            end
        end
        rx = (txSignal).*(H0_LOS(r_x,r_y,i_t)*G_Con*T_s); % Photocurrent produced by the PD without noise (A)
        rxSignal = awgn(rx,snrdB,'measured');
        rxSym = qamdemod(rxSignal,M);
        dataOut = de2bi(rxSym,k);
        nErrors = biterr(data,dataOut);
        
        numErrs = numErrs + nErrors;
        numBits = numBits + numSymPerFrame*k;
    end
        % Estimate the BER
    berEst(n) = numErrs/numBits;
 end

%Determine the theoretical BER curve by using the berawgn function.
berTheory = berawgn(EbNoVec,'qam',M);
P_r_LOS = P_t.*H0_LOS.*T_s.*G_Con;
P_rec_dBm = 10*log10(P_r_LOS*1000);

% Plot the estimated and theoretical BER data. The estimated BER data 
% points are well aligned with the theoretical curve.
figure
semilogy(EbNoVec,berEst,'*')
hold on
semilogy(EbNoVec,berTheory)
grid
legend('Estimated BER','Theoretical BER')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')

[kim zheng cho]

Hi may i know how you eventually solve the problem?

[Haitham AL Satai]

@kim zheng cho I did not solve it yet dear. Do you work on the same project?

[kim zheng cho]

ya i am working on applying fbmc modulation into vlc circuit hardware circuit but i havent apply yet but rxsignal from scope has high ber around 0.49 so i also wonder how to solve it. Then may i ask do you identify why the ber is high ?what is the reason?

[Haitham AL Satai]

@kim zheng cho This is my first year in my study, and I did not go further in my work. To be honest with you, I did not solve it and I could not find any solution in Internet. But I promise you, if I find any solution I will send it to you. Let us keep in touch.

[kim zheng cho]

that's very kind of you.Likewise if i have done i will send you as well. Wish you good luck......

[Haitham AL Satai]

@kim zheng cho You are most welcome. U can send me an email if you want [e-mail].