%% Code to generate bartlett beamformer spectrum

% This code generates the bartlett spectrum for AoA estimation (Fig5a)
% In case of any queries, please contact R.S. Prasobh Sankar (prasobhsankar1@gmail.com)

% Since this code does not involve monte-carlo runs, there is no
% corresponding code to plot Fig5a from data.

clear all
close all
clc


aspread_aoa = 60;% AoA is drawn randomly from -aspread_aoa to aspread_aoa
no_sim = 1; % Monte carlo runs
case_val = 1; %  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 =  0;
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


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

lev_s1 = 3*sqrt((SNR(1)+1)/2); 


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

doa_true = [10];


% 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';  % channel matrix      


% data transmission

S = w_opt*sqrt(SNR).*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); 

b2 = 1;  % arbitrary voltage level
Z = sigma_delta_ADC(X,psi,d_r,b2); 

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

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

    w_ref = sqrt(1/(N_r^2)).*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);
    cost_sd2(i) = abs(w_ref'*R_z*w_ref);

end



figure
plot(aoa_grid, 10.*log10(cost_a), 'b--','LineWidth',1.2) 
hold on
plot(aoa_grid, 10.*log10(cost_sd), 'r-.','LineWidth',1.2) 
hold on
plot(aoa_grid, 10.*log10(cost_sd2), 'k:','LineWidth',1.2) 
grid on
legend('UQ','SD c=proposed','SD c=1')
xlabel('\theta (degrees)')
ylabel(' J(\theta) (Bartlett spectrum) (dB)')




%%
% savefig('Bartlett_spec.fig')
% saveas(gcf,'Bartlett_spec','epsc')








