%% Code to generate quantization noise power spectrum
% This code generate data for Fig2b

% Performing large number of Monte Carlo runs may lead to excessive runtime.
% Please use Fig2b.m to plot the quantization noise power spectrum from
% provided data.

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

clear all
close all
clc
psi_set = [-15; 0; 15];
nosim = 100 ;

for outerloop = 1:length(psi_set)
    

OSR = 4; % spatial oversampling ration at BS
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 BS
d_t = 0.5; % inter-element spacing at UE

SNR_dB = 0;
psi = psi_set(outerloop) ; % steering angle
L=1; % single path scenario

SNR = 10.^(0.1.*SNR_dB);


aspread_aoa = 30;
spacing_aoa = 20;


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

%%%%%
w_opt = [1; zeros(N_t-1,1)]; % precoder for AoA estimation - isotropic transmission
T_1 = 10;
aoa_grid = -60:1:60;

q_noise = zeros(1,length(aoa_grid));

for iter=1:length(q_noise)        
    err= 0;   
    
    a_ref = gen_a(N_r,d_r,aoa_grid(iter));    
    lev_s1 =  3*(sqrt(SNR(1)+1)/sqrt(2) ) ;
   
    
    for loop_inner = 1:nosim
        
        % 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);
        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    
        
        S = w_opt*sqrt(SNR(1)).*ones(1,T_1);
        S_true = S;
        
        N_noise = sqrt(1/2).*(randn(N_r,T_1) + (1j).*randn(N_r,T_1));

        X = H_true*S + N_noise; % received signal (unquantized)
        b = lev_s1; % quantization voltage level (b)
        Y = sigma_delta_ADC(X,psi,d_r,b); % output of sigma-delta quantizer
        err = err +  var( a_ref'*(Y - X) ); % Variance of quantization noise
        
              
       [outerloop iter loop_inner] 
    end
    
    q_noise(iter) = err/nosim;
       
end

%%
figure
plot((aoa_grid), 10*log10(q_noise), 'r','LineWidth',1.2)
xlabel('Angle of arrival')
ylabel('Quantization noise (dB)')
grid on




end