%% Code to generate data for Fig4a, Fig4b, Fig5b and Fig5c
% This code illustrate the estimation of AoAs and complex path gains using
% the proposed algorithm

% Performing large number of Monte Carlo runs may lead to excessive runtime.
% Please use Fig4a.m, Fig4b.m, Fig5b.m, and Fig5c.m under "Codes for plotting" to obtain the plot from provided data.

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

aspread_aoa = 30;% AoA is drawn randomly from -aspread_aoa to aspread_aoa
no_sim = 100; % Monte carlo runs
case_val = 0; %  Choice of voltage levels
% case_val = 0 => use proposed voltage levels 
% case_val = 1 => use c = 1
% case_val = 2 => use c corresponding to 0dB


OSR = 4; % spatial oversampling ratio
N_t = 32; % number of antennas at UE
N_r = (128/OSR)*OSR;  % number of antennas at BS
d_r = 0.5/OSR;  % inter-element spacing at the BS
d_t = 0.5; % inter-element spacing at the UE

SNR_dB =  -16:4:16;
SNR = 10.^(0.1.*SNR_dB);

aoa_grid = -90:1:90; % Grid for AoA estimation
spacing_aoa = 20; % minumum spacing between AoAs

% codebook parameters
M = 128;
u_grid = linspace(-1,1,M);
% k th grid point corresponds to : sin^-1( -1 + (2/(M-1))(k-1)   ) degrees
spacing_aod = 20; % minimum spacing for AoDs

edge_skip = 3;
% edge_skip = 3 => considered grid points: 4 to M-3
% edge_skip = 3 => AoDs lie in  approximately -75^0 to 75^0 

l_eff = M - 2*edge_skip;

L = 1; % number of paths
T_1  = 10; % Number of snapshots for AoA estimation

tic              

MSE_a_aoa_det = zeros(1, length(SNR));
MSE_sd_aoa_det = zeros(1,length(SNR));


MSE_a_gain = zeros(1, length(SNR));
MSE_sd_gain = zeros(1,length(SNR));


w_opt = [1; zeros(N_t-1,1)]; % precoder for AoA estimation - isotropic transmission
psi = 0; % steering angle
U = tril(ones(N_r));
U_inv = inv(U);

for iter=1:length(SNR)
    
   
    err_a_aoa = 0;
    err_sd_aoa = 0;    
    gain_dr = 0;    
    err_a_gain = 0;
    err_sd_gain = 0;
    
    if(case_val == 0)
        lev_s1 = 3*sqrt((SNR(iter)+1)/2); 
    elseif(case_val == 1)
        lev_s1 = 1;
    elseif(case_val == 2)
        lev_s1 = 3*sqrt((1+1)/2);
    end
      
    R_n = eye(N_r) + ((2*lev_s1^2)/3).*U_inv*U_inv'; 
    R_n_pw = inv(sqrtm(R_n)); % pre-whitening matrix
    
    for loop_inner = 1:no_sim
        
        % channel realization
       
        % unit modulus path gains (Pure LoS scenario)
        path_gain_angle = 2*pi*rand(L,1) - pi;
        alpha = ones(L,1).*exp(1j.*path_gain_angle);         
        
        % generate Aoa
        n=2*aspread_aoa;
        k=L;
        b=spacing_aoa;  % Min. spacing 
        [as,is]=sort(randperm(n-(k-1)*(b-1),k)); 
        a = as + (0:k-1)*(b-1);
        a = a(is) - aspread_aoa;
        doa_true = sort(a);
        
        
        % generate AoD
        n=l_eff;
        k=L;
        b=spacing_aod;
        [as,is]=sort(randperm(n-(k-1)*(b-1),k)); 
        a = as + (0:k-1)*(b-1);
        aod_true_dc_ind = edge_skip + sort(a);
        aod_true_dc = u_grid(aod_true_dc_ind);
                
        A_TX = gen_a_v2(N_t,d_t, aod_true_dc);
        A_RX = gen_a(N_r,d_r, doa_true);     
        H_true = sqrt(1/L).*A_RX*diag(alpha)*A_TX';  % channel matrix      
        gain_dr = gain_dr + norm(alpha)^2;
        
        % data transmission
        
        S = w_opt*sqrt(SNR(iter)).*ones(1,T_1);
        S_true = S; % pilot              
        N_noise = sqrt(1/2).*(randn(N_r,T_1) + (1j).*randn(N_r,T_1));
        X = H_true*S + N_noise; % received unquantized signal
        X_true = X;
        
       % Step1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
           
        b = lev_s1;
        Y = sigma_delta_ADC(X,psi,d_r,b); 
        Y_true = Y;

        R_x = (1/T_1).*X*X';
        R_y = (1/T_1).*Y*Y';
        cost_a = zeros(1,length(aoa_grid));
        cost_sd = zeros(1,length(aoa_grid));

        % computing Bartlett beamformer spectrum
        for i=1:length(aoa_grid)

            w_ref = gen_a(N_r, d_r, aoa_grid(i));
            cost_a(i) = abs(w_ref'*R_x*w_ref);
            cost_sd(i) = abs(w_ref'*R_y*w_ref);

        end
 
        % UQ        
        [val pos] = local_max(cost_a,L); 
        aoa_est_a = sort(aoa_grid(pos));
        A_RX_est_a = gen_a(N_r, d_r, sort(aoa_est_a));        
        err_a_aoa = err_a_aoa + (nnz(sort(doa_true) ~= sort(aoa_est_a)) > 0);
         
        % SD         
       [val pos] = local_max(cost_sd,L);
       aoa_est_sd = sort(aoa_grid(pos));
       A_RX_est_sd = gen_a(N_r, d_r, sort(aoa_est_sd));            
       err_sd_aoa = err_sd_aoa + (nnz(sort(doa_true) ~= sort(aoa_est_sd)) > 0);
       
       % %%%%%  path gain estimation      
              
       % UQ     
       G1 = sqrt(SNR(iter))*ones(L,T_1);
       D = sqrt(1/L).*krp( transpose(G1),A_RX_est_a  );
       x = reshape(X_true, T_1*N_r,1);
       alpha_est_a = D\x;
       
       err_a_gain =  err_a_gain + norm(alpha - alpha_est_a)^2;
       
       % SD     
       D_sd = sqrt(1/L).*krp( transpose(G1),R_n_pw*A_RX_est_sd  );
       y_pw = reshape(R_n_pw*Y_true, T_1*N_r,1);
       alpha_est_sd = D_sd\y_pw; 
       err_sd_gain =  err_sd_gain + norm(alpha - alpha_est_sd)^2;
                   
              
      [iter loop_inner]
        
    end
   
    MSE_a_aoa_det(iter) = err_a_aoa/no_sim;
    MSE_sd_aoa_det(iter) = err_sd_aoa/no_sim;
        
    MSE_a_gain(iter) = err_a_gain/gain_dr;
    MSE_sd_gain(iter) = err_sd_gain/gain_dr;
    
    
end


toc

%% Plotting

figure
semilogy(SNR_dB, (MSE_a_aoa_det), 'r','LineWidth',1.2)
hold on
semilogy(SNR_dB, (MSE_sd_aoa_det), 'b','LineWidth',1.2)
grid on 
xlabel('SNR (dB)')
ylabel(' \theta_{err} ')
legend('UQ','SD')
title('Aoa estimation error')
ylim([10^(-2) 1.1])

figure
plot(SNR_dB, 10*log10(MSE_a_gain), 'r','LineWidth',1.2)
hold on
plot(SNR_dB, 10*log10(MSE_sd_gain), 'b','LineWidth',1.2)
grid on 
xlabel('SNR (dB)')
ylabel(' path gain NMSE (dB) ')
legend('UQ','SD')
title('Path gain estimation error')








