% This matlab code computes SNR loss as a function of bit rate b as in Fig 5.
% Yuan-Pei Lin, ypl@mail.nctu.edu.tw

clearvars;
close all;

cl=3;
ray=10;
path=cl*ray;
sigma_t=7.5*ones(1,cl)/180*pi;% standard deviation of phi_t
sigma_r=7.5*ones(1,cl)/180*pi;
t=16; %number of transmit antennas
r=16; %number of receive antennas
dt=1/2; %normalized transmit antenna spacing
dr=1/2; %normalized receive antenna spacing

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 phi_t
   sigma_te=7.5*ones(1,cl)/180*pi;% standard deviation of phi_t

   sigma_ra=7.5*ones(1,cl)/180*pi;
   sigma_re=7.5*ones(1,cl)/180*pi;

    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   




mu_itr=100;
itr=100;


af=0;
b=1:5; %average number of quantization bits
ratio1=zeros(size(b));

ratio2=ratio1;

ag1=ratio1;ag2=ratio1;

ratio1f=zeros(size(b));

ratio2f=ratio1f;
e=zeros(size(b));
eg=zeros(size(b));

snr0_ave=ratio1;
snr1_ave=ratio1;
snr2_ave=ratio1;
snr1f_ave=ratio1;

q=2.^(b);


rng(1970);
for n_mu=1:mu_itr
    mu_ta=2*pi*rand(cl,1);
    mu_ra=2*pi*rand(cl,1);
    mu_te=2*pi*rand(cl,1);
    mu_re=2*pi*rand(cl,1);
   
for n_itr=1:itr

   [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;
    af=af+norm(f0)^2;

    for n_b=1:size(b,2)
        N_theta=q(n_b);    N_phi=N_theta;
        N_theta2=N_theta*2; N_phi2=N_theta/2;

        f1=f0;
        for i=1:t  
            f1(i)=q1(f1(i),N_theta,N_phi);
            e(n_b)=e(n_b)+abs(f1(i)-f0(i))^2;
        end
        

        f2ba=f0/1.000001;

         for i=1:t    
             f2ba(i)=q2(f2ba(i),N_theta*2,N_phi/2);             
         end


        g0=(H*f1)';
        [~,i0]=max(abs(g0));
        g0=g0/g0(i0)*2;
        ag1(n_b)=ag1(n_b)+norm(g0)^2;

        g1=g0;
        for i=1:r 
            g1(i)=q1(g1(i),N_theta,N_phi);
            eg(n_b)=eg(n_b)+abs(g1(i)-g0(i))^2;
        end



            
        g2ba=(H*f2ba)';
        g2ba=g2ba/max(abs(g2ba))/1.0000001*2;              
        
        ag2(n_b)=ag2(n_b)+norm(g2ba)^2;
        for i=1:r          
            %scheme2, direct quantization, uniform

             g2ba(i)=q2(g2ba(i),N_theta*2,N_phi2/2);
                   
         end
snr0=la;
snr1=(abs(g1*H*f1))^2/(norm(f1))^2/(norm(g1))^2;

snr2ba=(abs(g2ba*H*f2ba))^2/(norm(f2ba))^2/(norm(g2ba))^2;

snr0_ave(n_b)=snr0_ave(n_b)+la;
snr1_ave(n_b)=snr1_ave(n_b)+snr1;

snr2_ave(n_b)=snr2_ave(n_b)+snr2ba;
snr1f_ave(n_b)=snr1f_ave(n_b)+(norm(H*f1))^2/(norm(f1))^2;

ratio1f(n_b)=ratio1f(n_b)+snr0/(norm(H*f1))^2*(norm(f1))^2;

ratio2f(n_b)=ratio2f(n_b)+snr0/(norm(H*f2ba))^2*(norm(f2ba))^2;
ratio1(n_b)=ratio1(n_b)+snr0/snr1;



ratio2(n_b)=ratio2(n_b)+snr0/snr2ba;

end

end

end


%%
af=af/itr/mu_itr/t;
ag1=ag1/itr/mu_itr/r;
ag2=ag2/itr/mu_itr/r;

ratio1=ratio1/itr/mu_itr;
ratio1f=ratio1f/itr/mu_itr;

ratio2f=ratio2f/itr/mu_itr;

ratio2=ratio2/itr/mu_itr;

e=e/itr/mu_itr/t;
eg=eg/itr/mu_itr/r;
snr0_ave=snr0_ave/itr/mu_itr;
snr1_ave=snr1_ave/itr/mu_itr;

snr2_ave=snr2_ave/itr/mu_itr;
snr1f_ave=snr1f_ave/itr/mu_itr;

%% 

crude1=(1+pi^2/3*2*2.^(-2*b)./af).*(1+pi^2/3*2*2.^(-2*b)./ag1);
crude1f=1+pi^2/3*2*2.^(-2*b)/af;

crude2f=1+(pi^2/3*2.^(-2*b).*af/4+1/3*2.^(-2*b)*4)./af;

crude2=1+(pi^2/3*2.^(-2*b).*af/4+1/3*2.^(-2*b)*4)./af+(pi^2/3*2.^(-2*b).*ag2/4+1/3*2.^(-2*b)*4)./ag2+(pi^2/3*2.^(-2*b).*af/4+1/3*2.^(-2*b)*4)./af.*(pi^2/3*2.^(-2*b).*ag2/4+1/3*2.^(-2*b)*4)./ag2;


%%
figure; 


plot(b,10*log10(ratio1),'b-*','Linewidth',1);hold on;
plot(b,10*log10(crude1),'b:*','Linewidth',1);

plot(b,10*log10(ratio1f),'k-+','Linewidth',1);
plot(b,10*log10(crude1f),'k:+','Linewidth',1);

plot(b,10*log10(ratio2),'m-o','Linewidth',1);
plot(b,10*log10(crude2),'m:o','Linewidth',1);
plot(b,10*log10(ratio2f),'r-d','Linewidth',1);

plot(b,10*log10(crude2f),'r:d','Linewidth',1);

legend('S1, E[D_{t,r}]','S1, est. of E[D_{t,r}]','S1, E[D_{t}]','S1, est. of E[D_{t}]','S2, E[D_{t,r}]','S2, est. of E[D_{t,r}]','S2, E[D_{t}]','S2, est. of E[D_{t}]');
xlabel('Average number of quantization bits b');
ylabel('SNR degradation (dB)')
set(gca,'xTick',b)
grid on;

