%% Code to generate data for Fig4c, Fig6b, and Fig6c
% This code illustrate the estimation AoDs using the proposed algorithm
% To compute the AoD estimation error, we have assumed that the estimates
% of AoAs and path gains are perfect.

% Performing large number of Monte Carlo runs may lead to excessive runtime.
% Please use Fig4c.m, Fig6b.m, and Fig6c.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
T_s = 1; % number of soundings (Sec. )
no_sim = 10^3; % 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 -10dB
% case_val = 3 => use c corresponding to Step 1

ang_steer = 1;
% ang_steer : enable/disable angle steering during AoD estimation
% 0 implies angle steering is disabled. 1 implies angle steering is enabled

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);

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

tic        
 
MSE_a_aod_det = zeros(1, length(SNR));
MSE_sd_aod_det = zeros(1,length(SNR));


for iter=1:length(SNR)
    
    err_a_aod = 0;
    err_sd_aod = 0;
    
      
        lev_s1 = 3*sqrt((SNR(iter)+1)/2); % Step 1
        if(case_val == 0)
               lev_s2 =  3*sqrt((SNR(iter)*N_t+1)/2);  
        elseif(case_val == 1)
            lev_s2 = 1;
        elseif(case_val == 2) 
            lev_s2 = 4.3474; % -10 dB            
        elseif(case_val == 3) 
            lev_s2 = lev_s1;
        end
       
    for loop_inner = 1:no_sim
        
        % channel realization
        % unit modulus path gains
        
        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);
        dod_true_dc_ind = edge_skip + sort(a);
        dod_true_dc = u_grid(dod_true_dc_ind);
        dod_true = asind(dod_true_dc);
        
        A_TX = gen_a_v2(N_t,d_t, dod_true_dc);
        A_RX = gen_a(N_r,d_r, doa_true);     
        H_true = sqrt(1/L).*A_RX*diag(alpha)*A_TX';        
        
       % Step 1 is assumed to be perfect.  
       %% Step 2: AoD estimation
   
       aod_est_a_dc_ind =  Dod_est_UQ(H_true, SNR(iter), u_grid, T_s,  A_RX, L, d_t);
       aod_est_sd_dc_ind = Dod_est_SD(H_true, SNR(iter), u_grid, T_s, doa_true, L, d_t, lev_s2, d_r, ang_steer);
       
       err_a_aod = err_a_aod + (nnz(sort(dod_true_dc_ind) ~= sort(aod_est_a_dc_ind)) > 0);
       err_sd_aod = err_sd_aod + (nnz(sort(dod_true_dc_ind) ~= sort(aod_est_sd_dc_ind)) > 0);            
              
       [iter loop_inner]
        
    end
    
    MSE_sd_aod_det(iter) = err_sd_aod/no_sim;
    MSE_a_aod_det(iter) = err_a_aod/no_sim;
  
end


toc

%% Plotting

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

