correlation.h

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

		        Rxy[i]  += Q[k]*result_y;
			}
		}
			
        for(i=1;i<=N;i++)
		{
			temp[i] = complex(Rxx[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] = conj(R[i][j]); 
			}

    }
	return R;
}




/**********************/
/*   Laplacian PAS    */
/**********************/


      /**********************************/
      /* Normalisation of laplacian PAS */
      /**********************************/

/** Compute the power normalisation coef Q such that the PAS can be regarded as probability distribution 
/*  normalisation of laplacian PAS, to derive Q **/
void normalisation_laplacian(int cluster_number, float power_lin[], float AS_deg[],float delta_phi_deg[],
							 float *Q,float *sigma_deg)
{
    int i,j,k;
    float *delta_phi_rad,*b,*AS_rad,*sigma_rad;
    float **A;

    AS_rad = vector(1,cluster_number);
    delta_phi_rad = vector(1,cluster_number);
    sigma_rad = vector(1,cluster_number);
    b = vector(1,cluster_number);
    A = matrix(1,cluster_number,1,cluster_number);

/** Computation of sigma_deg  **/
    for(i=1;i<=cluster_number;i++)
	{
        sigma_deg[i] = AS_deg[i-1];
	}


/** 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,1] **/
    b[cluster_number] = 1;
   
/** Convert AS from degree to rad **/
    for(i=1;i<=cluster_number;i++)
    {
	    AS_rad[i] = AS_deg[i-1]*PI/180;
	    delta_phi_rad[i] = (delta_phi_deg[i-1]*PI/180)*sqrt2;
	    sigma_rad[i] = sigma_deg[i]*PI/180;
    }

/** Computation of Q **/
    if(cluster_number==1)
    {
        Q[1] = 1/(1-exp((-1)*sqrt2*delta_phi_rad[1]/sigma_rad[1]));
    }
    else
    {
	    for(i=1;i<=cluster_number-1;i++)
		{
	        A[i][1] = 1/(sigma_rad[1]*power_lin[0]);
		}
	    for(k=2;k<=cluster_number;k++)
		{
	        A[k-1][k] = (-1)/(sigma_rad[k]*power_lin[k-1]);
		}
	    for(i=1;i<=cluster_number;i++)
		{
	        A[cluster_number][i] = 1-exp((-1)*sqrt2*delta_phi_rad[i]/sigma_rad[i]);
		}
		
/** Compute AQ=b, and Q is save in b **/
    inverse(A,cluster_number,b);
    for(i=1;i<=cluster_number;i++)
	{
	    Q[i] = b[i];
	}
	}

}

/**********************************************/
/*----------------Rxx_laplacian---------------*/
/**********************************************/

/** Compute the correlation of the real/imaginary parts of signals in the case of laplacian 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_laplacian(float d_norm,float phi_0_deg,float sigma_deg,float delta_phi_deg)
{
    float D,phi_0_rad,sigma_rad,result,delta_phi_rad,temp_xx_laplacian;
    int m;
  
/** Convert parameters from degree to rad **/
    D = 2*PI*d_norm;
	phi_0_rad = phi_0_deg*PI/180;
	sigma_rad = sigma_deg*PI/180;
	delta_phi_rad = delta_phi_deg*PI/180;
	
/** Initialisation **/
	result = 0.0;
    temp_xx_laplacian = 1.0;
   	m = 0;

/** Computation of Rxx **/ 

    while(m<10)
	{
        m++;
	    temp_xx_laplacian = 4*bessj(2*m,D)*cos(2*m*phi_0_rad)*(sqrt2/sigma_rad
			                +exp((-1)*sqrt2*delta_phi_rad/sigma_rad)*(2*m*sin(2*m*delta_phi_rad)-sqrt2
						    *cos(2*m*delta_phi_rad)/sigma_rad))/(4*sqrt2*sigma_rad*m*m+2*sqrt2/sigma_rad);
	    result           += temp_xx_laplacian;

	
	}
    return result;
}

/**********************************************/
/*----------------Rxy_laplacian---------------*/
/**********************************************/

/** Compute the correlation of the real and imaginary parts of signals in the case of laplacian 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_laplacian(float d_norm,float phi_0_deg,float sigma_deg,int delta_phi_deg)
{ 
    float D,phi_0_rad,sigma_rad,result,delta_phi_rad,temp_xx_laplacian;
    int m;
  
/** Convert parameters from degree to rad **/
    D = 2*PI*d_norm;
	phi_0_rad = phi_0_deg*PI/180;
	sigma_rad = sigma_deg*PI/180;
	delta_phi_rad = delta_phi_deg*PI/180;

/** Initialization **/
	m = 0;
	result = (4*bessj1(D)*sin(phi_0_rad)*((sqrt2/sigma_rad)-(exp((-1)*sqrt2*delta_phi_rad/sigma_rad)
		     *(sin(delta_phi_rad)+sqrt2*cos(delta_phi_rad)/sigma_rad))))
			 /(sqrt2*sigma_rad*((2/(sigma_rad*sigma_rad)+1)));
    temp_xx_laplacian = 1.0;
	
/** Computation of Rxy **/ 
    while(m<10)
    {
     	m++;
    	temp_xx_laplacian = (4*bessj(2*m+1,D)*sin((2*m+1)*phi_0_rad)*((sqrt2/sigma_rad)-(exp((-1)
			                 *sqrt2*delta_phi_rad/sigma_rad)*((2*m+1)*sin((2*m+1)*delta_phi_rad)
							 +sqrt2*cos((2*m+1)*delta_phi_rad)/sigma_rad))))
							 /(sqrt2*sigma_rad*((2/(sigma_rad*sigma_rad)+(2*m+1)*(2*m+1))));
	    result           += temp_xx_laplacian;

    }
    return result;
}


/**************************************************/
/*-----------correlation of real matrix-----------*/
/**************************************************/

float **laplacian_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);
	for(i=1;i<=N;i++)
	{
		d_norm[i] = (i-1)*spacing;
	}

	float *Rxx,*Rxy;
	float *temp;
	float result_x,result_y,*Q,*sigma_deg;
	float **R;
	
	Rxx = vector(1,N);
	Rxy = vector(1,N);
	temp = vector(1,N);
    Q = vector(1,cluster_number);
    sigma_deg = vector(1,cluster_number);
	R = matrix(1,N,1,N);
	 
/** Compute real correlation of laplacian 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_laplacian(cluster_number,amplitude_cluster,AS_deg,delta_phi_deg,Q,sigma_deg);
	for(i=1;i<=N;i++)
	{
		for(k=1;k<=cluster_number;k++)
		{
			result_x = Rxx_laplacian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
			result_y = Rxy_laplacian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_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 **laplacian_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,*sigma_deg;
	complex **R;
	Rxx = vector(1,N);
	Rxy = vector(1,N);
	temp = cvector(1,N);
    Q = vector(1,cluster_number);
    sigma_deg = vector(1,cluster_number);
	R = cmatrix(1,N*N,1,N*N);

/** Compute compelx correlation of laplacian PAS **/	
    if(N==1)
	{
        R[N][N] = complex(1.0,0.0);	
	}
    else
    {
	    for(i=1;i<=N;i++)
		{
		    Rxx[i] = bessj0(2*PI*d_norm[i]);
		    Rxy[i] = 0.0;
		}
        normalisation_laplacian(cluster_number,amplitude_cluster,AS_deg,delta_phi_deg,Q,sigma_deg);
	    for(i=1;i<=N;i++)
		{
		    for(k=1;k<=cluster_number;k++)
			{
	    	    result_x = Rxx_laplacian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
	    	    result_y = Rxy_laplacian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
		        Rxx[i] += Q[k]*result_x;
		        Rxy[i] += Q[k]*result_y;
			}
		}
			
	    for(i=1;i<=N;i++)
		{
			temp[i] = complex(Rxx[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] =conj(R[i][j]); 
			}
 
	}
    return R;
}