Ofdm multipath fading channel output

clc
close all
clear all
nFFTSize = 64;
% for each symbol bits a1 to a52 are assigned to subcarrier 
% index [-26 to -1 1 to 26] 
subcarrierIndex = [11:26];
nBit = 2500; 
input = round(rand(1,nBit));
nBitPerSymbol = length(subcarrierIndex);
nSymbol = ceil(nBit/nBitPerSymbol);
ipMod = 2*input - 1; 
ipMod = [ipMod zeros(1,nBitPerSymbol*nSymbol-nBit)];%zreros(48)
% size(ipMod)
ipMod = reshape(ipMod,nSymbol,nBitPerSymbol);
size(ipMod)
tx1=[];
for ii = 1:nSymbol
inputiFFT = zeros(1,nFFTSize);
% assigning bits to subcarriers
inputiFFT(subcarrierIndex+nFFTSize/2+1) = ipMod(ii,:);
inputiFFT = fftshift(inputiFFT); 
outputiFFT = ifft(inputiFFT,nFFTSize);
% adding cyclic prefix of 16 samples 
outputiFFT_with_CP = [outputiFFT(49:64) outputiFFT];
tx1=[tx1 outputiFFT_with_CP];
end
%% 2nd user
subcarrierIndex = [-26:-1 1:26];
nBit = 2500; 
input = round(rand(1,nBit));
nBitPerSymbol = 52;
nSymbol = ceil(nBit/nBitPerSymbol);
ipMod = 2*input - 1; 
ipMod = [ipMod zeros(1,nBitPerSymbol*nSymbol-nBit)];
ipMod = reshape(ipMod,nSymbol,nBitPerSymbol);
tx2=[];
for ii = 1:nSymbol
inputiFFT = zeros(1,nFFTSize);
% assigning bits a1 to a52 to subcarriers [-20:-10 -5:-1, 1 to 20]
inputiFFT(subcarrierIndex+nFFTSize/2+1) = ipMod(ii,:)
inputiFFT = fftshift(inputiFFT); 
outputiFFT = ifft(inputiFFT,nFFTSize);
% adding cyclic prefix of 16 samples 
outputiFFT_with_CP1 = [outputiFFT(49:64) outputiFFT];
tx2=[tx2 outputiFFT_with_CP1];
end
figure,
fsMHz = 1;
[Pxx1,W] = pwelch(tx1 ,[],[],512,1); 
[Pxx2,W] = pwelch(tx2 ,[],[],512,1); 
pxx=[fftshift(Pxx1);fftshift(Pxx2)];
% plot([-2048:2047]*fsMHz/4096,10*log10(pxx));
figure(1);
plot(10*log10(pxx));
xlabel('frequency, MHz')
ylabel('power spectral density');
% % % channel models
for i=1:5
ch(i) = exp(-i/5);
end
for j=1:5
ch_norm(j) = ch(j)/norm(ch);
ch_norm(j) = ch_norm(j)*(randn(1)-i*randn(1))
end
%%% fade variable
N=10000;%Number of correlated Rayleigh samples to be generated
fd=100;%doppler frequency
fs=10000;%sampling frequency
fm=fd/fs;%normalized doppler frequency
km=floor(fm*N);%A point nerear to doppler frequency
%filter F_M(k) to shape gaussian random process
F_M(1) = 0;
for k=2:km
F_M(k)= sqrt(1/(2*(sqrt(1-((k-1)/(N*fm)).^2))));
end
F_M(km+1)= sqrt((km/2)*((pi/2) - atan((km-1)/(sqrt((2*km)-1)))));
for k=km+2:N-km
F_M(k)=0;
end
F_M(N-km+1)=sqrt((km/2)*((pi/2) - atan((km-1)/(sqrt((2*km)-1)))));
for k= N-km+2:N
F_M(k)= sqrt(1/(2*(sqrt(1-((N-k+1)/(N*fm)).^2))));
end

for k=1:100
F_M(k)=F_M(k)*(randn(1)-i*randn(1));
end
f_m=ifft(F_M,N);
abs_f_m=abs(f_m);
% f_m_db=20*log(abs_f_m);
figure(1)
hist(abs_f_m,100);
figure(2);
% hist(real(f_m),100);
figure(3)
% plot(1:10000,f_m_db);

%%% fading channel=time invariant channel*fade variable

for j=1:5
for k=1:N
r_channel(k,j)=real(ch_norm(j))*real(f_m(k));
i_channel(k,j)=imag(ch_norm(j))*imag(f_m(k));
end
end

fade_channel=r_channel+i*i_channel;
figure(4)
plot(1:10000,fade_channel);
y=abs(r_channel).^2
z=abs(i_channel).^2
mag=sqrt(y+z);
var=mean(mean(mag));
sigma=sqrt(var);
v=sigma*(randn(1)+i*randn(1));
for j=1:49
sigma=sqrt(var);
for k=1:68
v(j,k)=sigma*(randn(1)+i*randn(1));
end
% y1(j,1:68)=cconv(fade_channel(j,1:5),outputiFFT_wi th_CP(j,1:68))+v(j,1:68);
% y2(j,1:68)=cconv(fade_channel(j,1:5),outputiFFT_wi th_CP1(j,1:68))+v(j,1:68);
y1(j,1:68)=cconv(fade_channel(1,1:5),outputiFFT_wi th_CP(j,1:68));
y2(j,1:68)=cconv(fade_channel(1,1:5),outputiFFT_wi th_CP1(j,1:68));
% y1(j,k)=cconv(outputiFFT_with_CP(j,k),v(j,k));
% y2(j,k)=cconv(outputiFFT_with_CP1(j,k),v(j,k));
end
fsMHz = 1;
[Pxx1,W] = pwelch(y1 ,[],[],512,1); 
[Pxx2,W] = pwelch(y2 ,[],[],512,1); 
pxx=[fftshift(Pxx1);fftshift(Pxx2)];
% plot([-2048:2047]*fsMHz/4096,10*log10(pxx));
figure;
plot(10*log10(pxx));
xlabel('frequency, MHz')
ylabel('power spectral density');
title('noise added signal');
% TX=[tx1;tx2]; 
% plot(y1);

%% channel estimation in receiver
% for snr=1:5:40
% for j=1:
% % var=10^(-snr/10);
% 
% sigma=sqrt(var);
% for k=1:80
% v(j,k)=sigma*(randn(1)+i*randn(1));
% end
% y(j,1:80)=cconv(fade_channel(j,1:5),tx2(j,1:80),80 )+v(j,1:80);
% end
% end
% plot(y)