% Program 
%
% Simulation program to realize DS-CDMA system
%
% dscdmabpsk.m
%
% Programmed by XIKAI HU
% 
clear;
%******************** Preparation part **********************

sr=256000.0; % Symbol rate
ml=1;        % Number of modulation levels
br=sr.*ml;   % Bit rate (=symbol rate in this case)
nd = 100;
N=nd;% Number of symbols that simulates in each loop
ebn0=10;     % Eb/N0
IPOINT=8;    % Number of oversamples

%******************* Filter initialization ********************

irfn=21;     % Number of filter taps          
alfs=0.5;    % Rolloff factor
[xh] = hrollfcoef(irfn,IPOINT,sr,alfs,1);   %Transmitter filter coefficients 
[xh2] = hrollfcoef(irfn,IPOINT,sr,alfs,0);  %Receiver filter coefficients 


%******************* 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/IPOINT; 

% 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];

% 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];

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

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

% Number of fading counter to skip 
itnd0=nd*IPOINT*100;

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

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

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

% 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;

m=8;    %擴(kuò)頻碼的碼片數(shù)目


%********************** Spreading code initialization **********************

user  = 1;                                                          % number of users
seq   = 3;                                                          % 1:M-sequence  2:Gold  3:Orthogonal Gold
stage = 3;                                                          % number of stages
ptap1 = [1 3];                                                      % position of taps for 1st
ptap2 = [2 3];                                                      % position of taps for 2nd
regi1 = [1 1 1];                                                    % initial value of register for 1st
regi2 = [1 1 1];                                                    % initial value of register for 2nd

%******************** Generation of the spreading code *********************

switch seq
case 1                                                              % M-sequence
    code = mseq(stage,ptap1,regi1,user);
case 2                                                              % Gold sequence
    m1   = mseq(stage,ptap1,regi1);
    m2   = mseq(stage,ptap2,regi2);
    code = goldseq(m1,m2,user);
case 3                                                              % Orthogonal Gold sequence
    m1   = mseq(stage,ptap1,regi1);
    m2   = mseq(stage,ptap2,regi2);
    code = [goldseq(m1,m2,user),zeros(user,1)];
end
code = code * 2 - 1;
clen = length(code);

%**************************** START CALCULATION ****************************

nloop = 1000;                                                       % simulation number of times
noe   = 0;
nod   = 0;

for iii=1:nloop
    
%****************************** Transmitter ********************************
    data = rand(user,nd*ml) > 0.5;
    
%******************** BPSK Modulation ***********************  

      data1=data.*2-1   %————————測(cè)試點(diǎn)————————

%********************擴(kuò)展原始碼元以實(shí)現(xiàn)擴(kuò)頻序列*************
%for i=1:N
    %for j=(m*(i-1)+1):m*i
      % data1((m*(i-1)+1):m*i)=data1(i)*code(1,:)      %第一路天線需要發(fā)射的信號(hào)————————測(cè)試點(diǎn)
%end
data1=spread(data1,code)
%data8=despread(data1,code)%-------可以用于檢測(cè)
[data2] = oversamp( data1, N*m , IPOINT);                 %128  -1,0,1
[data3] = conv(data2,xh); % conv: built in function--------測(cè)試點(diǎn)  295


%****************** Attenuation Calculation *****************
	
    spow=sum(data3.*data3)/nd;%————————————————————————？？？？
	attn=0.5*spow*sr/br*10.^(-ebn0/10);
	attn=sqrt(attn);
    
%********************** Fading channel **********************

% Generated data are fed into a fading simulator
% In the case of BPSK, only Ich data are fed into fading counter
  [ifade,qfade]=sefade(data3,zeros(1,length(data3)),itau,dlvl,th1,n0,itnd1,now1,length(data3),tstp,fd,flat);
  
% Updata fading counter
  itnd1 = itnd1+ itnd0;


%************ Add White Gaussian Noise (AWGN) ***************
	
    inoise=randn(1,length(ifade)).*attn;  % randn: built in function
	data4=ifade+inoise;  %-----------------------------------------------------295
	data5=conv(data4,xh2);  % conv: built in function--------------------------462

	sampl=irfn*IPOINT+1;
	data6 = data5(sampl:IPOINT:IPOINT*nd*clen+sampl-1);  %------------------------------------16
   
    
    
%****************************解擴(kuò)信號(hào)************************
             data7=despread(data6,code)

%******************** BPSK Demodulation *********************
	
    demodata=data7 > 0;

%******************** Bit Error Rate (BER) ******************
	
    % count number of instantaneous errors
    noe2=sum(abs(data-demodata));  % sum: built in function
	
    % count number of instantaneous transmitted data
    nod2=length(data);  % length: built in function
	
    noe=noe+noe2
	nod=nod+nod2

	fprintf('%d\t%e\n',iii,noe2/nod2);
%end % for iii=1:nloop    

%********************** Output result ***************************

ber = noe/nod;
fprintf('%d\t%d\t%d\t%e\n',ebn0,noe,nod,noe/nod);


%******************** end of file ***************************