problem with BER for OFDM estimation for wimax please help

N    = 128;   % fft size
Pil  = N/4; % 32
S    = N-Pil;   % number of data subcarriers  112
nIteration = 500;
L    = 30;
GI   = N/4; % 32
Ep   = 2;
SNR_V = [0:3:27];
ber = zeros(1,length(SNR_V));   % initializing bit error rate

M = 2; % modulation object

% Generating data

T_Data=randint(N,1,M); %  128  1

% modulating data

mod_data = pskmod(T_Data,M); 

% serial to parallel conversion
par_data = reshape(mod_data,8,16); %    8 16
Ip=1:4:N;
Is = setxor(1:N,Ip); 
% pilot insertion
 pilot_ins_data=[];
 for Ip=1:4:13
     pilots= Ep*par_data(:,Ip);
     xp=[pilots  par_data(:,Ip+1) par_data(:,Ip+2) par_data(:,Ip+3)];
     pilot_ins_data=[ pilot_ins_data  xp];
 
 end

 pilot_ins_data; % 8 16
 PilotData=reshape(pilot_ins_data,8*16,1); % 128 1
 LocPilot= [ PilotData(1:8,1); PilotData(33:40,1); PilotData(65:72,1); PilotData(97:104,1) ]; % Locations of pilots after reshape 32 1
 
 % fourier transform time doamain data and normalizing the data

no_of_error=[];
ratio=[];
for i = 1 : length(SNR_V)
    SNR = SNR_V(i)
    for k = 1 : nIteration
        
        IFFT_data = ifft(PilotData,N); % 128 1

        % addition cyclic prefix

        Cylic_add_Data = [IFFT_data(N- GI + 1 : N); IFFT_data]; % 96 1

        % generating random channel coefficients
        h(1:L,1)  =     random('norm',0,1,L,1) + ...
                        j * random('norm',0,1,L,1);  
        h  = h./sum(abs(h));    % normalization   30 1

        ChannEffect = filter(h,1,Cylic_add_Data);    % channel effect  160 1

        % adding awgn noise
        Signal_Noise = awgn(ChannEffect , SNR - db(std(ChannEffect))); % normalization to signal power


        Cyclic_pre_Rem =  Signal_Noise(GI+1:N+GI);  %cyclic prefix removal 
 
        FFT_recdata =fft(Cyclic_pre_Rem,N);    % freq domain transform  128 1 

       % LS Estimation & 1D interpolation
       
           F = exp(2*pi*sqrt(-1)/N .* meshgrid([0:N-1],[0:N-1])...
                           .* repmat([0:N-1]',[1,N])); % 128 128
      
       Tpilot= LocPilot;% transmitted pilots 32 1
       
       Rpilot=[ FFT_recdata(1:8,1);   FFT_recdata(33:40,1);   FFT_recdata(65:72,1);  FFT_recdata(97:104,1) ]; % received pilots 32 1
     
      % Hpls=inv(diag(LocPilot))*  Rpilot; % 32 1
      
       Ipp=[ 1:8 33:40 65:72 97:104 ]; 
       
       G = (Ep * length(Ipp))^-1 ...
                         * ctranspose(sqrt(Ep)*diag(Tpilot)*ctranspose(F(1:L,Ipp))); % 30 32
       Hpls = G*Rpilot; % 30 1
              
       Hh=interp(Hpls,4); % 120 1
 
       Demod_Data = pskdemod(  FFT_recdata./(fft(Hpls,N)),M);  %demodulating the data  128 1 
       % or  Demod_Data = pskdemod(  (fft(Hh,N)./(fft(Hpls,N)),M);

       DataEstiRec=[  Demod_Data(9:32,1); Demod_Data(41:64,1); Demod_Data(73:96,1); Demod_Data(105:128,1) ];
       DataTra=[  T_Data(9:32,1); T_Data(41:64,1); T_Data(73:96,1); T_Data(105:128,1) ];
       
       [nErr bErr(i,k)] = symerr(DataEstiRec ,  DataTra);

        
        ms_error_mat=mean(((abs(h-Hpls))/abs(h)).^2);
 for v=1:16
    if(ms_error_mat(v)~=0)
        ms_error=ms_error_mat(v);
    end
end
ls_mse(k,i)=ms_error;
ls_mse_ave=mean(ls_mse);
semilogy(SNR,ls_mse_ave,'--ro')
grid on
 end
 end



% plotting the result

xlabel('SNR in dB')
ylabel('mean squared error')
title('PLOT OF SNR V/S MSE FOR AN OFDM SYSTEM WITH LS ESTIMATOR BASEDRECEIVERS')
figure(2)
semilogy(SNR_V,mean(bErr'),'--or','linewidth',2);
hold on;
EbN0Lin = 10.^(SNR_V/10);
theoryBer = 0.5.*(1-sqrt(EbN0Lin./(EbN0Lin+1)));
semilogy(SNR_V,theoryBer,'--ob','linewidth',2);
legend('simulated','theoritical')
grid on

xlabel('SNR');
ylabel('BER')
title('Bit error probability curve for BPSK using OFDM');