Why am I getting a poor BER vs SNR curve?

close all ;
m =512; % Total number of OFDM symbols
N =1024; % Length of each OFDM symbol
M =4; % Size of the Constellation ( M can be 4 , 8 , 16 , 32 , 64 , 128 , 256)
pilotFrequency =8 ; % Pilot Symbol insertion frequency
E =2; % Energy of each Pilot Symbol
Ncp =256 ; % Length of the Cyclic Prefix
%Tx = qammod ('M',M ) ; % Choosing the modulation format as M - ary Quadrature Amplitude Modulation (M- QAM )
%Rx = qamdemod ('M',M) ; % Fixing up the demodulation format at the receiving end as M - ary QAM
% ACO - OFDM Transmitter
Data = randi ([0,M-1] ,m , N ); % Generation of Random bits Matrix of size m by N
for k1 =1: m
    for m1 =1: N
        if mod ( m1 ,2) ==0 % Performing Modulo operation to extract the even subcarriers
            Data ( k1 , m1 ) =0; % Setting the Even Subcarriers to Zero
        end
    end
end
DataMod = qammod( Data , M);
%DataMod = modulate ( Tx , Data ) ; % Performing Data Modulation
DataMod_serialtoparallel = DataMod .'; % Performing Serial to Parallel Conversion
PLoc = 1: pilotFrequency :N ; % Fixing the location of Pilot carrires
DLoc = setxor (1: N , PLoc ) ; % Fixing the location of Data subcarriers
DataMod_serialtoparallel ( PLoc ,:) = E *DataMod_serialtoparallel ( PLoc ,:) ; % Inserting Pilot carriers
datamat = DataMod_serialtoparallel ; % Assigning the total data including the pilots to a variable called datamat
% Computation of Hermitian Symmetry Criteria
datamat (1 ,:) =0; % Assigning the First subcarrier to Zero
datamat (513 ,:) =0; % Assigning the Middle Subcarrier to Zero
datamat (514:1024 ,:) = flipud ( conj ( datamat (2:512 ,:) )) ; %Illustrating that only half of the subcarriers are exploited for data transmission as the remaining half are flipped omplex conjugate versions of the previous ones .
d_ifft = ifft (( datamat ) ) ; % Computation of IFFT operation
% Ensuring that only the positive portion of the signal is transmitted
for k2 =1: N
    for m2 =1: m
        if( d_ifft ( k2 , m2 ) <0)
            d_ifft ( k2 , m2 ) =0;
        end
    end
end
d_ifft_paralleltoserial = d_ifft .'; % Parallel to Serial Conversion
CP_part = d_ifft_paralleltoserial (: ,end - Ncp +1: end) ; % Addition of Cyclic Prefix
ACOOFDM_CP =[ CP_part d_ifft_paralleltoserial ]; % Transmissin of ACO - OFDM signal
% VLC Channel Modeling
theta = 70; %LED semi - angle
ml = - log10 (2) / log10 ( cos ( theta ) ) ; % Computation of Lambertian Mode Number
APD =0.01; % Area of the photodiode
lx =5; ly =5; lz =3; % Size of the Dimensions of the Indoor Room Environment
h =2.15;
[ XT , YT ]= meshgrid ([ - lx /4 lx /4] ,[ - ly /4 ly /4]) ;
Nx = lx *5; Ny = ly *5;
x = linspace ( - lx /2 , lx /2 , Nx ) ;
y = linspace ( - ly /2 , ly /2 , Ny ) ;
[ XR , YR ]= meshgrid (x , y ) ;
D1 = sqrt (( XR - XT (1 ,1) ) .^2+( YR - YT (1 ,1) ) .^2+ h ^2) ; D2 = sqrt (( XR - XT (2 ,2) ) .^2+( YR - YT (2 ,2) ) .^2+ h ^2) ;
cosphi_A1 = h ./ D1 ;
receiver_angle = acosd ( cosphi_A1 ) ;
H_A1 =3600*(( ml +1) * APD .* cosphi_A1 .^( ml +1) ./(2* pi .* D1.^2) +( ml +1) * APD .* cosphi_A1 .^( ml +1) ./(2* pi .^2* D1.^2* D2 .^2) ) ; % Computation of VLC channel Impulse Response taking into consideration Non Line of Sight ( NLOS ) Environment
H_A2 = H_A1./ norm ( H_A1 );
d_channell = filter ( H_A2 (1 ,1:2) ,1 , ACOOFDM_CP .') .';
% Illustration of channel effect on the transmitted ACO - OFDM signal
count =0;
snr_vector =0:1:30; % size of signal to noise ratio (SNR ) vector
for snr = snr_vector
    SNR = snr + 10* log10 ( log2 (M ) ) ;
    count = count +1 ;
    ACOOFDM_with_chann = awgn(d_channell , SNR ,'measured') ; % Addition of AWGN
    % Receiver of ACO - OFDM
    ACOOFDM_removal_CP = ACOOFDM_with_chann (: , Ncp +1: N + Ncp );
    % Removal of Cyclic Prefix
    ACOOFDM_serialtoparallel = ACOOFDM_removal_CP .'; %Serial to Parallel Conversion
    ACOOFDM_parallel_fft = fft ( ACOOFDM_serialtoparallel) ; 
    % Computation of FFT operation 
    % Channel Estimation
    TransmittedPilots = DataMod_serialtoparallel ( PLoc,:) ; % Extracting the transmitted pilots
    ReceivedPilots = ACOOFDM_parallel_fft ( PLoc ,:) ; 
    %Extracting the received pilot tones effected by channel
    H_LS = ReceivedPilots ./ TransmittedPilots ; % Least Square Channel Estimation
    for r =1:m
        H_MMSE(: , r ) = MMSE(ReceivedPilots (: , r ) ,TransmittedPilots (: , r ), N, pilotFrequency, H_A2 (1 ,1:2), SNR);% Minimum Mean SquareError ( MMSE ) Channel Estimation
    end
    for q =1:m
        HData_LS(: , q ) = interpolate ( H_LS (: , q ) .', PLoc ,N ,'spline ') ; % Interpolation
    end
    HData_LS_parallel1 = HData_LS .'; % Parallel to Serial Conversion
    HData_MMSE_parallel1 = H_MMSE .';
    ACOOFDM_SERIAL_LS = demodulate ( Rx ,(ACOOFDM_parallel_fft .') ./ HData_LS_parallel1 ) ; % Demodulation
    ACOOFDM_SERIAL_MMSE = demodulate (Rx ,(ACOOFDM_parallel_fft .') ./( HData_MMSE_parallel1) ) ;
    % Recovery of Pilots from the Original Transmitted signal and Received Signal
    Data_no_pilots = Data (: , DLoc ) ;
    Recovered_Pilot_LS = ACOOFDM_SERIAL_LS (: , DLoc ) ;
    Recovered_Pilot_MMSE = ACOOFDM_SERIAL_MMSE (: , DLoc ) ;
    % Computation of Bit Error Rate
    [~ , recoveredLS( count ) ]= biterr ( Data_no_pilots(: ,2:255) , Recovered_Pilot_LS (: ,2:255) ) ;
    [~ , recoveredMMSE( count ) ]= biterr ( Data_no_pilots (: ,2:255) , Recovered_Pilot_MMSE (: ,2:255) ) ;
end

% Plotting the BER curves
semilogy ( snr_vector , recoveredLS ,'rd -','LineWidth ' ,2) ;
hold on ;
semilogy ( snr_vector , recoveredMMSE ,'gs -','LineWidth ',2) ;
axis ([0 30 10^ -4 1]) ;
grid on ;

function H_MMSE = MMSE(RxP,TxP,N,pilotFrequency,h_CIR,SNR)
% I modified the MMSE_CE function provided in :MIMO-OFDM Wireless 
% Communications with MATLAB¢ç   Yong Soo Cho, Jaekwon Kim, 
% Won Young Yang and Chung G. Kang
%2010 John Wiley & Sons (Asia) Pte Ltd
% The modification made the function more suitable with my OFDM Matlab code
noiseVar = 10^(SNR*0.1);
Np=N/pilotFrequency; % Number of Pilots
H_LS = RxP./TxP;  % LS estimate
k=0:length(h_CIR)-1; 
hh = h_CIR*h_CIR'; 
% tmp = h_CIR.*conj(h_CIR).*k;
r = sum(h_CIR.*conj(h_CIR).*k)/hh;  
r2 = (h_CIR.*conj(h_CIR).*k)*k.'/hh;
t_rms = sqrt(r2-r^2);     % rms delay
D = 1j*2*pi*t_rms/N; % Denomerator of Eq. (6.16) page 192
K1 = repmat([0:N-1].',1,Np);
K2 = repmat([0:Np-1],N,1);
rf = 1./(1+D*(K1-K2*pilotFrequency));
K3 = repmat([0:Np-1].',1,Np);
K4 = repmat([0:Np-1],Np,1);
rf2 = 1./(1+D*pilotFrequency*(K3-K4));
Rhp = rf;
Rpp = rf2 + eye(length(H_LS),length(H_LS))/noiseVar;
H_MMSE = transpose(Rhp*inv(Rpp)*H_LS);  % MMSE channel estimate

    function [H_interpolated] = interpolate(H,pilot_loc,Nfft,method)
        % Input: H = Channel estimate using pilot sequence
        % pilot_loc = Location of pilot sequence
        % Nfft = FFT size
        % method = ’linear’/’spline’
        % Output: H_interpolated = interpolated channel
        if pilot_loc(1)>1
            slope = (H(2)-H_est(1))/(pilot_loc(2)-pilot_loc(1));
            H = [H(1)-slope*(pilot_loc(1)-1) H]; pilot_loc = [1 pilot_loc];
        end
        if pilot_loc(end) <Nfft
            slope = (H(end)-H(end-1))/(pilot_loc(end)-pilot_loc(end-1));
            H = [H H(end)+slope*(Nfft-pilot_loc(end))];
            pilot_loc = [pilot_loc Nfft];
        end
        if lower(method(1))=='l', H_interpolated=interp1(pilot_loc,H,[1:Nfft]);
        else
            H_interpolated = interp1(pilot_loc,H,[1:Nfft],'spline');
        end
    end
end