%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                             %%%
%%%     The CS-CDMA system                                      %%%
%%%                                                             %%%
%%%     by: Anvar M. Narzullaev                                 %%%
%%%     Mobile Communication Lab., Yeungnam University          %%%
%%%     copyright, 2006                                         %%%
%%%                                                             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;
%====================================%
% S Y S T E M   P A R A M E T E R S  %
%====================================%
cs_bits = 3;        % Number of Code selecting bits
para = cs_bits+1;
ml = 1;             % Modulation level : BPSK
sr = 33333.33333;   % symbol rate
br   = sr * ml;     % bit rate
noe_total = 0;    % Number of error data
nod_total = 0;    % Number of transmitted data
noe2 = zeros(1,para);
noe = zeros(1,para);
nod = zeros(1,para);
code_table=((dec2bin((0:(2^cs_bits)-1)')*2)-97);

%=====================================%
%%====== SIMULATION PARAMETERS ======%%
%=====================================%
nd   = 100;          % number of symbol
nloop = 5000;  % Number of simulation loops
use_pn = 1; % choose between PN and Walsh codes; 0 is PN 1 is Walsh
use_casuh = 1;

if use_pn == 1 && use_casuh == 1
    fprintf('\nCasuh code simulation\n');
    sf = 16;
    pn_code1 = pnsequence([1 0 0 1 1], [1 1 1 1], sf-1)*2-1;
    set = makeset(pn_code1, 2^cs_bits, sf-1);
    spr_code = [set, -ones(2^cs_bits,1)];
end

if use_pn == 1 && use_casuh == 0
    fprintf('\nPN code simulation\n');
    sf = 15;
    pn_code1 = pnsequence([1 0 0 1 1], [1 1 1 1], sf)*2-1;
    spr_code = makeset(pn_code1, 2^cs_bits, sf);    
end

if use_pn ~= 1
    fprintf('\nWalsh code simulation\n');
    sf = 16;
    walsh = hadamard(sf);   % Spreading sequence (walsh)
    set = walsh(1:8,:);
    spr_code = set;
end
%=====================================%
%== end of SIMULATION PARAMETERS   ==%%
%=====================================%


for ebn0 = 0:4:40
    %************************** Fading initialization **************************
for rfade = [0 1]
for now1 = [2 6]
    %rfade  = 0;                                         % Rayleigh fading 0:nothing 1:consider
    %now1   = 2;                                         % number of directwave + delayed wave (2 or 6)
    tstp   = 1 / sr / sf;                               % time resolution
    fd     = 200;                                       % doppler frequency [Hz]
    flat   = 0;                                         % flat Rayleigh environment
    itndel = nd * sf * 30;                              % number of fading counter to skip
    if now1 == 2 %2 path rayleigh fading
        itau   = [0, 2];                                % delay time
        dlvl1  = [0.0, 25.0];                           % attenuation level
        n0     = [6, 7];                                % number of waves to generate fading
        th1    = [0.0, 0.0];                            % initial Phase of delayed wave
        itnd1  = [3001, 4004];                          % set fading counter
    elseif now1 == 6 % 6 path rayleigh fading
        itau   = [0, 1, 2, 3, 4, 5];                   % delay time
        dlvl1  = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0];       % attenuation level
        n0     = [6, 7, 6, 7, 6, 7];                   % number of waves to generate fading
        th1    = [0.0, 0.2, 0.4, 0.6, 0.8, 0.9];       % initial Phase of delayed wave
        itnd1  = [1000, 2000, 3001, 4004, 5005, 6006]; % set fading counter
    end
    
    %**************************** START CALCULATION ****************************
    
    for iii = 1:nloop 
    %**************************** TRANSMITTER ****************************
        seridata = rand(1,para*nd*ml)>0.5;          % Serial data
        seridata2 = seridata*2-1;                   % BPSK
        paradata = reshape(seridata2, para, nd*ml); % S-to-P convert
        sp_bit(1,:) = paradata(end,:);              % bits for spreading
       
    
        %%===============================================%%
        %%== Spreading Code Selection proccess           %%
        %%===============================================%%
        selecting_bits = paradata(1:cs_bits,:)';
        n=1;
        for x = 1:nd*ml
            selected(1,n) = code_select(selecting_bits(x,:), code_table, cs_bits);
            n=n+1;
        end

        m = 1;
        for x = 1:nd*ml
            code(m,:) = spr_code(selected(x),:);
            m=m+1;
        end
        %%***********************************************%%
        %%== END of Spreading Code Selection proccess    %%
        %%***********************************************%%
           
    
        %%===============================================%%
        %%== Spreading                                   %%
        %%===============================================%%
        n=1;
        for x=1:nd*ml
            spread_data(n:n+sf-1) = spread(sp_bit(x),code(x,:))';
            n=n+sf;
        end
        %%***********************************************%%
        %%== END of Spreading                            %%
        %%***********************************************%%
    
    
    
        %%===============================================%%
        %%== IFFT                                        %%
        %%===============================================%%
        %x = ifft(spread_data);
        %ich = real(x);
        %qch = imag(x);
        ich = spread_data;
        qch = zeros(size(spread_data));
        %%***********************************************%%
        %%== END of IFFT                                 %%
        %%***********************************************%%     

        
        %***************************** Fading channel ******************************
 
        if rfade == 0
            ich1 = ich;
            qch1 = qch;
        else
            [ich1,qch1] = sefade(ich,qch,itau,dlvl1,th1,n0,itnd1, ... % fading channel
                                 now1,length(ich),tstp,fd,flat);
            itnd1 = itnd1 + itndel;
        end
    
            
        
    %**************************** RECEIVER ****************************    
        
        %%================ AWGN
        spow = sum(rot90(ich.^2 + qch.^2)) / nd; % attenuation Calculation
        attn = sqrt(0.5 * spow * sr / br * 10^(-ebn0/10)); %attenuation calculation
        [ich1, qch1]=comb(ich1,qch1,attn);
        %%================ End of AWGN
            
        %%================ FFT
        %clear i;
        %rx = ich1 + qch1.*i;
        %ry = fft(rx);
        %ich2=real(ry);
        ich2 = ich1;
        %%================ End of FFT 
    
    
        %%===============================================%%
        %%== Desicion stage (correlation)                %%
        %%===============================================%%
        temp = zeros(1,nd*ml);
        m = 1;
        for x = 1:nd*ml
            s1 = 0;
            for n = 1:2^cs_bits
                s2 = abs(sum(ich2(m:m+sf-1).*spr_code(n,:)));
                if s2>s1
                    s1 = s2;
                    temp(x) = n;
                end
            end
            m = m+sf;
        end

        m = 1;
        for x = 1:nd*ml
            code2(m,:) = spr_code(temp(x),:);
            m=m+1;
        end       
        %%===============================================%%
        %%== END of Desicion stage (correlation)         %%
        %%===============================================%%
    
    
    
        %%===============================================%%
        %%== DEspreading                                 %%
        %%===============================================%%
        n=1;
        for x=1:nd*ml
            despread_data(x) = despread(ich2(n:n+sf-1),code2(x,:))';
            n=n+sf;
        end
        %%***********************************************%%
        %%== END of DEspreading                          %%
        %%***********************************************%%
        orig_data = zeros(cs_bits,nd*ml);
        for x = 1:nd*ml
            orig_data(:,x) = code_table(temp(x),:)';        
        end
    
        recovered_data = [orig_data; despread_data];
        seridata3 = reshape(recovered_data,1,nd*ml*para);
        demodata = seridata3>0;
        
        recovered_data = (recovered_data>0)*2-1;
  
        %========= BER calc. ===========%
        for m = 1:para
            for n = 1:length(paradata)
                if recovered_data(m,n) ~= paradata(m,n)
                    noe(m) = noe(m) + 1;
                end
            end
        end
        
        for n = 1:length(seridata)
            if seridata(n) ~= demodata(n)
                noe_total = noe_total + 1;
            end
        end
        %========= end of BER calc. ===========%
        
    end %% END for iii=nloop
    %ebn0;
    %noe;
    %noe_total;
    ber = noe./(nloop*length(paradata));
    ber_total = noe_total/(nloop*length(seridata));
    
    %********************** Output result ***************************
    fprintf('%d\t%e\t%e\t',ebn0,ber,ber_total);
    fprintf('\n');
    
    if rfade == 0
        fid = fopen('cs_cdma_awgn.dat','a');
        fprintf(fid,'%d\t%e\t%e\t',ebn0,ber,ber_total); % fprintf: built in function
        fprintf(fid,'\n'); % fprintf: built in function
        fclose(fid);
    end
    if rfade == 1 && now1 == 2
        fid = fopen('cs_cdma_1path.dat','a');
        fprintf(fid,'%d\t%e\t%e\t',ebn0,ber,ber_total); % fprintf: built in function
        fprintf(fid,'\n'); % fprintf: built in function
        fclose(fid);
    end
    if rfade == 1 && now1 == 6 
        fid = fopen('cs_cdma_6path.dat','a');
        fprintf(fid,'%d\t%e\t%e\t',ebn0,ber,ber_total); % fprintf: built in function
        fprintf(fid,'\n'); % fprintf: built in function
        fclose(fid);
    end
    
    
    noe = zeros(1,para);
    noe_total = 0;
end %% for now1 = [2 6]
end %% for rfade = [0 1]
end %% END for ebn0=0:4:20