lscorrelation.m

建模源程序,程序比较经典很实用.后面还有很多经典的程序.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generation of correlated large scale (LS) parameters DS,ASA,ASD
%  and SF for all links
%
%  If layout parameters (co-ordinates etc.) are given, auto-correlation
%  between channel segments is generated. If layout parameters are not
%  defined, only cross-correlation between LS parameters is generated.
%
%  Notes! Auto-correlation is generated for each BS separately. If two MSs
%  are linked to same BS, correlation is distance dependent. If one MS is
%  linked to two sectors of single BS, correlation is full.
%
% The procedure and notations are from [1]
% Ref: [1]
% Authors: Pekka Ky鰏ti (EBIT), Per Zetterberg (KTH)

% Bug fixes:
%   Lines with comment %PZ are modified by Per Zetterberg 28.4.2006

function sigmas = LScorrelation(wimpar,linkpar,fixpar)

% scenario parameters
Scenario = wimpar.Scenario;
PropagCondition = wimpar.PropagCondition;
DS_lambda = fixpar.DS_lambda; 
AS_D_lambda = fixpar.AS_D_lambda; 
AS_A_lambda = fixpar.AS_A_lambda;
SF_lambda = fixpar.SF_lambda;


if isfield(linkpar,'BsXY')  % if layout parameters given, generate auto-correlation
    % layout parameters
    BsXY = linkpar.BsXY;
    NofSect = linkpar.NofSect;
    BsOmega = linkpar.BsOmega;
    MsXY = linkpar.MsXY;
    MsOmega = linkpar.MsOmega;
    Pairing = linkpar.Pairing;
    
    % determine MS/BS/Sect indices from Pairing matrix
    K = sum(Pairing(:));    % number of links
    [r c] = find(Pairing);
    indMs = c';      % links coupling to MSs
    indSect = r';    % links coupling to Sectors
    % get links coupling to BSs
    tmpsum = cumsum(NofSect);
    for k=1:length(r)
        tmp = find(diff(r(k)<=tmpsum));
        if isempty(tmp)
            indBs(k)=1;
        else
            indBs(k)=tmp+1;
        end
    end
    
    
    % generate iid 4xK Gaussian random numbers
    ksi = randn(4,K);
    
    %%% Auto-correlation generation %%%
    % determine auto-correlation separately for each BS
    for k=1:max(indBs)
        indLink = find(indBs==k);   % indices to MSs linked to BS k
        tmpMs = indMs(indLink);     % indices to MSs linked to BS k
        if length(tmpMs)>1  % if only 1 MS linked to BS, no auto-correlation
            dist = sqrt((repmat(MsXY(1,tmpMs),length(tmpMs),1)- ...
                repmat(MsXY(1,tmpMs)',1,length(tmpMs))).^2+ ...
                (repmat(MsXY(2,tmpMs),length(tmpMs),1)- ...
                repmat(MsXY(2,tmpMs)',1,length(tmpMs))).^2);
           
            % Note! below real(sqrtm(exp(-dist/DS_lambda))) fullfills C=MM' with 
            %  symmetric C. real(.) is needed to clean up tiny imaginary
            %  component
            %  which results from numeric operations
            warning('off', 'MATLAB:sqrtm:SingularMatrix');
            % DS auto-correlation
            M = real(sqrtm(exp(-dist/DS_lambda))); 
            ksi(1,indLink) = (M*ksi(1,indLink)').';
            % ASD auto-correlation
            M = real(sqrtm(exp(-dist/AS_D_lambda))); 
            ksi(2,indLink) = (M*ksi(2,indLink)').';
            % ASA auto-correlation
            M = real(sqrtm(exp(-dist/AS_A_lambda))); 
            ksi(3,indLink) = (M*ksi(3,indLink)').';
            % SF auto-correlation
            M = real(sqrtm(exp(-dist/SF_lambda))); 
            ksi(4,indLink) = (M*ksi(4,indLink)').';
            warning('on', 'MATLAB:sqrtm:SingularMatrix');
        end % end if
        
    end % end for
    
else  % if layout parameters not given
    
    K = size(linkpar.MsBsDistance,2);   % number of links
    % generate non-correlated 4xK Gaussian random numbers
    ksi = randn(4,K);
    
end % if (layout parameters given)


%%% Cross-correlation generation %%%
% Extract cross correlation parameters from input
a = fixpar.asD_ds;          % departure AS vs delay spread
b = fixpar.asA_ds;          % arrival AS vs delay spread
c = fixpar.asA_sf;          % arrival AS vs shadowing std
d = fixpar.asD_sf;          % departure AS vs shadoving std
e = fixpar.ds_sf;           % delay spread vs shadoving std
f = fixpar.asD_asA;         % departure AS vs arrival AS

% Cross-correlation matrix R(0) from [1, eq. 6.6]
% Order of rows and columns is ds asD asA sf
A = [ 1  a  b  e  ;...
      a  1  f  d  ;...
      b  f  1  c  ;...
      e  d  c  1  ];
  
% get R^(0.5)(0) for [1, eq. 6.6]
[V,D] = eig(A);
R_sqrt = V*sqrt(D);
ksi = R_sqrt*ksi;


if ((upper(Scenario)=='A1') | (upper(Scenario)=='C1') | (upper(Scenario)=='C2') | (upper(Scenario)=='D1') | (upper(Scenario)=='B3'))%PZ
    variances = [fixpar.DS_sigma, fixpar.AS_D_sigma, fixpar.AS_A_sigma, fixpar.SF_sigma]; %PZ
    Av=diag(variances)*A*diag(variances);%PZ
    [Vv,Dv] = eig(Av);%PZ
    [dummy,order_eig]=sort(diag(Dv));%PZ
    Vv=Vv(:,order_eig);%PZ
    Dv=diag(Dv);%PZ
    Dv=diag(Dv(order_eig));%PZ
    R_sqrtv = Vv*sqrt(Dv);%PZ
    ksiv = R_sqrtv*ksi;%PZ
end;

if (upper(Scenario)=='B1') % PZ
    Av=diag([1,1,1,fixpar.SF_sigma])*A*diag([1,1,1,fixpar.SF_sigma]);%PZ
    [Vv,Dv] = eig(Av);%PZ
    [dummy,order_eig]=sort(diag(Dv));%PZ
    Vv=Vv(:,order_eig);%PZ
    Dv=diag(Dv);%PZ
    Dv=diag(Dv(order_eig));%PZ
    R_sqrtv = Vv*sqrt(Dv);%PZ
    ksiv = R_sqrtv*ksi;%PZ
end;



% transform Normal distributed random numbers to scenario specific distributions
switch upper(Scenario)      % See distributions in [1, table 3.1]
    
    %% A1 should have six LS paramaters
    %% Where are the D2 parameters coming from, can't find them in D5.4 v.1.4 ?
    
    case {'A1','C1','C2','D1','D2'}     % in A1,C1,C2 and D1 LOS/NLOS distributions are equal
        sigma_ds  = 10.^(ksiv(1,:).'   + fixpar.DS_mu);      % LN PZ
        sigma_asD = 10.^(ksiv(2,:).' + fixpar.AS_D_mu);      % LN PZ
        sigma_asA = 10.^(ksiv(3,:).' + fixpar.AS_A_mu);      % LN PZ
        sigma_sf  = 10.^(0.1*ksiv(4,:).');                   % LN dB PZ
        
    case {'B1'}
        % Note1! Q(x)=0.5*efrc(-x/sqrt(2)), see Matlab helpwin erf.
        % Note2! 10^(x) = exp(log(10)x) in [1, p.18]
        if strcmpi(PropagCondition,'los')        
            sigma_ds  = 10.^(fixpar.DS_sigma*log(log(1./(0.5*erfc(ksiv(1,:)'/sqrt(2)))))...
                        +fixpar.DS_mu);                                      % Gumb.  PZ
            sigma_asD = 10.^(fixpar.AS_D_sigma*log(1./(1-...
                        (0.5*erfc(ksiv(2,:)'/sqrt(2))))-1)+fixpar.AS_D_mu);  % Logist. PZ
            sigma_asA = 10.^(fixpar.AS_A_sigma*log(1./(1-...
                        (0.5*erfc(ksiv(3,:)'/sqrt(2))))-1)+fixpar.AS_A_mu);  % Logist. PZ
            sigma_sf  = 10.^(0.1*ksiv(4,:).');                               % LN dB. PZ
            
        elseif strcmpi(PropagCondition,'nlos')
            sigma_ds  = 10.^(fixpar.DS_sigma*log(log(1./(0.5*erfc(ksiv(1,:)'/sqrt(2)))))...
                        +fixpar.DS_mu);                                      % Gumb PZ
            sigma_asD = 10.^(fixpar.AS_D_sigma*log(log(1./(0.5*erfc(ksiv(2,:)'/sqrt(2)))))...
                        +fixpar.AS_D_mu);                                    % Gumb PZ
            sigma_asA = 10.^(fixpar.AS_A_sigma*log(log(1./(0.5*erfc(ksiv(3,:)'/sqrt(2)))))...
                        +fixpar.AS_A_mu);                                    % Gumb PZ
            sigma_sf  = 10.^(0.1*ksiv(4,:).');                               % LN dB PZ
            
        else
            error(['Propagation condition LOS/NLOS must be defined for scenario ' Scenario '.'])
        end     % end if

    case {'B3'}     % in B3 distributions are equal in LOS and NLOS
        sigma_ds  = ksiv(1,:).'   + fixpar.DS_mu;                      % Normal PZ
        sigma_asD = ksiv(2,:).' + fixpar.AS_D_mu;                      % Normal PZ
        sigma_asA = ksiv(3,:).' + fixpar.AS_A_mu;                      % Normal PZ
        sigma_sf  = 10.^(0.1*fixpar.SF_sigma*ksiv(4,:).');             % LN dB PZ
        
end % end switch

% output
sigmas = [sigma_asD sigma_asA sigma_ds sigma_sf];