[SOLVED] LS and MMSE channel estimation for ofdm

[Romance42932]

hi...
i want to pass tx (QPSK) signal thru a awgn channel.... let say i have provided snr of 30dB and channel pdp with L=4(tap length) is exp(-n/5) where n is from(0 to 4)...
so i want to create channel with say length 5 with noise which satisfies given snr in time domain.... so that y = hx + n then i wanna do channel estimation in freq domain Y = HX + N.....
tell me also what is the equivalent variance of that noise in freq domain? (for N).... i am doing this code.... i dont know what is wrong in it...

i am using pilots as +1 and -1 ......i am trying find H values on all subcarriers .....not only on pilot....is that make any difference


%%%%%%%%%%CODE%%%%%%%%%%%%%%%
clear all
clc
sub_car = 128; % number of data subcarriers
nbitpersym = 128;
len_fft = 128;
t_pilot = 128/8;
pilotinterval = 8;
nsym = 10^4; % number of symbols
%Create a QPSK modulator System object with bits as inputs and Gray-coded signal constellation
h = comm.QPSKModulator('PhaseOffset',pi/16);
% change the phase offset to avoid same types of values
% Generating data
t_data=randint(nbitpersym,1);
% modulating data
mod_data = step(h,t_data);
% pilot insertion
mod_data(1:2*pilotinterval:end)= 1; % pilot is taking as +1 or -1
mod_data(9:2*pilotinterval:end)= -1;
IFFT_data = ifft(mod_data.',len_fft);
a = IFFT_data([1:8:end]);
% addition cyclic prefix
cyclic_add_data = [];
cyclic_add_data = [IFFT_data,[113:128]) IFFT_data];
noise_power = [];
noise_var = [];
pdp = exp(-(0:5)/5);
snr = 30; % snr in dB
% for i= 1:length(snr)
noise_power = 2./10.^(snr./10);
noise_var = sqrt(noise_power);
% end
pdp = exp(-(0:4)/5);
hx = pdp*randn(1,5)';
channel = hx + sqrt(noise_var)*randn(1);
%passing thru the channel
rec_sig = conv(cyclic_add_data,channel);
%discarding 149-144 = 5 to get actual signal
rec_sig = [rec_sig,[6:149])];
%cyclic prefix removal
Y = [rec_sig,[17:144])];
% freq domain transform
FFT_data =fft(Y.').';
% demodulation
z= comm.QPSKDemodulator('PhaseOffset',pi/16);
demod_data = step(z,FFT_data');
% [no_of_error,ratio] = biterr(t_data,demod_data);
X=zeros(128,128);
for l = 1:1:128
X(l,l)=IFFT_data(l);
end
H = Y*inv(X);
G = H';
H = fft(G)
u=rand(128,128);
F=fft(u)*inv(u);% 'F' is the twiddle factor matrix..
gg=zeros(128,128);
for i=1:128
gg(i,i)=G(i);
end
gg_myu = sum(gg, 1)/128;
gg_mid = gg - gg_myu(ones(128,1),;
sum_gg_mid= sum(gg_mid, 1);
Rgg = (gg_mid' * gg_mid- (sum_gg_mid' * sum_gg_mid) / 128) / (128 - 1);
for m=1:200

for n=1:6.0

SNR_send=5*n;
XFG=X*H;
n1=ones(128,1);
n1=n1*0.000000000000000001i;
noise=awgn(n1,SNR_send);
variance=var(noise);
N=fft(noise);
Y=XFG+N;
%the mean squared error for the LS estimator..
mean_squared_error_ls=LS_MSE_calc(X,H,Y);
%the mean squared error for the MMSE estimator..
mean_squared_error_mmse=MMSE_MSE_calc(X,H,Y,Rgg,variance);
SNR=SNR_send;
mmse_mse(m,n)=mean_squared_error_mmse;
ls_mse(m,n)=mean_squared_error_ls;
end;

end;
ls_mse
mmse_mse
mmse_mse_ave=mean(mmse_mse)
ls_mse_ave=mean(ls_mse)
%plotting
semilogy(SNR,mmse_mse_ave,'k-');
grid on;
xlabel('SNR in DB');
ylabel('mean squared error');
title('PLOT OF SNR V/S MSE FOR AN OFDM SYSTEM WITH MMSE/LS ESTIMATOR BASED RECEIVERS');

hold on;
semilogy(SNR,ls_mse_ave,'b*');
semilogy(SNR,ls_mse_ave,'b-');
semilogy(SNR,mmse_mse_ave,'kv');
grid on;
xlabel('SNR in DB');
ylabel('mean squared error');
title('PLOT OF SNR V/S MSE FOR AN OFDM SYSTEM WITH MMSE/LS ESTIMATOR BASED RECEIVERS');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Function for LS estimation code%%%%%%5

function ms_error=LS_MSE_calc(X,H,Y);
Hls =(inv(X)) * Y;

ms_error_mat=mean(((abs(H-Hls))/abs(H)).^2);
for i=1:128
if(ms_error_mat(i)~=0)
ms_error=ms_error_mat(i);
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%function for mmse code%%%%%%%%%



function ms_error=MMSE_MSE_calc(X,H,Y,Rgg,variance);
%Hmmse=F*Rgg*inv(Rgy)*Y;
u=rand(128,128);
F=fft(u)*inv(u);%The 64 X 64 twiddle factor matrix..
I=eye(128,128);
Rgy=Rgg * F'* X';
Ryy=X * F * Rgg * F' *X' + variance * I;
for i=1:64
yy(i,i)=Y(i);
end
Gmmse=Rgy * inv(Ryy)* Y;
Hmmse=fft(Gmmse);

ms_error_mat=mean(((abs(H)-abs(Hmmse))/abs(H)).^2);
for i=1:128
if(ms_error_mat(i)~=0)
ms_error=ms_error_mat(i);
end
end