MMSE Channel Estimation for Comb type pilots

I am working on a project for minimum mean square error (MMSE) channel estimation using comb-type pilot aided for OFDM system. I wrote the following code but I could not get the required results. Where I get the BER of the LS identical to the MMSE BER!!!

clear; clc; 
Nsc = 256; % OFDM symbol size (Number of subcarriers).
M = 16; % Modulation order
Nsmb = 100; % Number of OFDM symbols to be simulated
Ne = 100; % Number of bits in error
Ncp = Nsc/4;% Cyclic Prefix Length
Lh = 4; % Channel Length
Lp = 4;  % pilot interval
Np  = Nsc/Lp; % total number of Pilots
% Pilot Locations
      PL = 1:Lp:Nsc;          % location of pilots
      DL = setxor(1:Nsc,PL);  % location of data
      Pp = 4;                 % Pilot Power
str=-10;
stp=4; %SNR starts at -10 dB, with step size [dB] = stp 
Esnr=40; % Last value of SNR [dB] = Esnr
c = 0; % steps counter
NSteps = length(str:stp:Esnr); 
EbNo_vector=str:stp:Esnr;
berRslt_S=zeros(size(EbNo_vector)); 
berRslt_LS=zeros(size(EbNo_vector)); 
berRslt_MMSE=zeros(size(EbNo_vector));
for snr=str:stp:Esnr
    c =c+1;
%%Monte Carlo Simulation loop
disp(['STEP ',num2str(c),' of ',num2str(length(str:stp:Esnr)),' :Processing SNR = ',num2str(snr)]);
nEr_S = 0; % Number of collected errors
nEr_LS = 0; % Number of collected errors
nEr_MMSE = 0; % Number of collected errors
nSmb = 0; % Number of simulated OFDM symbols
while ((nEr_S < Ne) && (nSmb < Nsmb))
 % Random Channel Generation
        h= randn(1,Lh) + 1j * randn(1,Lh);
        h = h./norm(h); % normalization 
        H = fft(h.',Np);
    % Transmitter
    Dg=randi([0 M-1],1,Nsc); % Data Generation
    Dmod=qammod(Dg,M,'UnitAveragePower',true); % Baseband modulation (mapping)
    % Serial To Parallel
    Dmod_P=Dmod.';
    Dmod_P2=Dmod_P;
    % Pilot Insertion
    Dmod_P2(PL) = Pp * Dmod_P2(PL);
    % Implementing IFFT 
    d_mod_P=(1/sqrt(Nsc))*ifft(Dmod_P,Nsc);
    d_mod_P2=(1/sqrt(Nsc))*ifft(Dmod_P2,Nsc);
    % Parallel To Serial
    dAM_mod_S=d_mod_P.';
    dAM_mod_S2=d_mod_P2.';
    % Cyclic Prefix Insertion (GII)
    d_mod_S_GII = [dAM_mod_S(Nsc- Ncp + 1 : Nsc),dAM_mod_S];
    d_mod_S_GII2 = [dAM_mod_S2(Nsc- Ncp + 1 : Nsc),dAM_mod_S2];
    % Passing Through the Channel
    d_mod_S_GII_ch = filter(h,1,d_mod_S_GII2);    % channel effect
    % Adding Noise using AWGN
    d_mod_noisy_S_GII=awgn(d_mod_S_GII,snr,'measured');
    d_mod_noisy_S_GII_ch=awgn(d_mod_S_GII_ch,snr,'measured');
    % Cyclic Prefix Removal (GIR)
    d_mod_noisy_S_GIR=d_mod_noisy_S_GII(Ncp+1:Nsc+Ncp);
    d_mod_noisy_S_GIR_ch=d_mod_noisy_S_GII_ch(Ncp+1:Nsc+Ncp);
    % Serial To Parallel
    d_mod_noisy_P_GIR=d_mod_noisy_S_GIR.';
    d_mod_noisy_P_GIR_ch=d_mod_noisy_S_GIR_ch.';
    % Amplitude demodulation (DFT using fast version FFT)
    amdemod_P_GIR=(sqrt(Nsc))*fft(d_mod_noisy_P_GIR,Nsc);
    amdemod_P_GIR_ch=(sqrt(Nsc))*fft(d_mod_noisy_P_GIR_ch,Nsc);
    %%Channel Estimation
    % Extracting received pilots
    TxPilots = Dmod_P2(PL); % trnasmitted pilots
    RxPilots = amdemod_P_GIR_ch(PL); % received pilots
    % LS Channel Estimation
    Hp_LS= RxPilots./TxPilots; 
    % MMSE Channel Estimation
    R_HH = H*H';%toeplitz(H); 
    XX = toeplitz(TxPilots);     
    powerDB = 10*log10(var(RxPilots)); % Calculate Tx signal power
    sigmI = 10.^(0.1*(powerDB)); % Calculate the noise variance
    G=(R_HH)*(R_HH+(1/sigmI)*(XX))^(-1);
    Hmmse=((G)*(Hp_LS));
    % Channel Estimation at the Data Subcarriers
    HData_LS=interp1(PL,Hp_LS,1:Nsc,'spline');
    HData_MMSE = interp1(PL,Hmmse,1:Nsc,'spline');
    %%Equalization
    Deq_LS =((amdemod_P_GIR_ch.')./HData_LS); 
    Deq_MMSE =((amdemod_P_GIR_ch.')./HData_MMSE);
    % Detection (De-Mapping)
    y_P=qamdemod(amdemod_P_GIR,M,'UnitAveragePower',true);
    y_P_Eqz_LS=qamdemod(Deq_LS.',M,'UnitAveragePower',true);
    y_P_Eqz_MMSE=qamdemod(Deq_MMSE.',M,'UnitAveragePower',true);
    % Parallel To Serial
    y_S=y_P.';
    y_S_Eqz_LS=y_P_Eqz_LS.';
    y_S_Eqz_MMSE=y_P_Eqz_MMSE.';
%%Removing Pilots from received data and original data 
    D_no_pilots=Dg(DL); % removing pilots from D
    Rec_d_LS=y_S_Eqz_LS(DL); % removing pilots from d_received_chann_LS    
    Rec_d_MMSE=y_S_Eqz_MMSE(DL); % removing pilots from d_received_chann_MMSE
    % BER Calculation
    [n_S, r_S]=biterr(Dg,y_S);
    [n_LS, r_LS]=biterr(D_no_pilots,Rec_d_LS);
    [n_MMSE, r_MMSE]=biterr(D_no_pilots,Rec_d_MMSE);
    nEr_S=nEr_S+r_S;
    nEr_LS=nEr_LS+r_LS;
    nEr_MMSE=nEr_MMSE+r_MMSE;
    nSmb=nSmb+1;
end
berRslt_S(c)=nEr_S/(log2(M)*nSmb);
berRslt_LS(c)=nEr_LS/(log2(M)*nSmb);
berRslt_MMSE(c)=nEr_MMSE/(log2(M)*nSmb);
end
figure;
snr=str:stp:Esnr;
semilogy(snr,berRslt_S,'-b','linewidth',1,'markerfacecolor','b','markersize',8,'markeredgecolor','b');grid;hold on;
semilogy(snr,berRslt_LS,'-k','linewidth',1,'markerfacecolor','k','markersize',8,'markeredgecolor','k');
semilogy(snr,berRslt_MMSE,'-g','linewidth',1,'markerfacecolor','g','markersize',8,'markeredgecolor','g');
title('OFDM Bit Error Rate vs SNR');
ylabel('Bit Error Rate');
xlabel('SNR [dB]');
legend('FFT','LS','MMSE');

Answer 1

[Abby__DSP]

Hello,

Can you explain this line : powerDB = 10*log10(var(RxPilots)); % Calculate Tx signal power

it is not ok with the comment you made

Answer 2

[Sami Hadeyya]

Does it work well this code?

please reply.