%% MU-MIMO channel estimation
% This function applies proposed Algorithm 2 for channel estimation in
% MU-MIMO systems without quantization. This is acts as a benchmark for the
% channel estimation in MU-MIMO systems with spatial sigma-delta
% quantizers.

% For any queries, please contact R.S. Prasobh Sankar (prasobhsankar1@gmail.com)


%%%%%%%%%%% Input parameters  %%%%%%%%%%%%%%%%%%%%

% H - true channel (To generate received signal)
% S - K \times K orthogonal pilot matrix with normalized columns
% T - number of times the pilot matrix to be transmitted (We use T=1 for
% the plots in the main paper)
% d - interelement spacing of antenna array at the BS (in terms of
% wavelengths)
% L - Number of paths per UE. We assume L_k = L for all UEs.
% SNR - Uplink SNR
% aoa_grid - search grid for AoAs


%%%%%%%%%%% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% H_est - estimated MIMO channel matrix  (N_r \times K)
% aoa_est - estimated angle of arrivals (L \times K)
% alpha_est - estimated path gains   (L \times K)



function [H_est aoa_est alpha_est] = MU_MIMO_est_UQ(H,S, T, d, L, SNR, aoa_grid )

[N_r K] = size(H);
aoa_est = zeros(L,K);
alpha_est = zeros(L,K);

for ind = 1:K 
    
    cost1 = zeros(1,length(aoa_grid));
    X_pg = zeros(N_r,1);  
    
    for k1 = 1:T  %multiple soundings


        N_noise = sqrt(1/2).*( randn(N_r,K) + 1j.*randn(N_r,K) );
        X = H*S + N_noise; % received signal without quantization

        % Computing beamformer spectrum
        for i=1:length(aoa_grid)
            a = sqrt(1/N_r).*gen_a(N_r, d, aoa_grid(i));
            cost1(i) = cost1(i) +  abs(  mean(a'*(   X*S(ind,:)' )  ) )^2;
        end
        X_pg = X_pg + X*S(ind,:)' ;
    end

    cost1 = cost1/T;
    X_pg = X_pg/T;
   
    [val pos] = local_max(cost1,L); % finding peaks of Bartlett spectrum
    aoa_est(:,ind) = sort(aoa_grid(pos)); % estimated AoAs
    A_a_rec = gen_a(N_r, d, aoa_est(:,ind)); % estimated array manifold
    B = (SNR*sqrt(1/L)).*A_a_rec; 
    alpha_est(:,ind) = B\X_pg;   % estimated path gain  
end

% Computing overall channel estimate
H_est= zeros(N_r,K);
for i=1:K
    H_est(:,i) = sqrt(1/L).*gen_a(N_r, d, aoa_est(:,i))*alpha_est(:,i);
end

end


