%% MU-MIMO channel estimation
% This function performs channel estimation in 
% mmWave MU-MIMO systems equipped with spatial sigma-delta ADCs at the BS.
% This code implements Algorithm 2 in [Add TSP citation]
% 
% 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
% b - Quantizer voltage level
% R_n_pw - Pre-whitening matrix, R_n_pw = inv(R_n), where R_n is the
% covariance matrix of the effective noise term [Add equation number]

%%%%%%%%%%% 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_SD(H, S, T, d, L, SNR, aoa_grid, b, R_n_pw )

[N_r K] = size(H);  % K single antenna UEs and N_r antenna BS
aoa_est = zeros(L,K);
alpha_est = zeros(L,K);

for ind = 1:K 
    
    cost1 = zeros(1,length(aoa_grid));
    Y_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 unquantized signal
        Y = sigma_delta_ADC(X,0,d,b); % output of sigma-delta quantizer

        % Bartlett beamformer spectrum computation
        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'*(   Y*S(ind,:)' )  ) )^2; 
        end
        Y_pg = Y_pg + R_n_pw*Y*S(ind,:)' ; 
    end

    cost1 = cost1/T;
    Y_pg = Y_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
    Psi = (SNR*sqrt(1/L)).*R_n_pw*A_a_rec; 
    alpha_est(:,ind) = Psi\Y_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


