How to solve this '''Input to EIG must not contain NaN or Inf'''. problem
This is the code: % 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
Prashant Kumar answered .
2025-11-20
num=1;
if num==1
range=[0:3:15]; % range of angular spread in deg
Therefore range will start with 0.
for nn=1:length(range)
if num==1
deg=range(nn); % angular spread in deg
num is still 1, so the angular speed case applies. range(nn) will be accessed, and as indicated earlier that starts with 0. So deg will be assigned 0.
sig_ta=deg*ones(1,path)/180*pi;
sig_te=deg*ones(1,path)/180*pi;
with deg being 0, this makes sig_ta an sig_te all 0.
b=stat_Rh(sector,ULA,t,r,cl,ray,dt,tH,tW,alpha,mu_ta,mu_te,sig_ta,sig_te,Ta,Te);
temp=temp*exp_lapsec(-(m-i)*dt*sin(mute(ll))+(n-k)*dt*cos(mute(ll))*sin(muta(ll)),sig_te(ll),Te(1)-mute(ll),Te(2)-mute(ll));
temp=temp*exp_lapsec((n-k)*dt*sin(mute(ll))*cos(muta(ll)),sig_ta(ll),Ta(1)-muta(ll),Ta(2)-muta(ll));
function y= exp_lapsec(c,s,tmin,tmax)
so those 0 will be received into the variable s
b0=s/sqrt(2);
Since s is 0, b0 will be 0.
if tmin>=0
y=1/c0/(2*b0)/(1i*c-1/b0)*(exp((1i*c-1/b0)*tmax)-exp((1i*c-1/b0)*tmin));
elseif tmax<=0
y=1/c0/(2*b0)/(1i*c+1/b0)*(exp((1i*c+1/b0)*tmax)-exp((1i*c+1/b0)*tmin));
else
y=1/c0/(2*b0)*((1-exp((1i*c+1/b0)*tmin))/(1i*c+1/b0)+(exp((1i*c-1/b0)*tmax)-1)/(1i*c-1/b0));
end