analysis_3dec2003.m

智能天线中的波束成型算法

function SINRth = analysis(Nc, Ns, P, M, delay_spread, L, sensors, snr, sir, theta, main_groupusr,minor_groupusr,fix_sir)

num_sensor = length(sensors); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TO DEFINE THE INCIDENT SIGNALS' MULTIPATH PARAMETERS:
%
[SrcAll, PathAll] = size(theta);
amp(1:SrcAll,1)  = [1   10^(-1*sir/20) * ones(1,main_groupusr-1) 10^(-1*fix_sir/20) * ones(1,minor_groupusr)].';  
% the time-delayed multipaths' amplitudes relative to the direct paths'
amp(:,2) = 0.5 * amp(:,1);

% the Input SNR of the SOI's direct path 
n_amp = amp(1,1) * (10.^(-1*snr/20)) / sqrt(2); 

C_M=kron(eye(Nc),ones(L,L));

% assuming that Ns=L*Nc

first_row= Nc*[L:-1:1 zeros(1,delay_spread-L)];
R_L=toeplitz(first_row',first_row);

%R_M=R_L; %%% this is correct only for L=1.....

Nl=delay_spread;
%F =eye(Nl);              %%%%% computing theoretical covariance matrix of SOI
                         %%%%% in the I+N part
                         %%% this is correct only for L=1.....
%J=[zeros(1,Nl); eye(Nl-1) zeros(Nl-1,1)];
%for ind=1:Nl-1
%  F=[F F(:,(ind-1)*Nl+1:ind*Nl)*J'; J*F((ind-1)*Nl+1:ind*Nl,:) eye(Nl)];
%end
%F=F*Ns; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% computation of R_M and F via Monte Carlo simulation
%
Npok=200;
n=floor(Ns/2);
R_M=zeros(Nl,Nl);  F=zeros(Nl^2,Nl^2);
for ind=1:Npok
  c=kron(sign(rand(1,Nc)-0.5),ones(1,L));
  Cl=[zeros(1,n-1) c(1,1:Ns-n+1)];
  Clb=[c(1,Ns-n+2:Ns) zeros(1,Ns-n+1)];
  Cl0=c;
  Cl0b=zeros(1,Ns);
  for k=1:Nl-1
    Cl=[Cl; zeros(1,n+k-1) c(1,1:Ns-k-n+1)];
    Clb=[Clb; c(1,Ns-n-k+2:Ns) zeros(1,Ns-n-k+1)];
    Cl0=[Cl0; zeros(1,k) c(1,1:Ns-k)];
    Cl0b=[Cl0b; c(1,Ns-k+1:Ns) zeros(1,Ns-k)];
  end  
  Rm=Cl*C_M*Cl'+Clb*C_M*Clb';
  u=stack(Cl0b*Cl');
  v=stack(Cl0b*Clb'+Cl0*Cl');
  w=stack(Clb*Cl0');
  Fi=u*u'+v*v'+w*w';
  R_M=R_M+Rm;
  F=F+Fi;
end
R_M=R_M/Npok;
F=F/Npok;
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Em=zeros(1,Nl);
Em(1,1)=1;  Em(2,delay_spread)=1;  %%%%% Auxiliary matrix E_M

  amp(1:SrcAll,1)  = [1 10^(-1*sir/20) * ones(1,main_groupusr-1) 10^(-1*fix_sir/20) * ones(1,minor_groupusr)].';  
  amp(:,2) = 0.5 * amp(:,1);
  n_amp = amp(1,1) * (10.^(-1*snr/20)) / sqrt(2);
  A_S = [exp(j*2*pi*sin(theta(1,1))*sensors)*amp(1,1) ...
         exp(j*2*pi*sin(theta(1,2))*sensors)*amp(1,2)];
  A_I = [exp(j*2*pi*sensors*sin(theta(2:SrcAll,1))').*(ones(num_sensor,1)*amp(2:SrcAll,1)') ...
         exp(j*2*pi*sensors*sin(theta(2:SrcAll,2))').*(ones(num_sensor,1)*amp(2:SrcAll,2)')];
  R_I = A_I*A_I';
  R_S = A_S*A_S';
  R_N = 2*n_amp^2*eye(num_sensor);
  
  K_IN=kron(R_M,R_I)+kron(R_L,R_N);
  AsE=A_S*Em;
  pom=kron(eye(Nl),AsE);
  K_S=pom*F*pom';
  
  d_0=stack(A_S*Em*R_L);
  w_0= (K_IN+K_S)\d_0;
  
  psi=sum(diag(A_S'*inv(R_S+R_I+R_N)*A_S));
  N_AF=floor(Nc * L /delay_spread)-1; % P-fold redundance in I+N delay segemnt
  
  MSINR_input_theor  = d_0'*d_0/sum(diag(K_IN)); 
  MSINR_out_WLLZ_theor = (w_0'*d_0)^2/(w_0'*K_IN*w_0);  
  MSINR_out_SD_theor = d_0'*(K_IN\d_0);  %%% = theory for "free", also
  MSINR_out_WLLZ_nonasym  = Nc^2*psi^2/...
    (Nc*psi*(1+num_sensor*Nl/(M*N_AF))+num_sensor*Nl/M*(3+num_sensor*Nl/M));
                 %%%%%%%%%%%% ANALYSIS-2 expression - non-asymptotic
  SINRth=real([MSINR_out_WLLZ_theor MSINR_out_SD_theor...
              MSINR_out_WLLZ_nonasym MSINR_input_theor]);