self_deco.m

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

function [SINR,BER] = self_deco(soi_code, PSC, P, M, ind_sin, ind_in, delay_spread, L, sensors, snr_input, sir_input, theta, BS, DD, TrialAll,main_groupusr,minor_groupusr,fix_sir)

  %additional variable declaration:
  mSINRin = 0;
  mSINRw  = 0;
  mSINRwh = 0;
  mSINRwd = 0;
  miSINRin = 0;
  miSINRw  = 0;
  miSINRwh = 0;
  miSINRwd = 0;
  total_error = zeros(3,1);
  
rand('state',sum(100*clock));
randn('state',sum(100*clock));

num_sensor = length(sensors); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TO DEFINE THE INCIDENT SIGNALS' MULTIPATH PARAMETERS:
%
[SrcAll, PathAll] = size(theta);
amp(1:SrcAll,1)  = [1   10^(-1*sir_input/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_input/20)) / sqrt(2); 

%Nc = length(code(1,:))/L; % chips per information symbol (# of chips in signature sequence)
Nc = length(soi_code)/L; % chips per information symbol (# of chips in signature sequence)

%h_mf1 = conj(code(1,Nc*L:-1:1));
h_mf1 = conj(soi_code(1,Nc*L:-1:1));

  
% each trial is one statistically independent Monte Carlo experiment
%ww = waitbar(0, ['SNR=', int2str(snr_input), 'dB,     SIR=', int2str(sir_input),'dB']);
for trial = 1:TrialAll
  clear signal sig d_PSC d_PS code;
  % waitbar(trial/TrialAll, ww);    
  
  % to generate interfering sources' PN codes (they are randomized from run to run for
  % each scenario)
  %
  code(1,:) = soi_code ; 
  rand('state',sum(100*clock));
  % the program code below are correct for L=1,2,3,4 ... 
  for src = 2:SrcAll
     temp_code = sign(rand(1,Nc)-0.5);  
     for i=1:Nc
       code(src,(i-1)*L+1:i*L) = temp_code(1,i)*ones(1,L); 
     end;
  end;
  %

  % Arrival time delay between direct path and 2nd path :
  % Unlike original codes: the delay is generally different for each multipath
  delay(1,1)         = 0;
  %delay(1,2:PathAll) = round(rand(1,PathAll-1) * (delay_spread) - 0.5); 
  delay(1,2:PathAll)  = floor(rand(1,PathAll-1) * (delay_spread+1)); 
  %delay(2:SrcAll,:)  = round(rand(SrcAll-1,PathAll) * (Nc*L) - 0.5);
  delay(2:SrcAll,:)   = floor(rand(SrcAll-1,PathAll) * (Nc*L));
  % to generate the incident signals' information symbol sequences:
  data = round(rand(SrcAll,M+1))*2 - 1;  
  for i = 1 : SrcAll
     signal = data(i,1)*code(i,:);
     for m = 2:M+1     % the (M+1)-th symbol is actually the 0th symbol
        signal = [signal      (data(i,m)*ones(size(code(i,:)))) .* code(i,:)];      
     end;
     sig(i,:) = signal;
  end;

  % Assumes SrcAll = 2 in the following block of MATLAB codes:
  % Unlike the original codes: the (M+1)-th symbol is unequal to the 1st symbol :     
  for i = 1 : SrcAll
    sig_m(i,delay(i,2)+1 : Nc*M*L)     = sig(i,1                       : Nc*M*L-delay(i,2));
    sig_m(i,1            : delay(i,2)) = sig(i,Nc*(M+1)*L-delay(i,2)+1 : Nc*(M+1)*L);
  end;
 
  % Unlike the original codes, frequency-nonselective time-invariant fading is randomized
  % for each multipath in each Monte Carlo run.
  % rho_mp = 0.5*exp(j*pi/4);   % needs to be randomized for each trial and each multipath
  c = amp .* exp(j * 2 * pi * rand(SrcAll,PathAll));
  c(1,1) = abs(c(1,1));
  
  % to form I+N part of data set :
  y_IN = n_amp * (randn(num_sensor,Nc*M*L) + j*randn(num_sensor,Nc*M*L));
  for i = 2 : SrcAll
    y_IN = y_IN + c(i,1)*exp(j*2*pi*sin(theta(i,1))*sensors) * sig(i,1:M*Nc*L); 
    y_IN = y_IN + c(i,2)*exp(j*2*pi*sin(theta(i,2))*sensors) * sig_m(i,:);  
  end;
    
  % to form S part of data set :
  y_S = c(1,1)*exp(j*2*pi*sin(theta(1,1))*sensors) * sig(1,1:M*Nc*L); 
  y_S = y_S + c(1,2)*exp(j*2*pi*sin(theta(1,2))*sensors) * sig_m(1,:);  
  
  x = y_S + y_IN;           % SOI + I + N = received baseband data

  % To pre-process received baseband data with self-decorrelator before despreading:
  d_PSC = [];
  % d_PS  = [];
  for m = 1:M
    % self-decorrelator to remove SOI from I+N segment
    d_PSC = [d_PSC, x(:,(m-1)*L*Nc+1:m*L*Nc) * PSC.']; 
  end;
 
  % despreading by single-user-type detector match-filtering:
  for sensor = 1:num_sensor
     mfx(sensor,:)     = conv(x(sensor,:),    h_mf1);
     mfy_S(sensor,:)   = conv(y_S(sensor,:),  h_mf1);
     mfy_IN(sensor,:)  = conv(y_IN(sensor,:), h_mf1);
     mfd_PSC(sensor,:) = conv(d_PSC(sensor,:),h_mf1);
  end;
												
  if BS == 1
     % to transform from element-space to reduced-rank beamspace:
     %[T_czr, T_sd, T_upb] = beam_transform_M(x,y_IN,M,L,Nc,PS,code1); 
     
     % CZR's original algorithm w/o self-decorrelator:
     %TT = T_czr;	% for Z&R algorithm   
     %[wy1, wd1, w, X_stack, Xs_stack]...
     %  = rake2dw(M, T_czr, L, Nc, mfx, mfy_IN, mfd_PSC, DD, delay_spread); 
  
     % CZR's algorithm w/ self-decorrelator:
     %TT = T_sd;  
     %[wy1, wd1, w, X_stack, Xs_stack]...
     %  = rake2dw(M,    TT, L, Nc, mfx, mfy_IN, mfd_PSC, DD, delay_spread);   
  
     % CZR's algorithm w/o self-decorrelator --- hypothetically w/o SOI in I+N delay segment:
     %TT = T_upb;		% for No SOI in K_{I+N};
     %[w, wy1, wd1, mfx_sin_stack, mfy_IN_sin_stack]...
     %  = rake2dw(M,    TT, L, Nc, mfx, mfy_IN, mfd_PSC, DD, delay_spread);  
  else
   [w, wh, wd]...
    = rake2d(M, eye(num_sensor), P, delay_spread, ind_sin, ind_in, mfx, mfy_IN, mfd_PSC, DD);
  end;
  
  % for computation of beamformers' output SINR"
  mfy_S_relevant   = mfy_S(:, ind_sin); % used later in place of mfy_S_all
  mfy_IN_relevant  = mfy_IN(:,ind_sin); % used later in place of mfy_IN_all
  mfy_S_all_stack  = [];
  mfy_IN_all_stack = [];
  for p = 1:M
     ind_s = (p-1)*delay_spread+1;
     ind_e = p    *delay_spread;
     mfy_S_all_stack   = [mfy_S_all_stack  stack(mfy_S_relevant(:,ind_s:ind_e))];
     mfy_IN_all_stack  = [mfy_IN_all_stack stack(mfy_IN_relevant(:,ind_s:ind_e))];
  end;

  temp_S = 0;
  temp_IN=0;
  for in=1:size(mfy_S_all_stack,2)
     temp_S = temp_S + (norm(mfy_S_all_stack(:,in)))^2;
     temp_IN=temp_IN+(norm(mfy_IN_all_stack(:,in)))^2;
  end;
%instanteous value of SINR
  SINRin = temp_S/temp_IN;
  SINRw  = (norm(w' * mfy_S_all_stack) / norm(w' * mfy_IN_all_stack))^2;
  SINRwh = (norm(wh' * mfy_S_all_stack) / norm(wh'* mfy_IN_all_stack))^2;
  SINRwd = (norm(wd' * mfy_S_all_stack) / norm(wd'* mfy_IN_all_stack))^2;
  
 %adaptively calculate the mean of SINR 
  mSINRin = ((trial-1)/(trial))*mSINRin + SINRin/(trial);
  mSINRw  = ((trial-1)/(trial))*mSINRw  + SINRw/(trial);
  mSINRwh = ((trial-1)/(trial))*mSINRwh + SINRwh/(trial);
  mSINRwd = ((trial-1)/(trial))*mSINRwd + SINRwd/(trial);
 
  %adaptively calculate the mean of 1/SINR 
  miSINRin = ((trial-1)/(trial))*miSINRin + (1/SINRin)/(trial);
  miSINRw  = ((trial-1)/(trial))*miSINRw  + (1/SINRw) /(trial);
  miSINRwh = ((trial-1)/(trial))*miSINRwh + (1/SINRwh)/(trial);
  miSINRwd = ((trial-1)/(trial))*miSINRwd + (1/SINRwd)/(trial);
 
  RakeSig=[w wh wd]'* (mfy_S_all_stack+mfy_IN_all_stack);
%  RakeSig2=[w wh wd]'* mfy_S_all_stack;    
  BitEst=sign(real(RakeSig(:,2:end)./RakeSig(:,1:end-1)));
%  BitEst2=sign(real(RakeSig2(:,2:end)./RakeSig2(:,1:end-1)));
  BitTheor=ones(3,1)*(data(1,2:end-1).*data(1,1:end-2));

  total_error = total_error + (-sum(BitEst==BitTheor,2)+M-1);
  %BER(:,trial)=(-sum(BitEst==BitTheor,2)+M-1)/(M-1);         PT version
  
end;
clear  w  wh  wd  mfy_S_all  mfy_IN_all  mfy_S_all_stack  mfy_IN_all_stack ;
%close(ww);   % to close waitbar
SINR  = [mSINRw;  mSINRwh;  mSINRwd;  mSINRin; 1/miSINRw;  1/miSINRwh;  1/miSINRwd;  1/miSINRin];
%BER=mean(BER,2);                                             PT version
BER = total_error./((M-1)*trial);