%% Correlation between input and quantization noise
% This code generate data for Fig2a

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

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


clear all
close all
clc

no_sim = 100; % number of Monte-Carlo runs
psi_set = [-15 0 15];

for outerloop = 1:length(psi_set)

OSR = 4; % spatial oversampling ration at BS
d_r = 0.5/OSR; % inter-element spacing at BS
d_t = 0.5; % inter-element spacing at UE
N_r = (128/OSR)*OSR;   %Number of receive antennas
N_t = (32/OSR)*OSR;    %Number of transmit antennas


N = N_t; %number of snapshots
SNR_dB = 0; %-15:5:20;
SNR = 10.^(0.1.*SNR_dB);

aspread_aoa = 30;% AoA is drawn randomly from  -aspread_aoa to aspread_aoa
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;

psi = psi_set(outerloop); %steering angle

tic
    
R_q = zeros(N_r);
R_s = zeros(N_r);



for iter=1:no_sim
   
        % 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      

        N_noise = sqrt(1/2).*(randn(N_r,N) + (1j).*randn(N_r,N));

        S = sqrt(SNR(1)).*ones(1,N_t); 
        w_opt = [1; zeros(N_t-1,1)];
        S = w_opt*S; % pilots for AoA estimation phase

        b = 3*(sqrt(SNR(1)+1)/sqrt(2) ) ; % quantization voltage level

        X = H_true*S + N_noise; % received unquantized signal
        Y_q = b.*sgn(X); % one-bit quantized signal

        Y = sigma_delta_ADC(X, psi, d_r, b);

        % computing the correlation between input and quantization noise
        R_q = R_q + (1/N).*X*(Y_q - X)' ; % regular 1-bit quantizer
        R_s = R_s + (1/N).*X*(Y - X)' ;  % spatial sigma-delta quantizer

        [outerloop iter]

end

R_q = R_q./no_sim;
R_s = R_s./no_sim;

toc

%%


figure
plot(abs(diag(R_q)), 'r--','LineWidth',1.4);
hold on
plot(abs(diag(R_s)), 'b','LineWidth',1.4);
grid on
legend('One-bit','Sigma-Delta');
xlabel('Antenna Index');
ylabel('Correlation');
% 
% mat_title = ['A7Fig2a_',num2str(psi),'.mat'];
% save(mat_title)

end






