Prosenjit Paul asked . 2023-10-19

How to solve this '''Input to EIG must not contain NaN or Inf'''. problem

% 0: optimal BF
% 0q: optimal BF


clearvars;
close all;
filename='stat';
 seed=1970;
 rng(seed);
% rng(107);
cl=3;
ray=10;
deg=5;
path=cl*ray;
t=16;  %number of transmit antennas
r=16; %number of receive antennas
dt=1/2*2*pi; %normalized transmit antenna spacing
dr=1/2*2*pi; %normalized receive antenna spacing
phi_itr=100;mu_itr=100;  itr=phi_itr*mu_itr;
ULA=0; %1 for ULA and 0 for UPA

sector=4
if sector==1
    Ta=[0 2*pi];Te=[0 2*pi];
elseif sector==2
    Ta=[0 pi];
elseif sector==3
    Ta=[pi/6 5*pi/6];
elseif sector==4
    Ta=[0 2*pi];Te=[0 2*pi];
elseif sector==5
    Ta=[pi/3 2*pi/3]; Te=[pi/2, pi/2+pi/4];
end


num=1;
if num==1
    
   range=[0:3:15]; % range of angular spread in deg
  %range=3; % range of angular spread in deg
   disp('~~SNR vs angular spread');
   
elseif num==2
    range=1:8; % number of cluster
    disp('~~SNR vs number of cluster');
elseif num==3
    range=1:10; % number of ray per cluster
    disp('~~SNR vs number of rays (per cluster)');
elseif num==4
    range=2:7; %number of antenna
    disp('~~SNR vs number of antennas');
    
elseif num==5
    range=1; %
    disp('~~BER vs Pt/N0');
    snr_range=-25:2:5; 
    %snr_range=0:2:25; 
    qam=4;
    pt=10.^(snr_range/10);
    pt=pt*3/(qam-1);qam_const=1;
    ber0=zeros(size(snr_range));ber0q=ber0;ber1=ber0;ber1w=ber0;ber1q=ber0;ber_ay=ber0;ber_ay1=ber0;

elseif num==6
    range=0;
    disp('~~capacity~~');
    pt=10.^(range/10);
    cp0=zeros(1,itr);cp0q=cp0;
    cp1=cp0;cp1q=cp0;
    cp_ay=cp0;cp_ay1=cp0;

elseif num==7
    range=0;
        disp('~~transmission rate~~');
        snr_range=-25:2:5;
        pt=10.^(snr_range/10);
        Rb0=zeros(size(snr_range));Rb0q=Rb0;Rb1=Rb0;Rb2=Rb0;Rb1q=Rb0;Rb1s=Rb0; Rb_ay=Rb0;Rb_ay_bf=Rb0;
end   
Nth=8; Bf=4;  Nrf=1;
Rh=zeros(t,t);asnr0=zeros(size(range)); asnr0q=asnr0; asnr1=asnr0;asnr2=asnr0; asnr2q=asnr0; asnr10=asnr0; asnr1q=asnr0;asnr1a=asnr0;asnr1b=asnr0;

if ULA==1
    tW=1;
    tH=t;
else
  if ceil(sqrt(t))==sqrt(t)       
       tW=sqrt(t);
       tH=tW;
  else
        tW=sqrt(2*t);
        tH=tW/2;
  end  
end
if ceil(sqrt(r))==sqrt(r)        
       rW=sqrt(r);
       rH=rW;
else
        rW=sqrt(2*r);
        rH=rW/2;
end 


index=0;
for nn=1:length(range)
    
    if num==1   
       deg=range(nn);  % angular spread in deg  
       disp(strcat('~~~~~angular spread=',num2str(deg)));
    elseif num==2
        cl=range(nn); % number of cluster
    elseif num==3
        ray=range(nn); % number of ray per cluster
    elseif num==4
        t=2^range(nn); %number of antenna
        r=t;
        if ceil(sqrt(t))==sqrt(t)       
           tW=sqrt(t);tH=tW;
        else
           tW=sqrt(2*t);tH=tW/2;
        end      
        if ceil(sqrt(r))==sqrt(r)        
           rW=sqrt(r);rH=rW;
        else
            rW=sqrt(2*r);rH=rW/2;
        end 
    end   
        
        sigma_ta=deg*ones(1,cl)/180*pi;% standard deviation of phi_t
       sigma_te=sigma_ta;% standard deviation of phi_t   
       sigma_ra=sigma_ta;
       sigma_re=sigma_ta;
 
        sig_ta=deg*ones(1,path)/180*pi;
        sig_te=deg*ones(1,path)/180*pi;

          
    for n_mu=1:mu_itr
        
        mu_ta=Ta(1)+(Ta(2)-Ta(1))*rand(cl,1);
        mu_ra=2*pi*rand(cl,1);
        mu_te=Te(1)+(Te(2)-Te(1))*rand(cl,1);
        mu_re=2*pi*rand(cl,1);

 

       alpha=ones(1,path); 

 
 
       b=stat_Rh(sector,ULA,t,r,cl,ray,dt,tH,tW,alpha,mu_ta,mu_te,sig_ta,sig_te,Ta,Te);
 
 
 
 
 [V,D]=eig(b);
        [la1,i0]=max(diag(D));
        f1=V(:,i0);
        
 
        if ULA
            At=exp(1i*dt*(0:t-1)'*cos(mu_ta)');
            [f10,~]=eigs(At*At',1);
         end
            A=zeros(t,t);
            f2=zeros(t,1);
            for i_cl=1:cl
                temp_t=exp(1i*dt*(0:tH-1)'*cos(mu_te(i_cl)))*exp(1i*dt*(0:tW-1)*sin(mu_te(i_cl))*sin(mu_ta(i_cl)));
                temp_t=reshape(temp_t,t,1);
                f2=f2+temp_t;
                A=A+temp_t*temp_t';
            end
            [f10,~]=eigs(A,1);
            f2=f2/norm(f2);
    
       
    
        

        asnr1a(nn)=asnr1a(nn)+la1;
        Rh=zeros(t,t);
        for n_itr=1:phi_itr
            
            index=index+1;
            
            if sector==1
               [H, At, Ar]=stat_mmCh(ULA,dr,dt,t,r,tW,tH,rW,rH, path, cl,ray,mu_ta, mu_te, mu_ra, mu_re, sigma_ta,sigma_te,sigma_ra,sigma_re);
            else
               [H, At, Ar,spc_phi_ta,spc_phi_te]=stat_mmCh_sector(Ta,Te,ULA,dr,dt,t,r,tW,tH,rW,rH, path, cl,ray,mu_ta, mu_te, mu_ra, mu_re, sigma_ta,sigma_te,sigma_ra,sigma_re);
            end
            [V,D]=eig(H'*H);
            [la,i0]=max(diag(D));
            f0=V(:,i0);
            
            Rh=Rh+H'*H;           
            
            snr0=la;
%            f0q=exp(1i*round(angle(f0)/2/pi*Nth)*2*pi/Nth);
             f0=f0/max(abs(f0))*2;  f0q=q1(f0,Nth,Nth);
            g0q=H*f0q;
 %           g0q=exp(1i*round(angle(g0q)/2/pi*Nth)*2*pi/Nth);
            g0q=q1(g0q/max(abs(g0q))*2,Nth,Nth);
            
            
            f1q=q1(f1/max(abs(f1))*2,Nth,Nth);
%             f1q=exp(1i*round(angle(f1)/2/pi*Nth)*2*pi/Nth);
            g1q=H*f1q;
%             g1q=exp(1i*round(angle(g1q)/2/pi*Nth)*2*pi/Nth);
 
            g1q=q1(g1q/max(abs(g1q))*2,Nth,Nth);
            snr0q=norm(g0q'*H*f0q,2)^2/norm(g0q,2)^2/norm(f0q,2)^2; %using quanized optimal bf

            snr1=norm(H*f1,2)^2/norm(f1,2)^2;;
            snr1q=norm(g1q'*H*f1q,2)^2/norm(g1q,2)^2/norm(f1q,2)^2; %quantized

            snr10=norm(H*f10,2)^2;


            snr2=norm(H*f2,2)^2;
            

    
    if num<=4
        asnr0(nn)=asnr0(nn)+snr0;
        asnr0q(nn)=asnr0q(nn)+snr0q;
        asnr1(nn)=asnr1(nn)+snr1;
        asnr2(nn)=asnr2(nn)+snr2;
        asnr10(nn)=asnr10(nn)+snr10;
        asnr1q(nn)=asnr1q(nn)+snr1q;

    end
    
 

if num==5    
    
    [snr_ay,snr_ay1]= Ayach(1,t,dt,tH,tW,Nrf,H,ULA,At,Ar,f0,Ta,Te, 4, Nth,spc_phi_ta,spc_phi_te);
    [snr_ay_bf,snr_ay1]= Ayach(1,t,dt,tH,tW,Nrf,H,ULA,At,Ar,f0,Ta,Te, 6, Nth,spc_phi_ta,spc_phi_te);
    ber0=ber0+qam_const*Q_fun(sqrt(pt*snr0));
    ber0q=ber0q+qam_const*Q_fun(sqrt(pt*snr0q));
    ber1=ber1+qam_const*Q_fun(sqrt(pt*snr1));
     ber2=ber2+qam_const*Q_fun(sqrt(pt*snr2));
    %   ber1w=ber1w+qam_const*Q_fun(sqrt(pt*snr1w));
    ber1q=ber1q+qam_const*Q_fun(sqrt(pt*snr1q));
    %  ber1s=ber1s+qam_const*Q_fun(sqrt(pt*snr1s));
    %   ber1sq=ber1sq+qam_const*Q_fun(sqrt(pt*snr1sq));

    ber_ay=ber_ay+qam_const*Q_fun(sqrt(pt*snr_ay));
    ber_ay1=ber_ay1+qam_const*Q_fun(sqrt(pt*snr_ay1));
end
 
if num==6
    [snr_ay,snr_ay1]= Ayach(1,t,dt,tH,tW,Nrf,H,ULA,At,Ar,f0,Ta,Te, Bf, Nth,spc_phi_ta,spc_phi_te);

    cp0(index)=snr0;
    cp0q(index)=snr0q;
    cp1(index)=snr1;
    cp1q(index)=snr1q;
    cp_ay(index)=snr_ay;
    cp_ay1(index)=snr_ay1;

end
    
if num==7  
    [snr_ay,snr_ay1]= Ayach(1,t,dt,tH,tW,Nrf,H,ULA,At,Ar,f0,Ta,Te, Bf, Nth,spc_phi_ta,spc_phi_te);
    [snr_ay_bf,snr_ay1]= Ayach(1,t,dt,tH,tW,Nrf,H,ULA,At,Ar,f0,Ta,Te, Bf+1, Nth,spc_phi_ta,spc_phi_te);

    
    Rb0=Rb0+log2(1+pt*snr0);
    Rb0q=Rb0q+log2(1+pt*snr0q);
    Rb1=Rb1+log2(1+pt*snr1);
    Rb2=Rb2+log2(1+pt*snr2);
    Rb1q=Rb1q+log2(1+pt*snr1q);
    Rb_ay=Rb_ay+log2(1+pt*snr_ay);
    Rb_ay_bf=Rb_ay_bf+log2(1+pt*snr_ay_bf);
    Rb1s=Rb1s+log2(1+pt*la1);
 end





        end

        Rh=Rh/phi_itr;
        temp=eig(Rh);
        asnr1b(nn)=asnr1b(nn)+max(temp);

    end


    

end





%%
if num<=4
    asnr0=10*log10(asnr0/itr);
    asnr0q=10*log10(asnr0q/itr);
    asnr1=10*log10(asnr1/itr);
    asnr1q=10*log10(asnr1q/itr);
    asnr2=10*log10(asnr2/itr);
    asnr2q=10*log10(asnr2q/itr);
    asnr10=10*log10(asnr10/itr);
    asnr1b=10*log10(asnr1b/mu_itr);
    asnr1a=10*log10(asnr1a/mu_itr);
    
    
    %% spawc
    figure;
    
plot(range, asnr0,'k--o'); hold on;
%plot(range, asnr0q,'k->');
%plot(range, asnr1b,'b-s');
plot(range, asnr1,'r--d');
%plot(range, asnr1q,'r--*');
%plot(range, asnr2q,'g--+');



%legend('full CSI','full CSI,quantized','opt stat','stat','stat,quantized');

if num==1    
   xlabel('Angular spread (deg)');% range of angular spread in deg   
elseif num==2
    xlabel('Number of clusters');% number of cluster
elseif num==3
    xlabel('Number of rays per cluster'); % number of ray per cluster
else
    xlabel('Number of antennas'); %number of antenna
end   
%axis([0 15 18 25]); %for t=16
axis([0 15 22 32]); %t=64, r=16
%axis([0 15 28 35]); %t=r=64
ylabel('Average SNR (dB)');
grid on;
%%
figure;
    
plot(range, asnr0,'k'); hold on;
plot(range, asnr1,'g--d');
plot(range, asnr2,'r--+');
plot(range, asnr10,'g--o');
plot(range, asnr1a,'r-*');
plot(range, asnr1b,'b-s');
legend('optimal','statistical','multibeam','stat,zero var','theoretical','opt stat');

if num==1    
   xlabel('Angular spread (deg)');% range of angular spread in deg   
elseif num==2
    xlabel('Number of clusters');% number of cluster
elseif num==3
    xlabel('Number of rays per cluster'); % number of ray per cluster
else
    xlabel('Number of antennas'); %number of antenna
end   


ylabel('Average SNR (dB)');
grid on;

     temp=clock;
    st=strcat(num2str(temp(2)),'.',num2str(temp(3)),'.',num2str(temp(4)),'.',num2str(temp(5)));   
    st=strcat(st,',stat-new.m, BER, t=',num2str(t),',r=',num2str(r));
    st=strcat(st,',cl=',num2str(cl));
    st=strcat(st,',ray=',num2str(ray(1)));
    st=strcat(st,',deg=',num2str(sigma_ta(1)/pi*180));
    st1=strcat('Ta=',num2str(Ta/pi),'\pi');
    st1=strcat(st1,',Te=',num2str(Te/pi),'\pi');
 

    st1=strcat(st1,',mu_itr=',num2str(mu_itr),',phi_itr=',num2str(phi_itr));
   if ULA==1        
       st1=strcat(st1,',ULA');
    else
        st1=strcat(st1,',UPA');
   end
   st1=strcat(st1,',Bf=',num2str(Bf),',Nth=',num2str(Nth));
   title({st,st1});
    grid on

end


%%
if num==4
    figure;
    range=2.^range;
plot(range, asnr0,'k-.o'); hold on;
plot(range, asnr0q,'k->'); 
plot(range, asnr1,'r-.d');
plot(range, asnr1q,'r-*');
legend('FULL CSI','FULL CSI,quantized','stat','stat,quantized');
xlabel('Number of antennas'); 
ylabel('Average SNR (dB)');
grid on;

end


%%
if num==5
   ber0=ber0/itr; 
   ber0q=ber0q/itr; 
    ber1=ber1/itr;
    ber2=ber2/itr;
    ber1q=ber1q/itr;
    ber_ay=ber_ay/itr;
    ber_ay1=ber_ay1/itr;
  %%
  figure;
    semilogy(snr_range, ber0, 'k-');hold on;
    semilogy(snr_range, ber0q, 'k-*');
    semilogy(snr_range, ber1, 'r-s');
    semilogy(snr_range, ber1q, 'r-*');
    semilogy(snr_range, ber_ay, 'g-d');
    axis([min(snr_range) max(snr_range) 10^(-5) 1]);
    
    legend('0','0q','1','1q','ay');
     temp=clock;
    st=strcat(num2str(temp(2)),'.',num2str(temp(3)),'.',num2str(temp(4)),'.',num2str(temp(5)));   
    st=strcat(st,',stat-new.m, BER, t=',num2str(t),',r=',num2str(r));
    st=strcat(st,',cl=',num2str(cl));
    st=strcat(st,',ray=',num2str(ray(1)));
    st=strcat(st,',deg=',num2str(sigma_ta(1)/pi*180));
    st1=strcat('Ta=',num2str(Ta/pi),'\pi');
    st1=strcat(st1,',Te=',num2str(Te/pi),'\pi');
 

    st1=strcat(st1,',mu_itr=',num2str(mu_itr),',phi_itr=',num2str(phi_itr));
   if ULA==1        
       st1=strcat(st1,',ULA');
    else
        st1=strcat(st1,',UPA');
   end
   st1=strcat(st1,',Bf=',num2str(Bf),',Nth=',num2str(Nth));
   title({st,st1});
    grid on
    
    
end

%% plot cdf

if num==6

    cp0=log2(1+cp0*pt);
    cp0q=log2(1+cp0q*pt);
    cp1=log2(1+cp1*pt);
    cp1q=log2(1+cp1q*pt);
    cp_ay=log2(1+cp_ay*pt);
    cp_ay1=log2(1+cp_ay1*pt);
    
    [pdf0,x0]=hist(cp0,100);
    [pdf0q,x0q]=hist(cp0q,100);
    [pdf1,x1]=hist(cp1,100);
    [pdf1q,x1q]=hist(cp1q,100);
    
    [pdf_ay,x_ay]=hist(cp_ay,100);
    [pdf_ay1,x_ay1]=hist(cp_ay1,100);
    pdf0=pdf0/itr;
    pdf0q=pdf0q/itr;
    pdf1=pdf1/itr;
    pdf1q=pdf1q/itr;
    
    pdf_ay=pdf_ay/itr;
    pdf_ay1=pdf_ay1/itr;
    
    
    cdf0=zeros(size(pdf0));cdf0q=cdf0;cdf1=cdf0; cdf1q=cdf0; cdf_ay=cdf0;cdf_ay1=cdf0;
    cdf0(1)=pdf0(1);
    cdf0q(1)=pdf0q(1);
    cdf1(1)=pdf1(1);
    cdf1q(1)=pdf1q(1);
    cdf_ay(1)=pdf_ay(1);cdf_ay1(1)=pdf_ay1(1);
    for ii=2:100
        cdf0(ii)=cdf0(ii-1)+pdf0(ii-1);
        cdf0q(ii)=cdf0q(ii-1)+pdf0q(ii-1);
        cdf1(ii)=cdf1(ii-1)+pdf1(ii-1);
        cdf1q(ii)=cdf1q(ii-1)+pdf1q(ii-1);
        cdf_ay(ii)=cdf_ay(ii-1)+pdf_ay(ii-1);
        cdf_ay1(ii)=cdf_ay1(ii-1)+pdf_ay1(ii-1);
    end
    %%
    figure;
    
    plot(x0,cdf0,'k'); hold;
    plot(x0q,cdf0q,'k:');
    plot(x1,cdf1,'r');
    plot(x1q,cdf1q,'--');
    plot(x_ay,cdf_ay,'-.m');
    plot(x_ay1,cdf_ay1,'--g');
    legend('cdf0','cdf0q','cdf1','cdf1q','SPC, N_{rf}=1','SPC, N_{rf}=2','Location','NorthWest');
    
    temp=clock;
    st=strcat(num2str(temp(2)),'.',num2str(temp(3)),'.',num2str(temp(4)),'.',num2str(temp(5)));   
    st=strcat(st,filename,'cdf, t=',num2str(t),',r=',num2str(r));
    st=strcat(st,',cl=',num2str(cl));
   % st=strcat(st,',ray=',num2str(ray));
    st=strcat(st,',deg=',num2str(sigma_ta(1)/pi*180));
    st=strcat(st,',pt(dB)=',num2str(range));
    st=strcat(st,',Bf=',num2str(Bf));
    st=strcat(st,',Nth=',num2str(Nth));
    st=strcat(st,',Ta=',num2str(Ta/pi),'\pi');
    st=strcat(st,',Te=',num2str(Te/pi),'\pi');
    if ULA==1        
       st1=strcat('ULA,','itr=');
    else
        st1=strcat('UPA,','itr=');
    end

    st1=strcat(st1,num2str(itr));
    title({st,st1});    

    grid on
    axis([2 10 0 1.1]);
    grid on;
    
end
if num==7
    Rb0=Rb0/itr; 
    Rb0q=Rb0q/itr;
    Rb1s=Rb1s/itr;
    Rb1=Rb1/itr;
    Rb2=Rb2/itr;
    Rb1q=Rb1q/itr;
    Rb_ay=Rb_ay/itr;
    Rb_ay_bf=Rb_ay_bf/itr;
 
    %%
    figure;
    plot(snr_range, Rb0, 'k-.');hold on;
    plot(snr_range, Rb0q, 'k-*');
    plot(snr_range, Rb1, 'r-d');
    plot(snr_range, Rb2, 'r-+');
    plot(snr_range, Rb1q, 'r-s');
    plot(snr_range, Rb1s, 'b-x')
    plot(snr_range, Rb_ay, 'm->');

    legend('optimal','0q','1','2multibeam','1q','1s','ay');
    %axis([min(snr_range) max(snr_range) 10^(]);
    
     temp=clock;
    st=strcat(num2str(temp(2)),'.',num2str(temp(3)),'.',num2str(temp(4)),'.',num2str(temp(5)));   
    st=strcat(st,'stat-new.m, Rb, t=',num2str(t),',r=',num2str(r));
    st=strcat(st,',cl=',num2str(cl));
    st=strcat(st,',ray=',num2str(ray));
    st=strcat(st,',deg=',num2str(sigma_ta(1)/pi*180));
    if ULA==1        
       st1=strcat('ULA,','itr=');
    else
        st1=strcat('UPA,','itr=');
    end

    st1=strcat(st1,num2str(itr));
    title({st,st1});    
    grid on
    
     %%
    figure;
    plot(snr_range, Rb0, 'k-.o');hold on;
    plot(snr_range, Rb0q, 'k->');
plot(snr_range, Rb_ay_bf, 'm:o');   
    plot(snr_range, Rb1, 'r-.d');
    plot(snr_range, Rb1q, 'r-*');
        
     plot(snr_range, Rb_ay, 'b-+');
    legend('Full CSI','Full CSI, quantized','SPC, 10 bits','stat','stat, quantized','SPC, 8bits');
    axis([-20 0 1 10]) ;
        grid on
        ylabel('Transmission rate (bits)');
        xlabel('P_t/N_0');

end

Radar , Phased Array System Toolbox , Phased Array Design and Analysis