correlation.h

测试信道相关性所需要的一些关于统计和矩阵运算的头文件

/* **************************************************************************************** */
/* correlation.h ,C++ header file used in correlation_gen.gcc                               */
/* 		                                                                                    */
/* This file computes the correlation matrices RTx at Tx side                               */
/* and RRx at Rx side according to the PAS type(uniform,                                    */
/* truncated Gaussian and truncated Laplacian)                                              */
/*                                                                                          */
/* The following functions are included in this file:                                       */
/*                                                                                          */
/* void normalisation_uniform(int cluster_number,float power_lin[],float AS_deg[],float *Q) */
/* void normalisation_gaussian(int cluster_number,float power_lin[],float AS_deg[],         */
/*                             float delta_phi_deg[], float *Q, float *sigma_deg)           */
/* void normalisation_laplacian(int cluster_number,float power_lin[],float AS_deg[],        */
/*                              float delta_phi_deg[], float *Q, float *sigma_deg)          */
/* cluster_number ,  the number of clusters                                                 */
/* power_lin[] ,  the power of clusters in linear values                                    */
/* AS_deg[] ,  the Angular Spread of each cluster                                           */
/* delta_phi_deg ,  the half domain definition of the PAS,being either 90 or 180 degrees    */
/* *Q ,  the normalization coeffecient of PAS                                               */
/* *sigma_deg[] ,  used later to derive Rxx_*() and Rxy_*()                                 */ 
/* The above functions are used to derive the normalization coefficients of the PAS.        */
/*                                                                                          */
/* float Rxx_uniform(float d_norm,float phi_0_deg,float AS_deg )                            */
/* float Rxx_gaussian(float d_norm,float phi_0_deg,float sigma_deg,float delta_phi_deg)     */
/* float Rxx_laplacian(float d_norm,float phi_0_deg,float sigma_deg,float delta_phi_deg)    */
/* float Rxy_uniform(float d_norm,float phi_0_deg,float AS_deg )                            */
/* float Rxy_gaussian(float d_norm,float phi_0_deg,float sigma_deg,float delta_phi_deg)     */
/* float Rxy_laplacian(float d_norm,float phi_0_deg,float sigma_deg,float delta_phi_deg)    */
/* d_norm ,  the spacing between antennas                                                   */ 
/* phi_0_deg ,  the AoA(Angle-of-Arrival)                                                   */
/* AS_deg ,  the Angular Spread                                                             */
/* sigma_deg ,  result returned by the normalisation_*() functions                          */
/* delta_phi-deg ,  the half domain definition of the PAS,being either 90 or 180 degrees    */                                                                 */
/* The functions Rxx_*() of the above are used to compute the homogeneous part of the       */
/* correlation function(between  real or imaginary part)                                    */
/* The functions Rxy_*() of the above are used to compute the cross part of the             */
/* correlation function(between  real and imaginary part)                                   */
/*                                                                                          */
/* float **uniform_real(int N,float spacing,int cluster_number,float amplitude_cluster[],   */
/*                      float phi_deg[],float AS_deg[],float delta_phi_deg[])               */
/* float **gaussian_real(int N,float spacing,int cluster_number,float amplitude_cluster[],  */
/*                       float phi_deg[],float AS_deg[],float delta_phi_deg[])              */
/* float **laplacian_real(int N,float spacing,int cluster_number,float amplitude_cluster[], */
/*                       float phi_deg[],float AS_deg[],float delta_phi_deg[])              */
/* complex **uniform_cmpx(int N,float spacing,int cluster_number,float amplitude_cluster[], */
/*                      float phi_deg[],float AS_deg[],float delta_phi_deg[])               */
/* complex **gaussian_cmpx(int N,float spacing,int cluster_number,float amplitude_cluster[],*/
/*                          float phi_deg[],float AS_deg[],float delta_phi_deg[])           */
/* complex **laplacian_cmpx(int N,float spacing,int cluster_number,float amplitude_cluster[]*/
/*                           ,float phi_deg[],float AS_deg[],float delta_phi_deg[])         */
/* N ,  the number of antennas                                                              */
/* spacing ,  the normalized spacing between antennas ,in wavelength                        */
/* cluster_number ,  the number of clusters                                                 */
/* amplitude_cluster[] ,  the power of each cluster in linear values                        */
/* phi_deg[] ,  AoAs at Rx side or AoDs at Tx side, in degrees                              */
/* AS_deg[] ,  the Angular Spread of each cluster                                           */
/* delta_phi_deg , the half domain definition of the PAS,being either 90 or 180 degrees     */
/* The *_real functions of the above are used to compute the  correlation matrix with       */
/* real power coeffecients                                                                  */
/* The *_cmpx functions of the above are used to compute the  correlation matrix with       */
/* complex field coeffecients                                                               */
/*                                                                                          */
/* Created by : Xu Xiaojun, Zeng Meng                                                       */
/*              National Communications Lab,UESTC        	                                */
/*              16 April 2005                                                               */
/* Revised on : 25 Aptil,10 May, 2005                                                       */
/* Version 1.0                                                                              */
/* (c) Xu Xiaojun,ZengMeng 2005, NCL-UESTC                                                  */      
/* **************************************************************************************** */

#include "float.h"             /* function _isnan() in it is used to judge if a number is NaN */
#include "complex.h"           /* operation of complex numbers  */
#include "numerical_recipes.h" /* functions of numerical recipes */

#define sqrt3  1.7320508
#define sqrt2  1.4142136
#define PI  3.1415926536


/*********************/
/*    Uniform PAS    */
/*********************/

      /********************************/
      /* Normalization of uniform PAS */
      /********************************/
      
/** compute the power normalisation coef Q such that the PAS can be regarded as probability distribution **/
void normalisation_uniform(int cluster_number, float power_lin[], float AS_deg[],float *Q)
{
    int   i,j,k;
    float *delta_phi_rad,*b;
    float **A;

    delta_phi_rad=vector(1,cluster_number);
    b = vector(1,cluster_number);

    A = matrix(1,cluster_number,1,cluster_number);
    
/** initialize A=zeros(cluster_number) **/
    for(i=1;i<=cluster_number;i++)
	{
        for(j=1;j<=cluster_number;j++)
            A[i][j] = 0.0;
	}
       
/** initialize b=zeros(cluster_number,1) **/
    for(j=1;j<=cluster_number;j++)
	{
        b[j]=0.0;
	}
     
/** set b=[0,0,0.5] **/
    b[cluster_number] = 0.5;
   
/** convert AS from degree to rad **/
    for(i=1;i<=cluster_number;i++)
	{
        delta_phi_rad[i] = (AS_deg[i-1]*PI/180)*sqrt3;
	}


    if(cluster_number==1)
    {
		Q[1] = 1/(2*delta_phi_rad[1]);
    }
    else
    {
	    for(i=1;i<=cluster_number-1;i++)
		{
	        A[i][1] = 1/power_lin[0];
		}
	    for(k=2;k<=cluster_number;k++)
		{
	        A[k-1][k] = (-1.0)/power_lin[k-1];
		}
    	for(i=1;i<=cluster_number;i++)
		{
	        A[cluster_number][i] = delta_phi_rad[i];
		}
	    
/** compute AQ=b, and Q is saved in b **/
    inverse(A,cluster_number,b);
    for(i=1;i<=cluster_number;i++)
	{
        Q[i] = b[i];
	}
	}
	
   

}

      /******************************/
      /* Computation of Rxx_uniform */
      /******************************/
      
/** compute the correlation of the real/imaginary parts of signals in the case of uniform PAS
/*  at the spacing given by d_norm. **/
/** The PAS is charaterised by the AOA phi_0_deg and by AS  AS_deg
/*  result is a 1*3 vector, and by toeplitz transformation, it is turned into the desired Rxx **/
float Rxx_uniform(float d_norm,float phi_0_deg,float AS_deg )
{
    float D,phi_0_rad,AS_rad,result,delta_phi_uniform_rad,temp_xx_uniform;
    int m;
    
/** Convert parameters from degree to rad **/
    D = 2*PI*d_norm;
    phi_0_rad = phi_0_deg*PI/180;
    AS_rad = AS_deg*PI/180;
    delta_phi_uniform_rad = AS_rad*sqrt3;
    
/** initialisation **/
    result = 0.0;
    temp_xx_uniform = 1.0;

/** Computation of Rxx **/ 
    m = 0;
    while(m<100)
    {
        m++;
    	temp_xx_uniform = 2*bessj(2*m,D)*sin(2*m*delta_phi_uniform_rad)*cos(2*m*phi_0_rad)/m;
        result += temp_xx_uniform;
    }
     
    return result;
    }

      /******************************/
      /* Computation of Rxy_uniform */
      /******************************/
      
/** compute the correlation of the real and imaginary parts of signals in the case of uniform PAS
/*  at the spacing given by d_norm. **/
/** The PAS is charaterised by the AOA phi_0_deg and by AS  AS_deg
/*  result is a 1*3 vector, and by toeplitz transformation, it is turned into the desired Rxy **/
float Rxy_uniform(float d_norm,float phi_0_deg,float AS_deg)
{
    float D,phi_0_rad,AS_rad,result,delta_phi_uniform_rad,temp_xx_uniform;
    int m;

/** Convert parameters from degree to rad **/
    D = 2*PI*d_norm;
    phi_0_rad = phi_0_deg*PI/180;
    AS_rad = AS_deg*PI/180;
    delta_phi_uniform_rad = AS_rad*sqrt3;
    
/** Initialization **/
    result = 4*bessj1(D)*sin(delta_phi_uniform_rad)*sin(phi_0_rad);;
    temp_xx_uniform = 1.0;

/** Computation of Rxy **/ 
    m = 0;
    while(m<100)
    {
        m++;
        temp_xx_uniform = 4*bessj(2*m+1,D)*sin((2*m+1)*delta_phi_uniform_rad)*sin((2*m+1)*phi_0_rad)/(2*m+1);
        result += temp_xx_uniform;
    }
    
    return result;
    }


      /******************************/
      /* correlation of real matrix */
      /******************************/
      
float **uniform_real(int N,float spacing,int cluster_number,float amplitude_cluster[],float phi_deg[],
                     float AS_deg[],float delta_phi_deg[])
{
    int i,j,k;
    float *d_norm;
    d_norm = vector(1,N);
    
/** Computation of d_norm in case of ULA **/
    for(i=1;i<=N;i++)
        d_norm[i] = (i-1)*spacing;
        
    float *Rxx,*Rxy;
    float *temp;
    float result_x,result_y;
    float *Q;
    float **R;
    
    Rxx = vector(1,N);
    Rxy = vector(1,N);
    temp = vector(1,N);
    Q = vector(1,cluster_number);
    R = matrix(1,N,1,N);
    
/** Computation of real correlation of uniform PAS **/
    if(N==1)
	{
        R[N][N] = 1;
	}
    else
    {
	    for(i=1;i<=N;i++)
		{
	        Rxx[i] = bessj0(2*PI*d_norm[i]);
	        Rxy[i] = 0.0;
		}
        normalisation_uniform(cluster_number,amplitude_cluster,AS_deg,Q);

    	for(i=1;i<=N;i++)
		{
	        for(k=1;k<=cluster_number;k++)
			{
				result_x = Rxx_uniform(d_norm[i],phi_deg[k-1],AS_deg[k-1]);
		        result_y = Rxy_uniform(d_norm[i],phi_deg[k-1],AS_deg[k-1]);
	     	    Rxx[i]  += Q[k]*result_x;
		        Rxy[i]  += Q[k]*result_y;
			}
		}
 			
    	for(i=1;i<=N;i++)
		{
	        temp[i] = Rxx[i]*Rxx[i]+Rxy[i]*Rxy[i];
		}
	    
/** Toeplitz operation **/
        for(i=1;i<=N;i++)
            for(j=i,k=1;j<=N;j++,k++)
			{
	            R[i][j] = temp[k];
	            R[j][i] = R[i][j];
			}
	}
    
    return R;
    
}

      /*********************************/
      /* Correlation of complex matrix */
      /*********************************/
      
complex **uniform_cmpx(int N,float spacing,int cluster_number,float amplitude_cluster[],float phi_deg[],
                       float AS_deg[],float delta_phi_deg[])
{
	int i,j,k;
	float *d_norm;
	d_norm = vector(1,N);
	for(i=1;i<=N;i++)
	{
	    d_norm[i] = (i-1)*spacing;
	}
	float *Rxx,*Rxy;
	complex *temp;
	float result_x,result_y,*Q;
	complex **R;
	Rxx = vector(1,N);
	Rxy = vector(1,N);
	temp = cvector(1,N);
    Q = vector(1,cluster_number);

	R = cmatrix(1,N*N,1,N*N);

/** Compute compelx correlation of uniform PAS **/	
    if(N==1)
	{
        R[N][N] = complex(1.0,0.0);	
	}

    else
    {