correlation.h

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

	    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] = 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;
}





/*********************/
/*   Gaussian PAS    */
/*********************/

      /*********************************/
      /* Normalization of Gaussian PAS */
      /*********************************/
      
/** compute the power normalisation coef Q such that the PAS can be regarded as probability distribution **/
void normalisation_gaussian(int cluster_number, float power_lin[], float AS_deg[],float delta_phi_deg[],
							float *Q, float *sigma_deg)
/** compute the power normalisation coef Q such that the PAS can be regarded as probability distribution **/
/** normalisation of gaussian PAS, to derive Q **/
{
    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.0/erff(delta_phi_rad[1]/(sqrt2*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] = erff(delta_phi_rad[i]/(sqrt2*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];
		}

	}

}

      /*******************************/
      /* Computation of Rxx_gaussian */
      /*******************************/
      
/** compute the correlation ofthe real/imaginary parts of signals in the case of gaussian 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_gaussian(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_gaussian;
    int m;
    complex A,operand1,operand2;
  
/** 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 **/
    m = 0;
    result = 0.0;
    temp_xx_gaussian = 1.0;
	
/** Computation of Rxx **/ 
    while(m<100)
    {
        m++;
	    operand1 = complex(delta_phi_rad/(sigma_rad*sqrt2),(-1)*m*sigma_rad*sqrt2);
    	operand2 = complex((-1)*delta_phi_rad/(sigma_rad*sqrt2),(-1)*m*sigma_rad*sqrt2);
	    A = erfcomp(operand1)-erfcomp(operand2);
	    if(_isnan(real(A))||_isnan(imag(A)))   /* break the loop if either the real or imaginary part is NaN */
		{
			break;
		}
    	temp_xx_gaussian = bessj(2*m,D)*cos(2*m*phi_0_rad)*exp((-2)*sigma_rad*sigma_rad*m*m)*real(A);
	    result += temp_xx_gaussian;
   	}
    return result;
}

      /*******************************/
      /* Computation of Rxy_gaussian */
      /*******************************/
/** Compute the correlation of the real and imaginary parts of signals in the case of gaussian 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_gaussian(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_xy_gaussian;
    int m;
    complex B,operand1,operand2;
  
/** 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 **/
    m = 0;
    operand1 = complex(delta_phi_rad/(sigma_rad*sqrt2),  (-1)*(m+0.5)*sigma_rad*sqrt2);
    operand2 = complex((-1)*delta_phi_rad/(sigma_rad*sqrt2),  (-1)*(m+0.5)*sigma_rad*sqrt2);
    B = erfcomp(operand1)-erfcomp(operand2);
       	   
    result = bessj1(D)*sin(phi_0_rad)*exp((-2)*sigma_rad*sigma_rad*0.25)*real(B);
	temp_xy_gaussian = 1.0;
	
/** Computation of Rxy **/ 
    while(m<100)
    {	
    	m++;
    	operand1=complex(delta_phi_rad/(sigma_rad*sqrt2),  (-1)*(m+0.5)*sigma_rad*sqrt2);
    	operand2=complex((-1)*delta_phi_rad/(sigma_rad*sqrt2),  (-1)*(m+0.5)*sigma_rad*sqrt2);
	    B=erfcomp(operand1)-erfcomp(operand2);
	    if(_isnan(real(B))||_isnan(imag(B)))  /* break the loop if either the real or imaginary part is NaN */
		{
			break; 
		}
	    temp_xy_gaussian = bessj(2*m+1,D)*sin((2*m+1)*phi_0_rad)
		                   *exp((-2)*sigma_rad*sigma_rad*(m+0.5)*(m+0.5))*real(B);
	    result += temp_xy_gaussian;
     
	}
    return result;
}


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

float **gaussian_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 gaussian 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_gaussian(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_gaussian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
		    	result_y = Rxy_gaussian(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 **gaussian_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 gaussian 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_gaussian(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_gaussian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
	            result_y = Rxy_gaussian(d_norm[i],phi_deg[k-1],sigma_deg[k],delta_phi_deg[k-1]);
		        Rxx[i]  += Q[k]*result_x;