% This matlab code computes rate versus angular spread performance as in Fig 6(b).
% Yuan-Pei Lin, ypl@mail.nctu.edu.tw
%
% subscript notation
% 0: the optimal unconstrained case
% 1: quantized version using Structure 1
% 2: quantized version using Structure 2


clearvars;
close all;

t=64;  %number of transmit antennas
r=64; %number of receive antennas
M=2;


cl=3;
ray=10;
path=cl*ray;

ULA=0;% ULA=1 for ULA antenna array and ULA=0 for UPA antenna array
    
    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   


%ga=6.95; %snr gap for 10^(-5)
ga=-log(5*10^(-5))/1.5;



dt=1/2; %normalized transmit antenna spacing
dr=1/2; %normalized receive antenna spacing
kr=(0:r-1)';
kt=(0:t-1)';


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_phi/2;% the quantization resolution for Structure 2

snr=0;
pt=10.^(snr/10);

phi_itr=100;mu_itr=100;
itr=phi_itr*mu_itr;

sp=0:2:12; % range of angular spread 
r0=zeros(1,length(sp));r1=r0; r2=r0;
r_ay=r0;r_ay1=r0;

sinr0=zeros(M,length(pt));sinr1=sinr0; 
sinr2=sinr0;
sinr_ay=sinr0;sinr_ay1=sinr0;




for n_spread=1:length(sp)
    
   disp(strcat('spread=',num2str(sp(n_spread))));
   sigma_ta=sp(n_spread)*ones(1,cl)/180*pi;% standard deviation of phi_t
   sigma_te=sigma_ta;
   sigma_ra=sigma_ta;
   sigma_re=sigma_ta;
   
   rng(1970);
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);



   [U,La,V]=svd(H);
   F0=V(:,1:M);
   for i=1:M
       s=max(abs(F0(:,i)))/2; % for scaling f0, extra 2 s.t. max(abs(f0))=2
       F0(:,i)=F0(:,i)/s;
   end


   % quantization of Frf using structure 1
   Ph1=zeros(t,M);
   for i2=1:M
       Ph1(:,i2)=V(:,i2)-Ph1*pinv(Ph1)*V(:,i2);
       s=max(abs(Ph1(:,i2)))/2; % for scaling f0, extra 2 s.t. max(abs(f0))=2
       Ph1(:,i2)=Ph1(:,i2)/s;
        for i1=1:t 
            Ph1(i1,i2)=q1(Ph1(i1,i2),N_theta,N_phi);
 
        end
   end
   
%   quantization of Frf using structure 2

   
   Ph2=zeros(t,M);
   % quantize the subsequent vectors with orthogonal projection
   for i2=1:M

       Ph2(:,i2)=V(:,i2)-Ph2*pinv(Ph2)*V(:,i2);
       Ph2(:,i2)=Ph2(:,i2)/max(abs(Ph2(:,i2)))*2/1.0000001;

        for i1=1:t 

            Ph2(i1,i2)=q2(Ph2(i1,i2),N_theta*2,N_phi/2); 

        end
   end




 % power normalization 
 F0=V(:,1:M)/sqrt(M);
   
% Fbb computed using the following paper
%  A. Alkhateeb and R. W. Heath Jr, ?Frequency Selective Hybrid
% Precoding for Limited Feedback Millimeter Wave Systems,? 
% IEEE Transactions on Communications, vol. 64, no. 5, pp.1801-1818, May 2016.

   [~,~,Veff]=svd(H*Ph1);
   F1=Ph1*(Ph1'*Ph1)^(-1/2)*Veff;   
   F1=F1/norm(F1,'fro');%power=1
   
   
   [~,~,Veff]=svd(H*Ph2);
   F2=Ph2*(Ph2'*Ph2)^(-1/2)*Veff;   
   F2=F2/norm(F2,'fro');%power=1
   
    
   G0=U(:,1:M)';

   
   G1zf=pinv(H*F1);

   G2zf=pinv(H*F2);
   G2bazf=pinv(H*F2);

   
   Phg1=zeros(r,M);
   G1zf=G1zf';

   for i2=1:M
       Phg1(:,i2)=G1zf(:,i2)-Phg1*pinv(Phg1)*G1zf(:,i2);
       s=max(abs(Phg1(:,i2)))/2; % for scaling f0, extra 2 s.t. max(abs(f0))=2
       Phg1(:,i2)=Phg1(:,i2)/s;
 
        for i1=1:r
            %[~,i0]=min(abs(tab2-Phg1(i1,i2)));Phg1(i1,i2)=tab1(i0);
            Phg1(i1,i2)=q1(Phg1(i1,i2),N_theta,N_phi);
        end
   end     
   Phg1=Phg1';


   
   Phg2=zeros(r,M);
   G2bazf=G2bazf';
   for i2=1:M
       Phg2(:,i2)=G2bazf(:,i2)-Phg2*pinv(Phg2)*G2bazf(:,i2);
       Phg2(:,i2)=Phg2(:,i2)/max(abs(Phg2(:,i2)))*2/1.000001;
        for i1=1:r 
            Phg2(i1,i2)=q2(Phg2(i1,i2),N_theta*2,N_phi/2);
        end
   end     
   Phg2=Phg2';
   



   Ag1=(Phg1*H*F1)^(-1);

 
   Ag2=(Phg2*H*F2)^(-1);
 

   G1=Ag1*Phg1;
  
 
   G2=Ag2*Phg2;


    %comparison
    %O. E. Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. W. Heath Jr, 
    %?Spatially sparse precoding in millimeter wave MIMO systems,? 
    %IEEE Transactions on Wireless Communications, vol. 99, pp. 1-15, Jan. 2014.
    Nrf=2*M;
   Fres=F0;
   Frf=[];
   for i1=1:Nrf
       Si=At'*Fres;
       [~,i0]=max(diag(Si*Si'));
        a=exp(1i*round(angle(At(:,i0))/2/pi*N_theta)*2*pi/N_theta);
        Frf=[Frf a];
        Fbb=pinv(Frf)*F0;
        temp=Frf*Fbb;
        Fres=F0-temp;
        if i1==M
            F_ay=temp/norm(temp,'fro');% Nrf=M case
        end
           
   end
   F_ay1=temp/norm(temp,'fro');  %Nrf=2M case
   
   G_ay0=pinv(H*F_ay)';
   Wrf=[];
   Wres=G_ay0;
   for i1=1:M
       Si=Ar'*(H*F_ay*F_ay'*H')*Wres;
       [~,i0]=max(diag(Si*Si'));
       a=Ar(:,i0);
        a=exp(1i*round(angle(Ar(:,i0))/2/pi*N_theta)*2*pi/N_theta);
        Wrf=[Wrf a];
        Wbb=pinv(F_ay'*H'*Wrf);
%        Wbb=pinv(Wrf)*G_ay0;
        G_ay=Wrf*Wbb;
        Wres=G_ay0-G_ay;
   end;
   
   G_ay=G_ay';
   %G_ay=G_ay0';
   
   G_ay0=pinv(H*F_ay1)';
   Wrf=[];
   Wres=G_ay0;
   for i1=1:Nrf
       Si=Ar'*(H*F_ay1*F_ay1'*H')*Wres;
       [~,i0]=max(diag(Si*Si'));
       a=Ar(:,i0);
        a=exp(1i*round(angle(Ar(:,i0))/2/pi*N_theta)*2*pi/N_theta);
        Wrf=[Wrf a];
        Wbb=pinv(F_ay1'*H'*Wrf);
%        Wbb=pinv(Wrf)*G_ay0;
        G_ay1=Wrf*Wbb;
        Wres=G_ay0-G_ay1;
   end;
   
   G_ay1=G_ay1';
   
 
   
   
   
   %transfer function
   B0=G0*H*F0;
   B1=G1*H*F1;
   B2=G2*H*F2;
   B_ay=G_ay*H*F_ay;
      B_ay1=G_ay1*H*F_ay1;

  for i=1:M
      sinr0(i,:)=abs(B0(i,i))^2./(norm(B0(i,:),2)^2-abs(B0(i,i))^2+1./pt*norm(G0(i,:),2)^2);
      sinr1(i,:)=pt/norm(G1(i,:),2)^2;
      sinr2(i,:)=pt/norm(G2(i,:),2)^2;  
      sinr_ay(i,:)=abs(B_ay(i,i))^2./(norm(B_ay(i,:),2)^2-abs(B_ay(i,i))^2+1./pt*norm(G_ay(i,:),2)^2);
      sinr_ay1(i,:)=abs(B_ay1(i,i))^2./(norm(B_ay1(i,:),2)^2-abs(B_ay1(i,i))^2+1./pt*norm(G_ay1(i,:),2)^2);
   
  end

    r0(n_spread)=r0(n_spread)+sum(log2(1+sinr0/ga),1);
    r1(n_spread)=r1(n_spread)+sum(log2(1+sinr1/ga),1);
    r2(n_spread)=r2(n_spread)+sum(log2(1+sinr2/ga),1);
    r_ay(n_spread)=r_ay(n_spread)+sum(log2(1+sinr_ay/ga),1);
    r_ay1(n_spread)=r_ay1(n_spread)+sum(log2(1+sinr_ay1/ga),1);
end

end
end
% average sum rate over the number of realizations
r0=r0/itr;
r1=r1/itr;
r2=r2/itr;
r_ay=r_ay/itr;
r_ay1=r_ay1/itr;



%%
figure

plot(sp, r0,'k--');hold
plot(sp,r2,'r-*');
plot(sp, r1,'b--o');
plot(sp,r_ay1,'m-.^');
plot(sp,r_ay,'g-.d');
legend('unconstrained','S2','S1','SPC, N_{rf}=4','SPC, N_{rf}=2');
ylabel('Transmission rate');xlabel('Angular Spread (degrees)');
axis([0 12 8 15]);
grid on;
