%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% 	Performance of Hybrid FH/DS CDMA Systems with Tranmitter Diversity
%	
% 	Copyright(c) 2009-2015 HanYang University Wireless Communication Lab
% 
% 	Date		Name            Version         Description
%   ----------  -----------     -----------     ---------------------------
% 	2009.05.22	JINMINGSONG     v1.00           Create
% 	2009.05.22	JINMINGSONG     v1.01           Upgrade
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

clc, clear;

%----------------------------- Initialization -----------------------------
format long;
sam = 1;
para=256;   % Number of parallel channel to transmit (points)
fftlen=256*sam; % FFT length
noc=256;    % Number of carrier
nd=1;       % Number of information OFDM symbol for one loop
ml=4;       % Modulation level : QPSK
ml_q=4;       % Modulation level : 16 QAM
sr=250000;  % Symbol rate
br=sr.*ml;  % Bit rate per carrier
br_q = sr.*ml_q;
gilen=64*sam;   % Length of guard interval (points)
berg = [];
ber_dg_QPSK = [];
ber_dg_QAM = [];
ber_cg1 = [];
ber_crg = [];
CR = 1.2;
PAR = 2;
countn = 0; % count the number of peak points
m2 = mseq(9, [4,9], [1,1,1,1,1,1,1,1,1]);           % WCL_mseq.m
pad = zeros(511,1);
mse = [m2 pad];
mse = [mse mse];
codesidelen = 0; % 4,8,16
codesideinf = [];
L = 4;
pw = 2;
qq = 1; % number of references
sideinf = (hadamard(32)+ones(32,32))/2;
codesideinf = sideinf;
mse = mse(1:32,1:(1024-64-codesidelen));
% Q matrix
Q = [];
Qs = [];
for ii = 1:16
    pr = (ii-1)*16+1;
    for jj = 1:L
        l = jj-1;
        Q(ii,jj) = exp(-j*2*pi*pr*l/256);
    end
end
Q_h = Q';

for ii = 1:256
    s = ii;
    for jj = 1:L
        l = jj-1;
        Qs(ii,jj) = exp(-j*2*pi*s*l/256);
    end
end

fftlen1 = fftlen+gilen;

%******************* Fading initialization ********************
% If you use fading function "sefade", you can initialize all of parameters.
% Otherwise you can comment out the following initialization.
% The detailed explanation of all of valiables are mentioned in Program 2-8.

% Time resolution

tstp=1/sr/(fftlen+gilen); 

% Arrival time for each multipath normalized by tstp
% If you would like to simulate under one path fading model, you have only to set 
% direct wave.

itau = [0,2,3,4];

% Mean power for each multipath normalized by direct wave.
% If you would like to simulate under one path fading model, you have only to set 
% direct wave.
dlvl = [0,10,20,25];

% Number of waves to generate fading for each multipath.
% In normal case, more than six waves are needed to generate Rayleigh fading
n0=[6,7,6,7];

% Initial Phase of delayed wave
% In this simulation four-path Rayleigh fading are considered.
th1=[0.0,0.0,0.0,0.0];

% Number of fading counter to skip 
itnd0=nd*(fftlen+gilen)*10;

% Initial value of fading counter
% In this simulation one-path Rayleigh fading are considered.
% Therefore one fading counter are needed.
  
itnd1=[1000,2000,3000,4000];

% Number of directwave + Number of delayed wave
% In this simulation one-path Rayleigh fading are considered
now1=4;        

% Maximum Doppler frequency [Hz]
% You can insert your favorite value
fd=200;       

% You can decide two mode to simulate fading by changing the variable flat
% flat     : flat fading or not 
% (1->flat (only amplitude is fluctuated),0->nomal(phase and amplitude are fluctutated)
flat = 1;
%*************************** reference signal*********************
for ebn0 = 0:5:20
%     ebn0 = esn0-3; % SNR is Es/N0; QPSK
%    Eb/N0
%************************** main loop part **************************

nloop=100;  % Number of simulation loops
nloop1 = nloop;
noe = 0;    % Number of error data
nod = 0;    % Number of transmitted data
eop=0;      % Number of error packet
nop=0;    % Number of transmitted packet

noe_c = 0;
nod_c = 0;
eop_c = 0;
nop_c = 0;

for iii=1:nloop1
k = zeros(para*sam,1);
kadd = zeros(para*sam,1);
pilot_rec_c = [];
pilot_rec_cr = [];
ebn0
iii
%************************** transmitter *********************************
%************************** Data generation **************************** 
hamweight=0.5;
seldata=rand(1,para*nd*ml-64-codesidelen)>hamweight;  %  rand : built in function
pilot = rand(1,64)>hamweight;
%%%%%%%%%%%%%%% Select M-sequence %%%%%%%%%%%%%%%%%%%
clip_ratio = [];
peak_ratio = [];

scrdata = xor(seldata,mse(1,:));
sidedata = scrdata;
pilotdata = [];
for dd = 1:64
    pilotdata = [pilotdata pilot(dd) sidedata(1,1+(dd-1)*15:15+(dd-1)*15)];
end
    paradata = reshape(pilotdata,para,nd*ml);
    [ich1 qch1] = WCL_qammod(paradata,para,nd,ml);              % WCL_qammod.m
x=ich1+qch1.*j;
for ii = 1:16
    x((ii-1)*16+1) = pw*x((ii-1)*16+1);
end
pilotsig = [];
for pp = 1:16
    pilotsig = [pilotsig x(1+16*(pp-1))];
end
y=ifft(x);      %  ifft : built in function
ich2 = real(y);
qch2 = imag(y);
rms = sqrt(1/(para*sam)*sum((abs(y).^2)));
limit = CR*rms;%(10.^(CR/20))*rms;

ich2_c = real(y);
qch2_c = imag(y);  

%********************* Filtering ****************************
ich3_c=ich2_c';
qch3_c=qch2_c';

%*******************  Guard Interval Insertion  ********************
fftlen2 = fftlen+gilen;
[ich3_c,qch3_c] = giins(ich3_c,qch3_c,fftlen,gilen,nd);     % WCL_giins.m

%********* Attenuation Calculation *********
spow_c=sum(ich3_c.^2+qch3_c.^2)/nd./para;  %  sum : built in function
attn_c=0.5*spow_c*1/4*10.^(-ebn0/10);
attn_c=sqrt(attn_c);

%***************************  Receiver  *****************************
%********************** Fading channel **********************

% Generated data are fed into a fading simulator
[ifade_c,qfade_c]=sefade(ich3_c,qch3_c,itau,dlvl,th1,n0,itnd1,now1,length(ich3_c),tstp,fd,flat);
% seffade.m
ich3_c = ifade_c;
qch3_c = qfade_c;
% Updata fading counter
itnd1 = itnd1+ itnd0;
% %***************** AWGN addition ********* 

% 4 is only AWGN signal
[ich4_c,qch4_c]=comb(ich3_c,qch3_c,attn_c);     % comb.m

% ****************** Guard interval removal *********
[ich4_c,qch4_c] = girem(ich4_c,qch4_c,fftlen2,gilen,nd);    % girem.m
ich4_c=ich4_c';
qch4_c=qch4_c';
ich5_c=ich4_c';
qch5_c=qch4_c';


for iiii = 1:1
%******************  FFT  ******************
rx_c=ich5_c+qch5_c.*j;
ry_c=fft(rx_c);   	% fft : built in function

ich6_c=real(ry_c);	% real : built in function
qch6_c=imag(ry_c);	% imag : built in function
ich7_c=ich6_c;
qch7_c=qch6_c;
%***************** demodulation *******************
%%%%%%%%%%%%%%%%%% iterative cancellation %%%%%%%%%%%%%%%%%%%%%%%%%%
irp_c = reshape(ich7_c,16,16);
qrp_c = reshape(qch7_c,16,16);

for i = 1:16
    i_rec_c(i) = irp_c(1,i);
    q_rec_c(i) = qrp_c(1,i);
end
    pilot_recsig_c = i_rec_c+j.*q_rec_c;   % received pilot signal (in frequency domain)
    channel_c_f = pilot_recsig_c./pilotsig;
    channel_c_t = inv(Q_h*Q)*Q_h*(real(channel_c_f)'+j*imag(channel_c_f)');

    echannel_c_f = Qs*channel_c_t;
    echannel_c_f = real(echannel_c_f)'+j*imag(echannel_c_f)';
    echannel_c_t = ifft(echannel_c_f);

ich8_c = real((ich7_c'+j*qch7_c')./echannel_c_f);
qch8_c = imag((ich7_c'+j*qch7_c')./echannel_c_f);
ich7_c = ich8_c;
qch7_c = qch8_c;
ich7_c = ich7_c';
qch7_c = qch7_c';

% REMOVE PILOT
ich8_c = reshape(ich7_c,16,16);
qch8_c = reshape(qch7_c,16,16);
ich9_c = [];
qch9_c = [];
for jj = 1:16
  for ii = 2:16
    ich9_c = [ich9_c ich8_c(ii,jj)];
  end
end
for jj = 1:16
  for ii = 2:16
    qch9_c = [qch9_c qch8_c(ii,jj)];
  end
end   
    
[demodata_c]=qamdemod(ich9_c,qch9_c,1,(para-16)*nd,ml);     % WCL_qamdemod.m
demodata_c = (reshape(demodata_c,4,para-16))';
end
%**************  Parallel to serial conversion  *****************
demodata1_c=reshape(demodata_c,1,para*nd*ml-64);
sd = [];
dmse = [];
rcodesideinf = demodata1_c(1:codesidelen);
rdemodata1_c = demodata1_c(codesidelen+1:para*nd*ml-64);
dmse = mse(1,:);
demodata1_c = xor(rdemodata1_c,dmse);
%************************** Bit Error Rate (BER) ****************************

% instantaneous number of error and data  %  sum : built in function
noe2_c = sum(abs(demodata1_c-seldata));
nod2_c=length(seldata);  %  length : built in function

% cumulative the number of error and data in noe and nod
noe_c = noe_c+noe2_c;
nod_c=nod_c+nod2_c;

% calculating PER
  
if noe2_c~=0  
   eop_c=eop_c+1;
else
   eop_c=eop_c;
end   
   eop_c;
   nop_c=nop_c+1;

end %  loop finish 

%------------------------------- Figure Out -------------------------------
per_c = eop_c/nop_c;
ber_c = noe_c/nod_c;
ber_cg1 = [ber_cg1 ber_c];
nloop = nloop1;
end
% EbN0 loop finish

figure;
semilogy(5:5:25, ber_cg1,'r*-');
axis([5 25 10^-6 0.5])
legend(['DS/FH Only'])
xlabel('Eb/No(dB)');
ylabel('Symbol Error Rate');
grid;