matlab code for frequency selective fading

I have this code- clc; clear all; close all; N = 128; % No of subcarriers Ncp = 16; % Cyclic prefix length Ts = 1e-3; % Sampling period of channel %Ts = 1e-6; Fd = 1; % Max Doppler frequency shift Np = 4; % No of pilot symbols M = 2; % No of symbols for PSK modulation Nframes = 10^3; % No of OFDM frames %D = round((M-1)*rand((N-2*Np),Nframes)); %D=rand(1,Nframes); const = pskmod([0:M-1],M); r=3.56995; %r=2.9; n=120; n1=1000; m0=.6; m=zeros(n,n1); m(1,1)=m0; for i=1:119999 m(i+1)=r*m(i)*(1-m(i)); end w0=.6; w(1)=w0; for k=1:127 w(k+1)=r*w(k)*(1-w(k)); end D = round((M-1)*rand((N-2*Np),Nframes));

Dmod2 = pskmod(D,M); Dmod1=real(Dmod2);

%m1=m.'; Dmod=Dmod1.*m;

Data = [zeros(Np,Nframes); Dmod ; zeros(Np,Nframes)]; % Pilot Insertion

%% OFDM symbol

IFFT_Data = (128/sqrt(120))*ifft(Data,N); TxCy = [IFFT_Data((128-Ncp+1):128,:); IFFT_Data]; % Cyclic prefix [r c] = size(TxCy); Tx_Data = TxCy;

%% Frequency selective channel with 4 taps

tau = [0 1e-5 3.5e-5 12e-5]; % Path delays %tau = [0 1e-2 3.5e-2 12e-2]; pdb = [0 -1 -1 -3]; % Avg path power gains h = rayleighchan(Ts, Fd); h.StoreHistory = 0; h.StorePathGains = 1; h.ResetBeforeFiltering = 1;

%% SNR of channel

EbNo = -7:20; EsNo= EbNo + 10*log10(120/128)+ 10*log10(128/144); % symbol to noise ratio snr= EsNo - 10*log10(128/144); Ec_N0_dB = EbNo - 10*log10(4);

%% Transmit through channel

berofdm = zeros(1,length(snr)); Rx_Data = zeros((N-2*Np),Nframes); for i = 1:length(snr) for j = 1:c % Transmit frame by frame hx = filter(h,Tx_Data(:,j).'); % Pass through Rayleigh channel a = h.PathGains; AM = h.channelFilter.alphaMatrix; g = a*AM; % Channel coefficients G(j,:) = fft(g,N); % DFT of channel coefficients y = awgn(hx,snr(i)); % Add AWGN noise

%% Receiver

        Rx = y(Ncp+1:r);                                % Removal of cyclic prefix 
        FFT_Data = (sqrt(120)/128)*fft(Rx,N)./G(j,:);   % Frequency Domain Equalization
        FFT=FFT_Data.*w;
        Rx_Data(:,j) = pskdemod(FFT(5:124),M);     % Removal of pilot and Demodulation 
        %Rx_Data1= Rx_Data.*m;

    end
    berofdm(i) = sum(sum(Rx_Data~=D))/((N-2*Np)*Nframes);
end

%% Plot the BER EbN0Lin = 10.^(EbNo/10); theoryBer = 0.5.*(1-sqrt(EbN0Lin./(EbN0Lin+1)));

figure; semilogy(EbNo,theoryBer,'bs-','linewidth',2); hold on semilogy(EbNo,berofdm,'mx-','linewidth',2); axis([-7 20 10^-3 2]) grid on; title('OFDM BER vs SNR in Frequency selective Rayleigh fading channel'); xlabel('EbNo'); ylabel('BER');

How to get the correct BER? Plz help me out.