scm_core.m

3GPP推荐的SCM信道模型的MATLAB源码

        d_s = antpar.BsElementPosition * wavelength;
        
        X_BS_v = reshape(BsGain(:,:,1,:,:),K,S,N,M); 
        X_BS_h = reshape(BsGain(:,:,2,:,:),K,S,N,M); 
        X_MS_v = reshape(MsGain(:,:,1,:,:),K,U,N,M); 
        X_MS_h = reshape(MsGain(:,:,2,:,:),K,U,N,M);     
        aod = bulkpar.aods(1:K,1:N,1:M)*pi/180;
        aoa = bulkpar.aoas(1:K,1:N,1:M)*pi/180;  
        phase_v_v = reshape(bulkpar.subpath_phases(1:K,1,1:N,1:M),K,N,M)*pi/180;
        phase_v_h = reshape(bulkpar.subpath_phases(1:K,2,1:N,1:M),K,N,M)*pi/180;
        phase_h_v = reshape(bulkpar.subpath_phases(1:K,3,1:N,1:M),K,N,M)*pi/180;
        phase_h_h = reshape(bulkpar.subpath_phases(1:K,4,1:N,1:M),K,N,M)*pi/180;
        %sq_r_n1 = reshape(sqrt(1/bulkpar.xpd(1:K,1,1:N)),K,N);
        %sq_r_n2 = reshape(sqrt(1/bulkpar.xpd(1:K,2,1:N)),K,N);
        r_n1 = 1 ./ reshape(bulkpar.xpd(1:K,1,1:N),K,N);
        r_n2 = 1 ./ reshape(bulkpar.xpd(1:K,2,1:N),K,N);    
        v = linkpar.MsVelocity(1:K);
        theta_v = linkpar.MsDirection(1:K)*pi/180;
        sq_Pn = sqrt(path_powers_all(1:K,1:N*L));
        XpdIndependentPower_mex = 0;
        if strcmpi(scmpar.XpdIndependentPower,'yes')
            XpdIndependentPower_mex = 1;
        end
        % the ANSI-C function call
        [H temp_output_SubPathPhases] = scm_mex_core(POLARIZED, X_BS_v, X_BS_h, X_MS_v, X_MS_h, aod, aoa, d_s, d_u, phase_v_v, phase_v_h, phase_h_v, phase_h_h, r_n1, r_n2, t, k_CONST, v, theta_v, sq_Pn, look_up_points, U, S, N, M, K, T, GainsAreScalar, XpdIndependentPower_mex);
        
        output_SubPathPhases(:,1,:,:) = temp_output_SubPathPhases;
        output_SubPathPhases(:,2,:,:) = temp_output_SubPathPhases;
        output_SubPathPhases(:,3,:,:) = temp_output_SubPathPhases;
        output_SubPathPhases(:,4,:,:) = temp_output_SubPathPhases;
        
        output_SubPathPhases = prin_value((output_SubPathPhases*180/pi + bulkpar.subpath_phases));
        
        
        
    end % is it polarized?
    
    %%%%%%
    if PROFILE_MODE_FLAG
        profile report
    end
    %%%%%%
    
    % LOS OPTION
    
    if strcmp(lower(scmpar.ScmOptions),'los')
        
        
        if DISPLAY_MODE_FLAG
            disp('entering LOS option...');
        end
        
        % Take the values of K factors and probability of having LOS case
        K_factors       = bulkpar.K_factors;
        
        %indx            = (K_factors.'~=0)';
        
        ThetaBs      = linkpar.ThetaBs; ThetaBs=ThetaBs(:).';
        ThetaMs      = linkpar.ThetaMs; ThetaMs=ThetaMs(:).';
        
        % if strcmp(str,'safe') % we only do it 'safe' for this option
        
        %adjusting parameters for calling the C-language routine
        output_Phi_LOS  = zeros(K,1);
        d_u = antpar.MsElementPosition * wavelength;
        d_s = antpar.BsElementPosition * wavelength;
        
        
        % the ANSI-C function call
        [H output_Phi_LOS] = scm_mex_core(LOS, BsGain_Theta_BS, MsGain_Theta_MS, ThetaBs*pi/180, ThetaMs*pi/180, d_s, d_u, bulkpar.Phi_LOS(:,1)* pi/180, t, k_CONST, linkpar.MsVelocity(:), linkpar.MsDirection(1,:)*pi/180, H, output_Phi_LOS, K_factors, U, S, N*L, K, T);
        
        % adjusting angles
        output_Phi_LOS = prin_value((output_Phi_LOS*180/pi + bulkpar.Phi_LOS));
        
    end % if 'LOS' option is activated   
    %%%%%%%%%%%%%%%%%%ANSI-C core part ends%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    
else % if ANSI-C is not used
    if ~strcmpi(scmpar.ScmOptions,'polarized')
        if DISPLAY_MODE_FLAG
            disp('entering main loop...');
        end
        
        output_SubPathPhases = zeros(K,N,M);
        
        for u = 1:U % cycles (MS) antennas
            %du = antpar.MsElementSpacingULA * (u-1) * wavelength;
            du = antpar.MsElementPosition(u)*wavelength;
            for s = 1:S % cycles Tx (BS) atennas
                %ds = antpar.BsElementSpacingULA * (s-1) * wavelength;
                ds = antpar.BsElementPosition(s)*wavelength;
                for k = 1:K % cycles links
                    for n = 1:N % cycles paths
                        
                        LM_index = 0; % 
                        
                        for km = 1:L % cycles midpaths
                            
                            temp                 = zeros(M,T);   % oversized, just to keep it always the same size
                            
                            for m=1:LN(km) % cycles subpaths
                                
                                LM_index = LM_index+1;
                                
                                temp(m,:) =  BsGain(k,s,n,LM(LM_index)) *...
                                    exp(j * (...
                                    k_CONST * ds * sin((bulkpar.aods(k,n,LM(LM_index)))*pi/180) +...
                                (bulkpar.subpath_phases(k,n,LM(LM_index))*pi/180)+...
                                    k_CONST * du * sin((bulkpar.aoas(k,n,LM(LM_index)))*pi/180)...
                                )) *...
                                    MsGain(k,u,n,LM(LM_index)) * ...
                                    exp(j * k_CONST * linkpar.MsVelocity(k) * cos((bulkpar.aoas(k,n,LM(LM_index)) - linkpar.MsDirection(k))*pi/180) * t(k,:));
                                
                            end % subpaths
                            
                            H(u,s,(n-1)*L+km,:,k) =  sqrt(path_powers_all(k,(n-1)*L+km) / LN(km)) * sum(temp,1);
                            
                        end % midpaths
                        
                    end % paths 
                    
                end % links 
            end % Tx antennas 
        end % Rx antennas 
        
        
        for k = 1:K % cycles links  % Of course the for loop could be avoided. 
            for n = 1:N % cycles paths
                for m=1:M % cycles supaths
                    
                    % SIMPLE LOOP
                    output_SubPathPhases(k,n,m) =  k_CONST * linkpar.MsVelocity(k) * cos((bulkpar.aoas(k,n,m) - linkpar.MsDirection(k))*pi/180) * (delta_t(k)*T);   
                    
                end % subpaths
            end % paths 
        end % links 
        
        output_SubPathPhases = prin_value((output_SubPathPhases*180/pi + bulkpar.subpath_phases));
        
    else % it's polarized!
        if DISPLAY_MODE_FLAG
            disp('entering polarized option...');
        end    
        
        output_SubPathPhases = zeros(K,4,N,M);
        
        
        % BsGain has must have size: [K x S x 2 x N x M]
        % the first dimension in the polarization must be vertical
        %
        % MsGain has must have size: [K x U x 2 x N x M]
        % the first dimension in the polarization must be vertical
        %
        % subpath_phases has size: [K x 4 x N x M]
        % bulkpar.xpd has size: [K,2,N]
        temp = zeros(M,T);    
        for u = 1:U % cycles (MS) antennas
            %du = antpar.MsElementSpacingULA * (u-1) * wavelength;
            du = antpar.MsElementPosition(u)*wavelength;
            for s = 1:S % cycles Tx (BS) atennas
                %ds = antpar.BsElementSpacingULA * (s-1) * wavelength;
                ds = antpar.BsElementPosition(s)*wavelength;
                for k = 1:K % cycles links
                    for n = 1:N % cycles paths
                        
                        LM_index = 0; % 
                        
                        for km = 1:L % cycles midpaths
                            
                            temp                 = zeros(M,T);   % oversized, just to keep it always the same size
                            
                            for m=1:LN(km) % cycles subpaths
                                
                                LM_index = LM_index+1;
                                
                                if strcmpi(scmpar.XpdIndependentPower,'yes')
                                
                                    temp(m,:) =  [BsGain(k,s,1,n,LM(LM_index)) BsGain(k,s,2,n,LM(LM_index))] * ...
                                        [sqrt(bulkpar.xpd(k,1,n)/(1+bulkpar.xpd(k,1,n))) *exp(j*bulkpar.subpath_phases(k,1,n,LM(LM_index))*pi/180)     , sqrt(1/(1+bulkpar.xpd(k,1,n))) * exp(j*bulkpar.subpath_phases(k,2,n,LM(LM_index))*pi/180);...
                                            sqrt(1/(1+bulkpar.xpd(k,2,n))) * exp(j*bulkpar.subpath_phases(k,3,n,LM(LM_index))*pi/180),      sqrt(bulkpar.xpd(k,2,n)/(1+bulkpar.xpd(k,2,n))) * exp(j*bulkpar.subpath_phases(k,4,n,LM(LM_index))*pi/180)] *...
                                        [MsGain(k,u,1,n,LM(LM_index)); MsGain(k,u,2,n,LM(LM_index))] * ...
                                        exp(j * (k_CONST * ds * sin((bulkpar.aods(k,n,LM(LM_index)))*pi/180))) * ...
                                        exp(j * (k_CONST * du * sin((bulkpar.aoas(k,n,LM(LM_index)))*pi/180))) * ...                    
                                        exp(j * k_CONST * linkpar.MsVelocity(k) * cos((bulkpar.aoas(k,n,LM(LM_index)) - linkpar.MsDirection(k))*pi/180) * t(k,:));
                                    
                                else 
                                    
                                    temp(m,:) =  [BsGain(k,s,1,n,LM(LM_index)) BsGain(k,s,2,n,LM(LM_index))] * ...
                                        [ exp(j*bulkpar.subpath_phases(k,1,n,LM(LM_index))*pi/180)     , sqrt(1/bulkpar.xpd(k,1,n)) * exp(j*bulkpar.subpath_phases(k,2,n,LM(LM_index))*pi/180);...
                                            sqrt(1/bulkpar.xpd(k,2,n)) * exp(j*bulkpar.subpath_phases(k,3,n,LM(LM_index))*pi/180),      exp(j*bulkpar.subpath_phases(k,4,n,LM(LM_index))*pi/180)] *...
                                        [MsGain(k,u,1,n,LM(LM_index)); MsGain(k,u,2,n,LM(LM_index))] * ...
                                        exp(j * (k_CONST * ds * sin((bulkpar.aods(k,n,LM(LM_index)))*pi/180))) * ...
                                        exp(j * (k_CONST * du * sin((bulkpar.aoas(k,n,LM(LM_index)))*pi/180))) * ...                    
                                        exp(j * k_CONST * linkpar.MsVelocity(k) * cos((bulkpar.aoas(k,n,LM(LM_index)) - linkpar.MsDirection(k))*pi/180) * t(k,:));
                                    
                                end
                                
                            end % subpaths
                            
                            H(u,s,(n-1)*L+km,:,k) =  sqrt(path_powers_all(k,(n-1)*L+km) / LN(km)) * sum(temp,1);
                            
                        end % midpaths
                    end % paths 
                end % links 
            end % Tx antennas 
        end % Rx antennas 
        
        for k = 1:K % cycles links  % Of course the for loop could be avoided. 
            for n = 1:N % cycles paths
                for m=1:M % cycles supaths
                    
                    
                    output_SubPathPhases(k,:,n,m) =  (k_CONST * linkpar.MsVelocity(k) * cos((bulkpar.aoas(k,n,m) - linkpar.MsDirection(k))*pi/180) * (t(k,end)+delta_t(k))) * ones(1,4);
                    
                end % subpaths
            end % paths 
        end % links 
        
        output_SubPathPhases = prin_value((output_SubPathPhases*180/pi + bulkpar.subpath_phases));
        
        
    end % is it polarized?
    
    
    
    %%%%%%
    if PROFILE_MODE_FLAG
        profile report
    end
    %%%%%%
    
    
    % LOS OPTION
    
    if strcmpi(scmpar.ScmOptions,'los')
        
%         if ~ strcmpi(scmpar.Scenario,'urban_micro')
%             error('LOS option is possible only for URBAN MICRO scenario') % assumption-all the users are in the same scenario
%         end
        if DISPLAY_MODE_FLAG
            disp('entering LOS option...');
        end
        
        % Take the values of K factors and probability of having LOS case
        K_factors       = bulkpar.K_factors;
        
        indx            = (K_factors.'~=0)';
        
        ThetaBs      = linkpar.ThetaBs; ThetaBs=ThetaBs(:).';
        ThetaMs      = linkpar.ThetaMs; ThetaMs=ThetaMs(:).';
        
        
        %temp = zeros(M,T);    
        
        
        
        for k = 1:K % cycles links
            if indx(k)~=0
                
                for u = 1:U % cycles (MS) antennas
                    %du = antpar.MsElementSpacingULA * (u-1) * wavelength;
                    du = antpar.MsElementPosition(u)*wavelength;
                    for s = 1:S % cycles (BS) atennas
                        %ds = antpar.BsElementSpacingULA * (s-1) * wavelength;
                        ds = antpar.BsElementPosition(s)*wavelength;
                        temp =  BsGain_Theta_BS(k,s) * exp(j * k_CONST * ds * sin( ThetaBs(k)*pi/180))* ...
                            MsGain_Theta_MS(k,u) * exp(j * (k_CONST * du * sin( ThetaMs(k)*pi/180  ) + bulkpar.Phi_LOS(k,1) * pi/180 )) * ...
                            exp(j * k_CONST * linkpar.MsVelocity(k) * cos((ThetaMs(k) - linkpar.MsDirection(k))*pi/180) * t(k,:));
                        
                        % all the parameters are fixed within one drop
                        
                        H(u,s,1,:,k) = (sqrt(1/(K_factors(k)+1)) * squeeze(H(u,s,1,:,k)) + sqrt(K_factors(k)/(K_factors(k)+1)) * temp.').';
                        
                    end % Rx antennas
                end % Tx antennas
                
                H(:,:,2:end,:,k)= sqrt(1/(K_factors(k)+1)) *  H(:,:,2:end,:,k);    
                
            end % if there is LOS
            
        end % links
        
        output_Phi_LOS       = zeros(K,1);
        
        for k = 1:K % cycles links
            if indx(k)~=0      
                output_Phi_LOS(k,1) = k_CONST * linkpar.MsVelocity(k) * cos((ThetaMs(k) - linkpar.MsDirection(k))*pi/180) * (t(k,end)+delta_t(k));                    
            end % if there is LOS        
        end % links 
        
        output_Phi_LOS = prin_value((output_Phi_LOS*180/pi + bulkpar.Phi_LOS));
        
    end % if 'LOS' option is activated
    
end

%%%%%%%%%%%%%%%%
%%%%%%%%%%%
%%%%%%%%
%%%%%        That's all folks !!!
%%
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A function that maps inputs from (-inf,inf) to (-180,180)
function y=prin_value(x)
y=mod(x,360);
y=y-360*floor(y/180);