Trying to analyze the BER for DCO-OFDM for Underwater VLC. The resulting curve is not convincing, something wrong with the way I have integrated the channel model, I guess.

clc ;

clear all ;

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)

Ncp =256 ; % Length of the Cyclic Prefix

Data = randi ([0 M-1] ,m , N ); % Generation of Random bits Matrix of size m by N

DataMod = qammod( Data , M);

DataMod_serialtoparallel = DataMod .'; % Performing Serial to Parallel Conversion

datamat = DataMod_serialtoparallel ; % Assigning the total data 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 complex conjugate versions of the previous ones .

d_ifft = ifft (( datamat ) ) ; % Computation of IFFT operation

d_ifft_paralleltoserial = d_ifft .'; % Parallel to Serial Conversion

CP_part = d_ifft_paralleltoserial (: ,end - Ncp +1: end) ; % Addition of Cyclic Prefix

DCOOFDM_CP =[ CP_part d_ifft_paralleltoserial ]; %Transmissin of DCO - OFDM signal

%% VLC channel modelling

theta=70;

% semi-angle at half power

m=-log10(2)/log10(cosd(theta));

%Lambertian order of emission

P_total=20;

%transmitted optical power by individual LED

Adet=1e-2;

%assumed a custom made PD with a detector physical area of 1cm^2

Ts=1;

%assumed an ideal optical filter

index=1.5;

%refractive index of a lens at a PD, value supported by paper 3.

C = 0.1162;

%underwater attenuation coefficient, value was taken from paper 3.

FOV=60*pi/180;

%FOV of a receiver

G_Con=(index^2)/sin(FOV);

%gain of an optical concentrator; ignore if no lens is used

lx=5; ly=5; lz=10;

% Underwater zone dimension in metre

h=10;

%the distance between source and receiver plane

XT=0; YT=0;

% position of LED;

Nx=lx*10; Ny=ly*10;

% number of grid in the receiver plane

x=-lx/2:lx/Nx:lx/2;

y=-ly/2:ly/Ny:ly/2;

[XR,YR]=meshgrid(x,y);

% receiver plane grid

D1=sqrt((XR-XT(1,1)).^2+(YR-YT(1,1)).^2+h^2);

% distance vector from source 1

cosphi_A1=h./D1;

% angle vector

H_A1=(m+1)*Adet.*cosphi_A1.^(m)./(2*pi.*D1.^2);

% channel DC gain for source 1

H_uw=H_A1.*exp(-C*D1);

% underwater channel DC gain for source 1; equation was supported 1,2,3

H_A2 = H_uw ./ norm ( H_uw );

d_channell = filter( H_A2(1,1:2),1,DCOOFDM_CP .').';

%%

bdc =7; %DC BIAS

clip = sqrt ((10.^( bdc /10) ) -1) ; % clipping factor k

bdcc = clip * sqrt (d_channell .* d_channell ) ;%Computation of DC bias

DCOOFDM_BIAS = bdcc + d_channell ; % Addition of DC bias to the cyclic prefix added 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 ;

DCOOFDM_with_chann = awgn ( DCOOFDM_BIAS , SNR ,'measured'); % Addition of AWGN

% Receiver of DCO - OFDM

DCOOFDM_with_chann1 = DCOOFDM_with_chann - bdcc ; %Removal of DC bias

DCOOFDM_removal_CP = DCOOFDM_with_chann1 (: , Ncp +1: N + Ncp ) ; % Removal of Cyclic Prefix

DCOOFDM_serialtoparallel = DCOOFDM_removal_CP .'; %Serial to Parallel Conversion

DCOOFDM_parallel_fft = fft (DCOOFDM_serialtoparallel) ; % Computation of FFT operation

DCOOFDM_Demodulation = qamdemod (DCOOFDM_parallel_fft .', M );

[~ , s_e1(count)]= symerr ( Data (: ,2:512),DCOOFDM_Demodulation (: ,2:512)) ;

end

% Plotting the BER curves

semilogy ( snr_vector , s_e1 ,'rd -','LineWidth' ,2);

legend('FFT - based DCO -OFDM -7 dB of Bias') ;

axis ([0 30 10^-4 1]) ;

xlabel ('SNR in dB');

ylabel ('BER');

grid on;