function out=mcds_cdma_BER()

%%%*****************************************************************************************
%%% Multiuser coherent MC/DS-CDMA system
%%% BER simulation with frequency interleaving, no time interleaving and coding 
%%% k : bit number in each signal
%%% Ps : average power of the transmitted signal (have been normalized to be 1, i.e.Ps==1)
%%% R : channel coding rate (if no channel coding, R=1)
%%% K : number of transmitted signals in each subchannel
%%% Lb : bit number in each subchannel
%%% N : number of subcarriers
%%% Lp : number of prefix points
%%% Ls : number of points in each OFDM signal,Ls=Lp+N
%%% P : number of spreading code chips (SF)
%%% M : number of data branches (SF=M*P)
%%% nloop : number of simulation loops
%%% Ku : number of total users in a cell ( It should be an even number when channel estimation is used)
%%% scid_u : spread code (the walsh sequence) index of the user
%%% scid_i : spread code (the walsh sequence) index vector of the interfer users, 1*(Ku-1)
%%% Lc : symbol number in each subchannel
%%%******************************************************************************************

Eb_N0=26;
k=6;
Ps=1;
R=1;
K=2;
Lb=K*k;
N=512;
Lp=128;
Ls=N+Lp;
P=32;
nloop=1;

scid_u=2;
scid_i=3:32;
Ku=length(scid_i)+1;
%****************** generate the user data *******************************
%%%******* generate spreading code *************************
scode=Hadamard(P);
%%%*********************************************************
%%%****** generate the information bit sequence ************
ss_code=scode(scid_u,:);                           %%% user spreading code,row vector            
sm=[];
for i=1:N                                          %%% in each subchannel
    ss(i,:)=source(Lb);                            %%% generate the bits in each subchannel
    %%%******* convolutional encoding****************
    %T=poly2trellis(9,[557 663 711]);              %%% IS-95 uplink ecoding
    %ss_enco(i,:)=convenc(ss(i,:),T);
    ss_enco(i,:)=ss(i,:);                          %%% (for no channel coding scenario)
    %%%**********************************************
    %%%******* time interleaving**********************
    %ti_ss(i,:)=Block_interleaving(20,36,ss(i,:));  
    ti_ss(i,:)=ss_enco(i,:);                       %%% (for no interleaving scenario)
    %%%**********************************************
    %%%********* modulation and spreading ***********
    ss_modu=sq64QAM_Gray_map(ss_enco(i,:));        %%% modulation, column vector 
    sp_modu=ss_modu*ss_code;                       %%% spreading, Lc*P matrix of signals for a subchannel
    sp_sig=reshape(sp_modu.',1,[]);                %%% spreaded signals in each subchannel, row vector (1*Lc)
    sm=[sm
        sp_sig];                                   %%% sm is a N*Lc matrix 
    %%%**********************************************
end
clear ss_enco ss_modu sp_modu sp_sig
[A,Lc]=size(sm);                                   %%% The symbol number in each subchannel is Lc
                                                   %%% (Lc=K*P for coherent,and Lc=(K+1)*P for incoherent)
%***************************************************************************

%***************** generate multiple access interference********************
MAI=zeros(N,Lc);
for i=1:Ku-1
    si_code=scode(scid_i(i),:);                    %%% row vector
    infm=[];
    for i=1:N                                      %%% in each subchannel
        ssi=source(Lb);
        %%%******* time interleaving**********************
        %ti_ssi(i,:)=Block_interleaving(20,36,ss(i,:));      
        %%%**********************************************
        %%%******* convolutional encoding***********************
        %T=poly2trellis(9,[557 663 711]);           %%% IS-95 uplink ecoding
        %ssi_enco=convenc(ti_ssi,T);
        ssi_enco=ssi;                               %%% (for no channel coding scenario)
        %%%*****************************************************
        spi_modu=sq64QAM_Gray_map(ssi_enco)*si_code;    %%% Lc*P matrix
        spi_sig=reshape(spi_modu.',1,[]);
        infm=[infm
              spi_sig];                    
    end
    MAI=infm+MAI;
    clear si_code ssi ssi_enco spi_modu spi_sig infm
end
%****************************************************************************

%*************************** add MAI  ***************************************
DATA=sm+MAI;                                       %%% (P*M)*Lc matrix
%DATA=sm;                                          %%% (for no MAI scenario)
%clear sm MAI
%****************************************************************************
%length(find(DATA==0))
%aa
%*********************frequency interleaving*********************************
%fi=frequency_interleaving(M,P,DATA);
fi=DATA;                                           %%% (for no interleaving scenario)
%****************************************************************************
%***************************IFFT, add prefix*********************************
for i=1:Lc                               
    mcds_modu(:,i)=sqrt(N)*ifft(fi(:,i));          %%% multicarrier modulation, N*Lc matrix
    apre(:,i)=prefix(mcds_modu(:,i),Lp);           %%% Add cyclic prefix; Ls*Lc matrix
end                                                %%% Note that the power is reduced to 1/N after IFFT transform
%****************************************************************************

%***************************** P/S ******************************************
trans=reshape(apre,1,[]);                          %%% transmitterd signal
%clear mcds_modu apre
%****************************************************************************

for si=1:length(Eb_N0)
    L_P=0;                                         %%% initial counter                      
    for i=1:nloop
        %************************* channel ***********************************
        %%%********************** fading channel *****************
        fout=sfading_channel(trans,L_P);           %%% frequency selective Rayleigh fading channels
        L_P=L_P+length(trans);                     %%% make the channel samples continuous
        %fout=trans;                               %%% (for no fading scenario)
        %%%*******************************************************
        
        %%%************* ideal channel estimation **************** 
        tc=reshape(fout,Ls,[]);                    %%% S/P, tc is a Ls*Lc matrix 
        for m=1:Lc
            rece(:,m)=deprefix(tc(:,m),Lp);        %%% remove prefix; rece is a N*Lc matrix
            dere(:,m)=fft(rece(:,m))./sqrt(N);
        end 
        H=dere./DATA ;                             %%% Note:If Ku is even, DATA may have zero
        clear tc rece dere
        %%%*******************************************************
        
        %%%*********************** AWGN channel ******************
        [cout,Pn]=AWGN(fout,Eb_N0(si),Ps,k,R,P);           %%% AWGN channels
        %cout=fout;                                        %%% (for no AWGN noise scenario)
        %%%*******************************************************
        %***********************************************************************
        
        %**************** remove prefix,FFT ************************************
        sout=reshape(cout,Ls,[]);                          %%% S/P, tc is a Ls*Lc matrix     
        for m=1:Lc
            pout(:,m)=deprefix(sout(:,m),Lp);              %%% remove prefix; pout is a N*Lc matrix
            mcds_demodu(:,m)=fft(pout(:,m))./sqrt(N);               %%% mc demodulation
        end                                                %%% ,each row is symbols for a subchannel
        %***********************************************************************
        
        %*********************frequency deinterleaving**************************
        %defi=frequency_deinterleaving(M,P,mcds_demodu);
        defi=mcds_demodu;                                    %%% ( for no interleaving scenario)
        %***********************************************************************
        
        for t=1:N
            branch=defi(t,:);                  %%% 1*Lc vector, signal matrix
            h=H(t,:);                          %%% 1*Lc vector, channel matirx
            %********* de-spreading and combination ****************************
            comb=mcds_cdma_comb_ideal(branch,ss_code,h,P);
            %*******************************************************************
            %********************** demodulation *******************************
            demap(t,:)=sq64QAM_Gray_demap(comb);           %%% demodultion in each subchannel
            %*******************************************************************
            %************ time deinterleaving***********************************
            %td_demap=Block_deinterleaving(20,36,demap);
            %*******************************************************************
            %***********************Viterbi decoding*****************************
            %decode(t,:)=vitdec(td_demap(t,:),T,9,'trunc','hard');    %%% decoding
            decode(t,:)=demap(t,:);                        %%% ( for no coding/decoding)
            %********************************************************************
        end
        errnum(si,i)=length(find(decode-ss));              %%% count the number of eror bits
        clear mcds_demodu demap
    end
    
     BER(si)=sum(errnum(si,:))/(N*Lb*nloop);
 end
 
BER
%save mcds_cdma_BER Eb_N0 BER
semilogy(Eb_N0,BER,'r-*')
grid on 
hold on