close all
clear all

%-----------------------------------------------------------------
%Constants
%-----------------------------------------------------------------
%def_f and B depends on each other
%Subcarrier scacing (valid values are 15, 30, 60, 120)
del_f     = 15; %kHz        
B         = 5;  %MHz      %Channel bandwidth
M_mod     = 16;  %modulační řád - 16qam
k         = log2(M_mod);  %počet bitů na symbol
SNR       = 0:1:50;
SNR2      = 30;
BER       = zeros(1, numel(SNR));
num_err   = zeros(1, numel(SNR));
N_symb    = 4096;
%Carrier frequncy for converion to band signal
fc = 10e6;


%-----------------------------------------------------------------
%Configuration flags for simulation
%-----------------------------------------------------------------
%Draw or not graphs - should be 0 when computing FFT quantization error
sim_graph_draw = 1;
%Select between matlab FFT function and Xilinx MEX FFT function
xfft_mex_use = 0;

%--------------Data format----------------------
%Select between floating point and fixed point data format
fxpt_use = 0;
%If fixed point, select signess
fxpt_signed = 1;
%Select data format width in bits
data_width = 16;
%Range of number data bits, minumum 16, maximum 32
%If both same, 
input_width_start = 16;
input_data_stop = 16;

%-----------------------------------------------------------------
%SIMULATION
%-----------------------------------------------------------------

for data_width = input_width_start:1:input_data_stop
	
end %input_width = input_width_start:1:input_data_stop

%For data generation
%RB_num = 275;   %Number of Resource blocks in range 27-275
DMRS_conf = 1;  %DRMS symbols configuration

RB_subcar = 12;   %12 subcariers in one Resource Block
%OFDM symbols in one Resource Block
if del_f == 15
  RB_symb   = 14;
elseif del_f == 30
  RB_symb   = 28;
elseif del_f == 60
  RB_symb   = 56;
elseif del_f == 120
  RB_symb   = 128;
elseif del_f == 240
  RB_symb   = 256;
end
%FFT_data = 3300;  %Data subcariers in FFT size
%FFT_symb = 7;     %Data symbols in FFT size

%-----------------------------------------------------------------
%Calculation of parameters
%-----------------------------------------------------------------
[M_data, N_fft, cp_1st, cp_oth, RB_num, Ts] = f_calc_params(del_f, B)
N_fft = 4096;
%N = M_mod*M_data; %for 2 bits per symbol (qpsk)
%t=-0.5e3:Ts:0.5e3;
%-----------------------------------------------------------------
%Generating data and DMRS, mapping DMRS and data into RB, mapping 
%RBs into matrix ready for zeros subars and IFFT and so on
%-----------------------------------------------------------------
if DMRS_conf==1
  %----%Data and DMRS generation for given DMRS configuration--------------
  RB_dmrs_subcar = 12;   %Number of dmrs subcariers in Besource Block
  RB_dmrs_symb = 1;     %Number of mdrs symbols in Resource Block
  RB_data_subcar = 12;  %Number of data subcariers in Resource Block
  RB_data_symb = RB_symb-RB_dmrs_symb;     %Number of data symbols in Resource Block

  
  max_data_size = RB_data_subcar*(RB_data_symb)*RB_num;
  max_dmrs_size = RB_dmrs_subcar*RB_dmrs_symb*RB_num;
  
  %Data generation
	data_in = randi([0 M_mod-1], max_data_size, 1);

  %Data Modulation
	if fxpt_use == 1 
% 		data_mod=qammod(data_in, M_mod, 'OutputDataType', ...
% 																				numerictype(fxpt_signed,input_width))./k;
																			
 		data_mod = (1/k).*qammod(data_in, M_mod);
%		data_mod2 = (1/4).*qammod(data_in, M_mod);
	else
		data_mod=qammod(data_in, M_mod);
	end
%  data_mod=qammod(data_in, M_mod);

  if sim_graph_draw == 1
		figure;
		subplot(3,3,1)
		plot(data_mod, '.', 'MarkerSize', 10)
		title('I/Q diagram 16QAM modulovaných dat', 'FontSize', 14)
		xlim([-4 4]);
		ylim([-4 4]);
		xlabel('I', 'FontSize', 13)
		ylabel('Q', 'FontSize', 13)
	end
  
  %DMRS generation
	val = 1/sqrt(2);
  dmrs_symb = randsrc(max_dmrs_size, 1, [val+val*j val-val*j -val+val*j -val-val*j]);
  
  %----%Resource Block filling by data and DMRS----------------------------  
  %Reshaping modulated data to (RB_data_subcar*RB_num) x (RB_data_symb) matrix
  data_mod_resh = reshape(data_mod, RB_data_subcar*RB_num, RB_data_symb);
  
  %Concatenating DMRS symbols and modulated data
  data_dmrs = [dmrs_symb,data_mod_resh];
  
end

%-----------------------------------------------------------------
%Inserting zeros 
%-----------------------------------------------------------------
data_zeros = [data_dmrs(1:M_data/2,:); ... 
              zeros(N_fft-M_data,RB_symb); ... 
              data_dmrs(M_data/2+1:end,:)   ];

%-----------------------------------------------------------------
%IFFT
%-----------------------------------------------------------------
%[M_ifft] = f_calc_ifft(data_mod, M_data, N_fft, N_fft_symb);

if xfft_mex_use == 1
  %For MEX simulation
	%Pre-allocating
	M_ifft = data_zeros;

	for m=1:1:numel(data_zeros(1,:))
		%Make IFFT over all subcariers and one symbol
		M_ifft(:, m) = xfft_mex(data_zeros(:, m), 0, 4096, fxpt_use, data_width);
	end
else
	%For clasic simulation
	M_ifft = ifft(data_zeros);
end

data_ifft_resh = reshape(M_ifft, 1, numel(M_ifft));

if sim_graph_draw == 1
	subplot(3,3,2)
	pwelch(data_ifft_resh, [], [], [], del_f*N_fft, 'centered');
	title('Spektrální hustota výkonu signálu před přenosovým kanálem', 'FontSize', 14)
	ylabel('Spektrální hustota výkonu [dB/Hz]', 'FontSize', 13)
	xlabel('Frekvence [MHz]', 'FontSize', 13)
	%pwelch(M_ifft);
	%title('Spectral estimate')
end

%-----------------------------------------------------------------
%Insert of Cyclic Prefix
%-----------------------------------------------------------------
[data_cp] = f_insert_cp(data_ifft_resh, N_fft, cp_1st, cp_oth, del_f);

%-----------------------------------------------------------------
%Convert to band signal
%-----------------------------------------------------------------
%data_band = real(data_cp).*cos(2*pi*fc*t)- ...
						%imag(data_cp).*sin(2*pi*fc*t);

%figure(2)
%plot(t,data_band)
					
%-----------------------------------------------------------------
%ISI - Multipath Propagation
%-----------------------------------------------------------------
% for i=1:3
%   i
%  if i==1
      h = [1 zeros(1,31) 0.5 zeros(1,123) 0.1];
      %h = 1;
%  elseif i==2
%    h = [1 zeros(1,200) 0.5 zeros(1,600) 0.1];
%  else
%    h = 1;
%  end
data_isi=conv(data_cp, h);
data_isi=data_isi(1, 1:end-length(h)+1);

%figure();
%pwelch(data_cp, [], [], [], 1/del_f*1000, 'centered');
%title('Spectral estimate with CP (centered)')

%-----------------------------------------------------------------
%Transmission
%-----------------------------------------------------------------
%chahnge M_ifft with data_cp after implementation of CP remove 
% for j=1:1:numel(SNR)
%   j
%rx_sig = awgn(data_isi, SNR(j), 'measured');
rx_sig = awgn(data_isi, SNR2, 'measured');

%-----------------------------------------------------------------
%Convert to complex signal
%-----------------------------------------------------------------
%data_hilbert = hilbert(rx_sig);

%data_cmplx = data_hilbert.*exp(-1i*2*pi*fc)


%-----------------------------------------------------------------
%Remove of Cycle Prefix
%-----------------------------------------------------------------
[data_cp_rm] = f_remove_cp2(rx_sig, RB_symb, N_fft, cp_1st, cp_oth, del_f);
%num_err_cp = data_cp_rm - M_ifft;

if sim_graph_draw == 1
	subplot(3,3,3)
	pwelch(data_cp_rm, [], [], [], del_f*N_fft, 'centered');
	title('Spektrální hustota výkonu signálu po průchodu přenosovým kanálem', 'FontSize', 14)
	ylabel('Spektrální hustota výkonu [dB/Hz]', 'FontSize', 13)
	xlabel('Frekvence [MHz]', 'FontSize', 13)
end

%-----------------------------------------------------------------
%FFT
%-----------------------------------------------------------------
%[M_fft] = f_calc_fft(rx_sig, M_data, M, N_fft_symb);
%[M_fft] = f_calc_fft(data_cp_rm, M_data, N_fft, N_fft_symb);

data_cp_rm_resh = reshape(data_cp_rm, N_fft, numel(data_cp_rm)/N_fft);

if xfft_mex_use == 1
  %For MEX simulation
	%Pre-allocation
	M_fft = data_cp_rm_resh;

	for m=1:1:numel(data_zeros(1,:))
		%Make FFT over all subcariers and one symbol
		M_fft(:, m) = xfft_mex(data_cp_rm_resh(:, m), 1, 4096, fxpt_use, data_width);
	end
else
	%For clasic simulation
	M_fft = fft(data_cp_rm_resh);
end

data_fft_resh = reshape(M_fft, ...
                        1, numel(M_fft));

if sim_graph_draw == 1
	subplot(3,3,4)
	plot(data_fft_resh, '.', 'MarkerSize', 10)
	title('IQ diagram modulovaných dat po průchodu přenosovým kanálem', 'FontSize', 14)
	xlim([-4 4]);
	ylim([-4 4]);
	xlabel('I', 'FontSize', 13)
	ylabel('Q', 'FontSize', 13)
end

%-----------------------------------------------------------------
%Removing zeros
%-----------------------------------------------------------------
data_fft_nozeros  = [M_fft(1:M_data/2,:); ...
                    M_fft(M_data/2+N_fft-M_data+1:end,:) ];

%-----------------------------------------------------------------
%Equalization, demapping and demodulation
%-----------------------------------------------------------------
data_idx = 1;   
RB_beg_idx = 1;

%data_eq = zeros(RB_data_subcar, numel(data_fft_nozeros)/RB_data_subcar)

if DMRS_conf == 1
  %prenosova funkce kanalu
  H = data_fft_nozeros(:, 1)./dmrs_symb(:, 1);
	
	if sim_graph_draw == 1
		subplot(3,3,5)
		%plot(abs(D_graphs./data_dmrs_graphs))
		plot([-numel(H)/2:numel(H)/2-1],[abs(H(numel(H)/2:end,:)); abs(H(1:numel(H)/2-1, :))])
		title('Modulová část přenosové funkce kanálu po průchodu kanálem', 'FontSize', 14)
		xlabel('Datová subnosná m [-]', 'FontSize', 13)
		ylabel('H(m)', 'FontSize', 13)
		xlim([-numel(H)/2 numel(H)/2-1])
	end
		
  
  %equalization function 
  Z = 1./H;
	
	if sim_graph_draw == 1
		subplot(3,3,6)
		plot([-numel(Z)/2:numel(Z)/2-1], [abs(Z(numel(Z)/2:end,:)); abs(Z(1:numel(Z)/2-1, :))]);
		title('Modulová část vyrovnávací funkce kanálu', 'FontSize', 14)
		xlabel('Datová subnosná m [-]', 'FontSize', 13)
		ylabel('1/H(m)', 'FontSize', 13)
		xlim([-numel(Z)/2 numel(Z)/2-1])
	end
		
  %matrix for equalization
  U = repmat(Z, 1, RB_symb-RB_dmrs_symb);
  %Data equalization
  data_eq = data_fft_nozeros(:, 2:end).*U;
  
  %reshaping data_mod to numel(data_eq) rows and 1 column
  data_for_demod = reshape(data_eq, numel(data_eq), 1);
  
	if sim_graph_draw == 1
		subplot(3,3,8)
		plot(data_for_demod, '.', 'LineWidth', 4)
		title('IQ diagram opravených modulovaných dat', 'FontSize', 14)
		xlim([-4 4]);
		ylim([-4 4]);
		xlabel('I', 'FontSize', 13)
		ylabel('Q', 'FontSize', 13)
	end
  
end

%-----------------------------------------------------------------
%Demodulation
%-----------------------------------------------------------------
%data_out = qamdemod(data_for_demod, M_mod)*k;

if fxpt_use == 1 
% 		data_mod=qammod(data_in, M_mod, 'OutputDataType', ...
% 																				numerictype(fxpt_signed,input_width))./k;
																			
	data_out = qamdemod(data_for_demod.*k, M_mod);
else
	data_out = qamdemod(data_for_demod, M_mod);
end


%comapare of input and output signal of modulator
[num_err, BER] = biterr(data_in, data_out);
%[num_err(i,j), BER(i,j)] = biterr(data_in, data_out);

% end
% end
% 
% figure;
% semilogy(SNR, BER(1, :), '-b+', 'LineWidth',2);
% hold on;
% semilogy(SNR, BER(2, :), '-go', 'LineWidth',2);
% hold on;
% semilogy(SNR, BER(3, :), '-r.', 'LineWidth',2);
% title('Závislost bitové chybovosti na SNR pro 16QAM, FFT 4096', 'FontSize', 14)
% xlabel('SNR [dB]', 'FontSize', 13);
% ylabel('BER [-]', 'FontSize', 13);
% legend('Tmax<TCP', 'Tmax>TCP', 'AWGN');
% grid on;






