%% Rake receiver of CDMA
%% may need call lcode.mat to load PN sequence (lcode)
%% Copyright: Xiaohua(Edward) Li, Assistant Professor%            Department of Electrical and Computer Engineering%            State University of New York at Binghamton%            http://ucesp.ws.binghamton.edu/~xli% June 2003%P=64; K=5; Nd=500; % spreading gain, # of users, symbol amount,
Lh=6;              % channel length,
L1=3*Lh; fd=round(L1/2); % equalizer length, equalization delay
dB=20; Pow=[1 ones(1,K-1)*10^(0/20)];  % SNR, Near-far ratio
tao=[0 floor(rand(1,K-1)*P)];   % asynchronous or user delay

%load lcode.mat;           % load the pre-computed code sequence
disp('Calculate PN code. To speed up, you can save lcode and reload')  
% or you can compute it here 
lcode=longcode(50000);     % it is better to "save lcode.mat lcode"

b=(round(rand(K,Nd))*2-1)+sqrt(-1)*(round(rand(K,Nd))*2-1); % symbols
sf=(0:1:K-1)*100;   % users use the same long code with various off-set

h=randn(K,Lh)+sqrt(-1)*randn(K,Lh);       % channels for K users
for i=1:K, h(i,:)=h(i,:)/norm(h(i,:));  end  % normalize channels
Mh=zeros(L1,Lh+L1-1);                   % User 1's Channel matrix H1
for i=1:L1,   
    Mh(i,i:i+Lh-1)=h(1,:);  
end
   
%%%%%%%%%% Construct received signal samples in CDMA  
x=lcode(1:Nd*P); x1=zeros(size(x));      % construct received samples
for i=1:Nd,               % signal from user 1: spreading part
    x((i-1)*P+1:i*P)=b(1,i)*x((i-1)*P+1:i*P); 
end 
x=filter(h(1,:),1,x);     % channel filtering part
for j=2:K                 % signal from other users
   for i=1:Nd, 
      x1((i-1)*P+1:i*P)=...  % spreading & near-far power
         Pow(j)*b(j,i)*lcode(sf(j)+(i-1)*P+1:sf(j)+i*P);
   end
   x1=filter(h(j,:),1,x1);   % channel filtering
   x(tao(j)+1:Nd*P)=x(tao(j)+1:Nd*P)+x1(1:Nd*P-tao(j));  % asynchrous delay
end

vn=randn(size(x))+sqrt(-1)*randn(size(x));    % AWGN
vn=vn/norm(vn)*10^(-dB/20)*norm(x);
SNR=20*log10(norm(x)/norm(vn)),              % check SNR of recieved samples
x=x+vn;                           % received noisy signals

%%%%%%%%%%%%% Construct sample matrices X and despreaded matrices Y
Y=zeros(L1,Nd-3); X=zeros(L1,P); 
for i=1:Nd-3
   for j=1:P,
      if i*P+j+fd-L1>0         % need carefully adjust the timing
          X(:,P-j+1)=x(i*P+j-1+fd:-1:i*P+j+fd-L1).'; 
      else X(1:i*P+j-1+fd,P-j+1)=x(i*P+j-1+fd:-1:1).';
      end
   end
  
   u=X*lcode((i+1)*P:-1:i*P+1).'/P;  % desreading by detectors from codesp
   Y(:,i)=u;                         % save to a matrix of detector output
%   if i/100==round(i/100), i/10000, end
end


%%%%%%% Rake Receiver: channel and symbol estimation
R=Y*Y'/(Nd-3);    % channel estimation by SVD 
[Ur,Sr,Vr]=svd(R);
hb1=Ur(:,1); hb1=hb1/(hb1'*Mh(:,fd)/abs(hb1'*Mh(:,fd)));
Channel_MSE=norm(hb1-Mh(:,fd))   % MSE of channel estimation

sy=hb1'*Y;                   % symbol estimation by equalizer

%%%%%%% Check performance 
fh=hb1'*Mh; i=find(abs(fh)==max(abs(fh))); 
sy1=sy/fh(i);          % scaling and ambiguity phase removing

sy2=sign(real(sy1))+sqrt(-1)*sign(imag(sy1));
temp=sy2-b(1,2:Nd-2); 
SER=length(find(temp~=0))/length(temp)    % calculate SER (i.e., BER)

if 1
    subplot(221), plot(b(1,:),'o'), grid, axis([-2 2 -2 2])
    title('Transmitted Symbols of User 1'), 
    xlabel('Real'), ylabel('Imag')
    
    subplot(222), plot(x,'o'), grid
    title('Received Noisy Signals')
    xlabel('Real'), ylabel('Imag')
    
    subplot(223), plot(sy,'o'), title('Equalized symbols of User 1'),grid
    xlabel('Real'), ylabel('Imag')
    
    subplot(224), plot(sy1,'o'), grid
    title('Scaled Equalized Symbols'),
    xlabel('Real'), ylabel('Imag')
end