%%%%%%%%%%%% Joint Channel Estimation and Multiuser Detection for Multipath CDMA Systems %%%%%%%%%%%%
% This algorithm based on the MOE criteria was developed by Zhengyuan Xu,
% Michail K. Tsatsanis in the following literatures:
%[1] On minimum output energy CDMA receivers in the presence of multipath.
%    Proc. CISS-97, Baltimore, MD, Mar. 19-21, 1997.
%[2] Performance analysis of minimum variance CDMA receivers.
%    IEEE Trans. Signal Processing, vol. 46, No. 11, Nov. 1998.
%[3] Blind adaptive algorithms for minimum variance CDMA receivers.
%    IEEE Trans. Commun., vol. 49, No. 1, Jan. 2001.
%[4] Blind multiuser detection: from MOE to subspace methods.
%    IEEE Trans. Signal Processing, vol. 52, No. 2, Feb. 2004.

clear
load channelinva_2
% initialize paramters
UserNum = 10;                       % the number of active users
PathNum = 4;                        % the number of resolvable paths, here all the users have the same PathNum
DataNum = 3000;                     % the number of simulation iteration
LoopNum = 5;                       % the number of loop
ProcGain = 31;                      % the length of spreadingcode
SourceGain = [0 5*ones(1,UserNum-1)];       % the power of each user from which we can see the desired user's power is 5dB less than the other user's
SNR = 15;

% generate the spreadingcode (Gold Code) for each user
SpreadingCode = goldcode(log2(ProcGain+1));            % this spreadingcode has been normalized to unit power
SpreadingCode = SpreadingCode(:,2:UserNum+1)/sqrt(ProcGain);

% begin simulation
for nn = 1:LoopNum
    % generate QPSK source signal 
    for n = 1:UserNum
        Source(n,:) = make_qam(4,DataNum,1)*exp(j*pi/4);
    end
    Source = diag(10.^(SourceGain/20))*Source;
  
    % generate complex additive Gaussian noise for Zhengyuan's algorithm
    % and mine, the difference here lies in the length of noise vector
%     Noise = (10^(-SNR/20))*(randn(DataNum*ProcGain,1)+j*randn(DataNum*ProcGain,1))/sqrt(ProcGain);
    Noise = (10^(-SNR/20))*(randn(DataNum*(ProcGain+PathNum-1),1)+j*randn(DataNum*(ProcGain+PathNum-1),1))/sqrt(ProcGain);
    
    % generate Rayleigh fading channel: the pre-defined channel and
    % stotistich one. the later one was used in Zhengyuan's paper
    Channel = H(:,:,nn)/1.3705;
%     Channel = (randn(PathNum,UserNum)+j*randn(PathNum,UserNum))/sqrt(2*PathNum);
    
    % delays for all users: stotistich delays and no delay
    delay = [0 randsrc(1,UserNum-1,0:1:30)];
%     delay = zeros(1,UserNum);
    
    % compose the Spreading Code Matrix and Compound Channel response vector,
    % which means the convolution effect of Channel and Spreading Code
    for n = 1:UserNum
        for m = 1:PathNum
            CodeMatrix(:,m,n) = [zeros(m-1,1);SpreadingCode(:,n);zeros(PathNum-m,1)];
        end
        ComChan(:,n) = CodeMatrix(:,:,n)*Channel(:,n);
    end
    
    % algorithms' parameters for Zhengyuan's model
    lambda = 0.998;
    mu_u = 0.02;
    mu_f = 0.02;
    mu_g = 0.004;
    
    f = zeros(ProcGain+PathNum-1,1);
    g = [1 0 0 0]';
    u = zeros(ProcGain-1,1);
    P = 1000*eye(ProcGain+PathNum-1);
    C = CodeMatrix(:,:,1);
    PI = eye(ProcGain+PathNum-1)-C*inv(C'*C)*C';
    B = inv(C'*C);
    [U,sigma,V] = svd(C);
    Cn = U(:,PathNum+1:end);
        
    for i = 2:DataNum-1
        for n = 1:UserNum
            if (delay(n) >= PathNum-1)
                Recei(:,n) = [zeros(delay(n),1);ComChan(1:ProcGain+PathNum-1-delay(n),n)]*Source(n,i) ...
                    +[ComChan(ProcGain-delay(n)+1:end,n);zeros(ProcGain-delay(n),1)]*Source(n,i-1);
            else
                Recei(:,n) = [zeros(delay(n),1);ComChan(1:ProcGain+PathNum-1-delay(n),n)]*Source(n,i) ...
                    +[ComChan(ProcGain-delay(n)+1:end,n);zeros(ProcGain-delay(n),1)]*Source(n,i-1) ...
                    +[zeros(ProcGain+delay(n),1);ComChan(1:PathNum-1-delay(n),n)]*Source(n,i+1);
            end
        end

    % algorithms for my model
%     f = zeros(ProcGain,1);
%     g = [1 0 0 0]';
%     u = zeros(ProcGain-PathNum,1);
%     P = 1000*eye(ProcGain);
%     C = CodeMatrix(1:ProcGain,:,1);
%     PI = eye(ProcGain)-C*inv(C'*C)*C';
%     B = inv(C'*C);
%     [U,sigma,V] = svd(C);
%     Cn = U(:,PathNum+1:end);

%     for i = 3:DataNum
%         for n = 1:UserNum
%             if (delay(n) <= ProcGain-PathNum+1)
%                 Recei(:,n) = [ComChan(ProcGain-delay(n)+1:end,n);zeros(ProcGain-delay(n)-PathNum+1,1)]*Source(n,i-1) ...
%                     +[zeros(delay(n),1);ComChan(1:ProcGain-delay(n),n)]*Source(n,i);
%             else
%                 Recei(:,n) = [ComChan(end-(delay(n)-ProcGain+PathNum-1)+1:end,n);zeros(2*ProcGain-delay(n)-PathNum+1,1)]*Source(n,i-2) ...
%                     +[ComChan(ProcGain-delay(n)+1:2*ProcGain-delay(n),n)]*Source(n,i-1) ...
%                     +[zeros(delay(n),1);ComChan(1:ProcGain-delay(n),n)]*Source(n,i);
%             end
%         end
        r = sum(Recei,2)+Noise((i-1)*(ProcGain+PathNum-1)+1:i*(ProcGain+PathNum-1));                % the received data vector composed of all the active user's information
            %%%%%%%%%%%%% Method I %%%%%%%%%%%%%%%%%%%%
%         f = PI*(f-mu_f*r*r'*f)+C*inv(C'*C)*g;
%         x = g-mu_g/mu_f*inv(C'*C)*(C'*(f-mu_f*r*r'*f)-g);
%         a1 = mu_g^2*g'*g;
%         a2 = mu_g*(g'*x+x'*g);
%         a3 = x'*x-1;
%         if (a2^2-4*a1*a3 >= 0)
%             rho = 1/(2*a1)*(-a2-sqrt(a2^2-4*a1*a3));
%             g = mu_g*rho*g+x;
%         else
%             g = g;
%         end
%         %%%%%%%%%%%%% Method II %%%%%%%%%%%%%%%%%%%%
%         f = PI*(f-mu_f*r*r'*f)+C*inv(C'*C)*g;
%         g = g+mu_g/mu_f*(eye(PathNum)-g*g'/(g'*g))*inv(C'*C)*(mu_f*C'*r*r'*f+g-C'*f);
%         g = g/norm(g,'fro');
        %%%%%%%%%%%%% Method III %%%%%%%%%%%%%%%%%%%%
%         y1 = C'*r;
%         y2 = Cn'*r;
%         w = y1'*B*g-y2'*u;
%         u = u+w*mu_u*y2;
%         g = g+w*mu_g*B*y1;
%         g = g/norm(g,'fro');
%         f = C*B*g-Cn*u;
%         %%%%%%%%%%%%% Method IV %%%%%%%%%%%%%%%%%%%%
        P = 1/lambda*(P-P*r*r'*P/(lambda+r'*P*r));
        [a,b] = eig(C'*P*C);
        g = a(:,find(diag(b) == min(diag(b))));
% %         f = min(diag(b))*P*CodeMatrix(1:ProcGain,:,1)*g;
%         f = min(diag(b))*P*CodeMatrix(:,:,1)*g;
%         SINR(nn,i) = 10*log10(abs(2*f'*ComChan(1:ProcGain,1)*ComChan(1:ProcGain,1)'*f/(abs(f'*(r-Source(1,i)*ComChan(1:ProcGain,1)))^2)));
        EstChan = g;                          % the estimated channel response vector of the user of interest
        error(nn,i) = norm(EstChan*exp(-j*angle(EstChan(1)))*exp(j*angle(Channel(1,1)))-Channel(:,1)/norm(Channel(:,1),'fro'),'fro')^2;
    end
end
% Estimation Error Figure
semilogy(mean(error),'y')
xlabel('Iteration Number')
ylabel('MSE of Channel Estimation')
title('50 Average Result of Blind Adaptive Channel Estimation ')
% figure
% semilogy(mean(SINR),'g')