% This matlab code computes BER performance as in Fig 4.
% Yuan-Pei Lin, ypl@mail.nctu.edu.tw
%
% notation
% f0: the optimal unconstrained case
% f1: quantized version using Structure 1

% f1less: quantized version using Structure 1 with 0.5 fewer average bits
% f1less2: quantized version using Structure 1 with 1 fewer average bit
% f2: quantized version using Structure 2
% f2less: quantized version using Structure 2 with 0.5 fewer average bits
% f2less2: quantized version using Structure 2 with 1 fewer average bit

clearvars;
close all;

cl=3; %number of clusters
ray=10; %number of rays in each cluster
path=cl*ray;

t=16;  %number of transmit antennas
r=t; %number of receive antennas
ULA=0; % ULA=1 for ULA antenna array and ULA=0 for UPA antenna array


   sigma_ta=7.5*ones(1,cl)/180*pi;% standard deviation of AoD in azimuth
   sigma_te=7.5*ones(1,cl)/180*pi;% standard deviation of AoD in elevation

   sigma_ra=7.5*ones(1,cl)/180*pi;% standard deviation of AoA in azimuth
   sigma_re=7.5*ones(1,cl)/180*pi;% standard deviation of AoA in elevation


   % compute the dimension of the antenna array  in z-axis (heightt) and y-axis (Wt)
   % if the total number of antennas is not a square, the Y dimension is
   % two times the Z dimension
    if ceil(sqrt(t))==sqrt(t)        
       Wt=sqrt(t);
       heightt=Wt;
    else
        Wt=sqrt(2*t);
        heightt=Wt/2;
    end   
   

    if ceil(sqrt(r))==sqrt(r)        
       Wr=sqrt(r);
       heightr=Wr;
    else
        Wr=sqrt(2*r);
        heightr=Wr/2;
    end   
    
   
dt=1/2; %normalized transmit antenna spacing
dr=1/2; %normalized receive antenna spacing
kr=(0:r-1)';
kt=(0:t-1)';
qam=4; %QPSK symbol transmitted


N_theta=4; N_phi=N_theta; % the quantization resolution for Structure 1, N_theta=2^b_theta
N_theta2=N_theta*2; N_phi2=N_theta/2; % the quantization resolution for Structure 2
snr_range=(-20:2:10);
pt=10.^(snr_range/10);
pt=pt*3/(qam-1);
qam_const=1;
if qam>4
   qam_const=4/(1-1/sqrt(qam));
end



ber0=zeros(length(pt),1);ber1=zeros(length(pt),1); ber1less=ber1;ber1less2=ber1; ber2=ber0;
ber2_less=ber0;ber2_less2=ber0;


phi_itr=100; % number of realizations for a given set of stats
mu_itr=100; % number of different set of stats
itr=phi_itr*mu_itr; % total number of iterations

for n_mu=1:mu_itr
    mu_ta=2*pi*rand(cl,1); % mean of azimuth angle at tx
    mu_ra=2*pi*rand(cl,1);% mean of elevation angle at tx
    mu_te=2*pi*rand(cl,1); % mean of azimuth angle at rx
    mu_re=2*pi*rand(cl,1);%mean of elevation angle at rx

for n_itr=1:phi_itr

      % channel generation  
   [H, At, Ar]=mmwCh(ULA,dr,dt,t,r,Wt,heightt,Wr,heightr, path, cl,ray,mu_ta, mu_te, mu_ra, mu_re, sigma_ta,sigma_te,sigma_ra,sigma_re);

    


[V,D]=eig(H'*H);
[la,i0]=max(diag(D));
f0=V(:,i0);
[~,i0]=max(abs(f0)); % for scaling f0, extra 2 s.t. max(abs(f0))=2
f0=f0/f0(i0)*2;



f1=zeros(t,1);  f1less=zeros(t,1);f1less2=zeros(t,1);
f2=f0/1.00001;  f2_less=f2; f2_less2=f2; 
for i=1:t 
    
    f1(i)=q1(f0(i),N_theta,N_phi);  %quantization using S1 
    f1less(i)=q1(f0(i),N_theta,N_phi/2);%quantization using S1 with 0.5 fewer average bits
    f1less2(i)=q1(f0(i),N_theta/2,N_phi/2);%quantization using S1 with 1 fewer average bit
    f2(i)=q2(f2(i),N_theta2,N_phi2);  %quantization using S2
    f2_less(i)=q2(f2_less(i),N_theta,N_phi/2);%quantization using S2 with 0.5 fewer average bits
    f2_less2(i)=q2(f2_less2(i),N_theta,N_phi/4);%quantization using S1 with 1 fewer average bit
end




g1=(H*f1)';  % the optimal combiner given f1 is the transmit beamformer
g1=2*g1/max(abs(g1)); %normalize the coefficient to have a maximum magnitude of 2

g1less=(H*f1less)';% the optimal combiner given f1less is the transmit beamformer
g1less=2*g1less/max(abs(g1less)); 
g1less2=(H*f1less2)';  % the optimal combiner given f1less2 is the transmit beamformer
g1less2=2*g1less2/max(abs(g1less2)); 

g2=(H*f2)';
g2=2*g2/max(abs(g2))/1.000001; 

g2_less=(H*f2_less)';
g2_less=2*g2_less/max(abs(g2_less))/1.000001; 
g2_less2=(H*f2_less2)';
g2_less2=2*g2_less2/max(abs(g2_less2))/1.000001; 



for i=1:r 

    g1(i)=q1(g1(i),N_theta,N_phi);
    g1less(i)=q1(g1less(i),N_theta,N_phi/2);
    g2(i)=q2(g2(i),N_theta2,N_phi2);
    g1less2(i)=q1(g1less2(i),N_theta/2,N_phi/2);    
    g2_less(i)=q2(g2_less(i),N_theta,N_phi/2);
    g2_less2(i)=q2(g2_less2(i),N_theta,N_phi/4);

end


        
snr0=la;

snr1=abs(g1*H*f1)^2/abs(f1'*f1)/norm(g1)^2;  % the snr when f1 and g1 are used
snr1less=abs(g1less*H*f1less)^2/abs(f1less'*f1less)/norm(g1less)^2;
snr1less2=abs(g1less2*H*f1less2)^2/abs(f1less2'*f1less2)/norm(g1less2)^2;
snr2_ba=abs(g2*H*f2)^2/abs(f2'*f2)/norm(g2)^2;
snr2_less=abs(g2_less*H*f2_less)^2/abs(f2_less'*f2_less)/norm(g2_less)^2;
snr2_less2=abs(g2_less2*H*f2_less2)^2/abs(f2_less2'*f2_less2)/norm(g2_less2)^2;




ber0=ber0+qam_const*Q_fun(sqrt(pt'*snr0));
ber1=ber1+qam_const*Q_fun(sqrt(pt'*snr1));
ber1less=ber1less+qam_const*Q_fun(sqrt(pt'*snr1less));
ber1less2=ber1less2+qam_const*Q_fun(sqrt(pt'*snr1less2));
ber2=ber2+qam_const*Q_fun(sqrt(pt'*snr2_ba));
ber2_less=ber2_less+qam_const*Q_fun(sqrt(pt'*snr2_less));
ber2_less2=ber2_less2+qam_const*Q_fun(sqrt(pt'*snr2_less2));


end

end
%%  BER averaged over the number the iteration
ber0=ber0/itr;
ber1=ber1/itr;
ber1less=ber1less/itr;
ber1less2=ber1less2/itr;
ber2=ber2/itr;
ber2_less=ber2_less/itr;
ber2_less2=ber2_less2/itr;


%% 
figure; 
semilogy(snr_range, ber1less2,'b--d');hold on;
semilogy(snr_range, ber1less,'b--s');
semilogy(snr_range, ber1,'b--o');
semilogy( snr_range,ber2_less2,'r->');
semilogy( snr_range,ber2_less,'r-x');
semilogy( snr_range,ber2,'r-*');
semilogy(snr_range, ber0,'k--');

xlabel('P_t/\sigma_n^2');
ylabel('BER');
legend('S1, b=1','S1, b=1.5','S1, b=2','S2, b=1','S2, b=1.5','S2, b=2','unconstrained');
grid on
axis([-20 5 10^(-5) 1]);

