biintpmimo_channel.c

2发2收的MIMO信道模拟程序

/*******************************************************************************************************
2发2收的MIMO信道模拟程序，信道满足瑞利衰落，采用最新的改进的Jakes模型，解决了传统算法中的统计性缺陷。
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define	PI		4*atan(1)

typedef struct
{
	double real;
	double image;
}COMPLEX;

//Add two complex a and b.
COMPLEX add_cmx(COMPLEX a, COMPLEX b)
{
	COMPLEX result;

	result.real  = a.real  + b.real;
	result.image = a.image + b.image;

	return result;
}

//Multiple complexs a and b.
COMPLEX mul_cmx(COMPLEX a, COMPLEX b)
{
	COMPLEX result;

	result.real  = a.real*b.real  - a.image*b.image;
	result.image = a.real*b.image + a.image*b.real;

	return result;
}
// generate channel data for current time
int GenCurCHData(double * alpha, double * beta, double * p_Time, COMPLEX * p_CHData, int PathNum, double* p_Phase, double* p_CorrMatrix)  //PathNum=TotalPathNum
{
	int n;
	int n_Path;
	int N0 = 25;
	//int N  = 4*N0;
	
    double SQRT_N0 = sqrt(N0);
  
	COMPLEX	 temp;
	COMPLEX  CHtemp;
	double  *p_curCorrMatrix;
	
	COMPLEX *p_IID_CHData;	// i.i.d channel data
	p_IID_CHData = (COMPLEX *)malloc(PathNum*sizeof(COMPLEX));
	
	//Generate i.i.d channel data for current TimeInstant
	for(n_Path=0; n_Path<PathNum; n_Path++)
	{
		CHtemp.real = CHtemp.image = 0;
		for(n=0; n<N0; n++)
		{
			temp.real   = cos( (*(alpha + n*PathNum + n_Path))*(*p_Time) + p_Phase[0]);   // alpha=0-191 192-383 ... 24*192-(25*192-1) 
			temp.image  = sin( (*(beta + n*PathNum + n_Path))*(*p_Time) + p_Phase[1]);    // beta=0-191 192-383 ... 24*192-(25*192-1) 
			
			CHtemp.real  = CHtemp.real+temp.real;
			CHtemp.image = CHtemp.image+temp.image;
		}
		
		p_IID_CHData[n_Path].real  = (CHtemp.real)/SQRT_N0;
		p_IID_CHData[n_Path].image = (CHtemp.image)/SQRT_N0;
	}
	
	//Add correlation for i.i.d channel data
	p_curCorrMatrix = p_CorrMatrix;
	for(n_Path=0; n_Path<PathNum; n_Path++) 
	{
		p_CHData[n_Path].real = 0;
		p_CHData[n_Path].image = 0;		
		for(n=0; n<PathNum; n++)
		{
			//p_CHData[n_Path].real  += p_CorrMatrix[n_Path*PathNum+n] * p_IID_CHData[n].real;
			//p_CHData[n_Path].image += p_CorrMatrix[n_Path*PathNum+n] * p_IID_CHData[n].image;
			p_CHData[n_Path].real  += p_curCorrMatrix[n] * p_IID_CHData[n].real;  // p_curCorrMatrix 是一维
			p_CHData[n_Path].image += p_curCorrMatrix[n] * p_IID_CHData[n].image;
		}
		p_curCorrMatrix += PathNum;         ///
	}
	free(p_IID_CHData);
	
	return n_Path;	
}

// intepolate for real, y = a*x.^2 + b*x + c
void dblBiInterp(double* p_IntpBegin, double* p_IntpMiddle, double* p_IntpEnd, int VecLen, double* p_Out, int n_OutLen, int BiInterpCoef)
{
	int i;
	int j;
	double a;
	double b;
	double c;
	
	for(i=0; i<VecLen; i++)
	{
		a =  2*p_IntpBegin[i] - 4*p_IntpMiddle[i] + 2*p_IntpEnd[i];
		b = -3*p_IntpBegin[i] + 4*p_IntpMiddle[i] -   p_IntpEnd[i];
		c =    p_IntpBegin[i];
		
		//for(j=0; j<BiInterpCoef; j++)
		for(j=0; j<n_OutLen; j++)
		{			
			//pp_Out[j][i] = a*j*j/BiInterpCoef/BiInterpCoef + b*j/BiInterpCoef + c;
			p_Out[j*VecLen+i] = a*j*j/BiInterpCoef/BiInterpCoef + b*j/BiInterpCoef + c;//  0-191 ... 192*255-192*256-1
		}
	}	
}

// intepolate for complex
void cmxBiInterp(COMPLEX* p_IntpBegin, COMPLEX* p_IntpMiddle, COMPLEX* p_IntpEnd, int VecLen, double* p_Out_real, double* p_Out_image, int n_OutLen, int BiInterpCoef)
{
	int i;
	
	double *p_IntpBegin_real;
	double *p_IntpBegin_image;
	double *p_IntpMiddle_real;
	double *p_IntpMiddle_image;
	double *p_IntpEnd_real;
	double *p_IntpEnd_image;
	
	p_IntpBegin_real   =  (double*)malloc(VecLen*sizeof(double));
	p_IntpBegin_image  =  (double*)malloc(VecLen*sizeof(double));
	p_IntpMiddle_real  =  (double*)malloc(VecLen*sizeof(double));
	p_IntpMiddle_image =  (double*)malloc(VecLen*sizeof(double));
	p_IntpEnd_real     =  (double*)malloc(VecLen*sizeof(double));
	p_IntpEnd_image    =  (double*)malloc(VecLen*sizeof(double));
	
	// get real part & image part of input vector
	for(i=0; i<VecLen; i++)
	{
		p_IntpBegin_real[i]		= p_IntpBegin[i].real;
		p_IntpBegin_image[i]	= p_IntpBegin[i].image;
		p_IntpMiddle_real[i]	= p_IntpMiddle[i].real;
		p_IntpMiddle_image[i]	= p_IntpMiddle[i].image;
		p_IntpEnd_real[i]		= p_IntpEnd[i].real;
		p_IntpEnd_image[i]		= p_IntpEnd[i].image;
	}
	
	// interpolation of real part & image part
	dblBiInterp(p_IntpBegin_real, p_IntpMiddle_real, p_IntpEnd_real, VecLen, p_Out_real, n_OutLen, BiInterpCoef);//?? why do so VecLen=TotalPathnum
	dblBiInterp(p_IntpBegin_image, p_IntpMiddle_image, p_IntpEnd_image, VecLen, p_Out_image, n_OutLen, BiInterpCoef);
	
	// free!
	free(p_IntpBegin_real);
	free(p_IntpBegin_image);
	free(p_IntpMiddle_real);
	free(p_IntpMiddle_image);
	free(p_IntpEnd_real);
	free(p_IntpEnd_image);

}

// pass through multiple path Rayleigh fading channel
void MIMO_Channel(double *p_txSignal_r, double *p_txSignal_i, int SigLenPerTx, int TxAntennaNum,
				  int RxAntennaNum, double *p_PathDelay, int maxDelay, int PathNum,
				  double *p_CorrMatrix, double fc, double speed, double Tc, 
				  double *p_Phase, double TimeBegin, double *p_rxSignal_r, double *p_rxSignal_i
				  )
{
	int i = 0;	
	int n_Sample= 0;
	int n_Rx	= 0;
	int n_Path	= 0;
	int n_Tx	= 0;
	int Index	= 0;		// index of transmit signal
	int n_rxSignal= 0;		// index of receive signal	
	int n_Segment = 0;		// index for segment
	
	int n = 0;
	int N0 = 25;
	int N = 4*N0;
	double Fd = speed*fc/(3.0e8)/3.6;	
	double Wd = 2*PI*Fd;
	double tt;
	double *alpha,*beta;
	double gamma;

	int tmpInterpCoef	= (int)(1.0/8.0/Fd/Tc);					// calculational interpolate factor ? 1.645714285714286e+004
	int InterpCoef		= (tmpInterpCoef > 128 ? 128:tmpInterpCoef);	// real interpolate factor
	int BiInterpCoef	= 2*InterpCoef;							// for bi_interpolation	
	int TotalPathNum	= RxAntennaNum*TxAntennaNum*PathNum;	// total path number
	int SegChDataLen	= TotalPathNum*BiInterpCoef;			// channel data number of one segment	 
	int txSigLen		= SigLenPerTx*TxAntennaNum;				// data number of all Txs
	int SampleLen		= SigLenPerTx + maxDelay;				// sample number
	int rxSigLen		= RxAntennaNum*SampleLen;				// data number of all Rxs
	int SegNum			= (SampleLen -1) / BiInterpCoef;					// segment number of interpolation
	int lastSegLen		= (SampleLen -1) % BiInterpCoef;				// last segment length for special process
	int lastSegChDataLen= TotalPathNum*(lastSegLen+1);

	COMPLEX *p_CHData;		// pointer to channel datas of one segment
	double  *p_CHData_r;	// real part of p_CHData
	double  *p_CHData_i;	// image part of p_CHData
	COMPLEX *p_curCHData;	// pointer to channel datas of current sample
	COMPLEX *p_txSignal;	// complex transmit signals
	COMPLEX  tmp;

	COMPLEX *p_IntpBegin;	// 1st sample channel datas of all paths for interpolation
	COMPLEX *p_IntpMiddle;	// 2nd sample channel datas of all paths for interpolation
	COMPLEX *p_IntpEnd;		// 3rd sample channel datas of all paths for interpolation
	COMPLEX *p_Tmp;	

	//for interpolate	
	p_IntpBegin = (COMPLEX*)malloc(TotalPathNum*sizeof(COMPLEX));
	p_IntpMiddle= (COMPLEX*)malloc(TotalPathNum*sizeof(COMPLEX));
	p_IntpEnd   = (COMPLEX*)malloc(TotalPathNum*sizeof(COMPLEX));
	
	p_CHData	= (COMPLEX*)malloc(SegChDataLen*sizeof(COMPLEX));
	p_CHData_r	= (double* )malloc(SegChDataLen*sizeof(double ));	
	p_CHData_i	= (double* )malloc(SegChDataLen*sizeof(double ));													
	p_txSignal	= (COMPLEX*)malloc(txSigLen*sizeof(COMPLEX));
	alpha		= (double* )malloc(N0*TotalPathNum*sizeof(double));
	beta		= (double* )malloc(N0*TotalPathNum*sizeof(double));

	// combine real & image part of transmit signal
	for(i=0; i<txSigLen; i++)
	{
		p_txSignal[i].real	= p_txSignal_r[i];
		p_txSignal[i].image = p_txSignal_i[i];
	}

	// generate alpha,beta for all samples
	for(n_Path=0; n_Path<TotalPathNum; n_Path++)	
	{
		for(n=0; n<N0; n++)
		{
			gamma = PI*( (double)2*n + (double)2*n_Path/TotalPathNum + (double)1/(2*TotalPathNum) )/N;// 0-191 192-383 ... 24*192-(25*192-1)
			alpha[n*TotalPathNum + n_Path] = Wd*cos(gamma);         ///M*N*L*N0=4800?
			beta[n*TotalPathNum + n_Path]  = Wd*sin(gamma);
		}
	}
	
	// tt = (1 + TimeBegin + 1e8)*Tc;
	tt = TimeBegin;
	
	//GenCurCHData(alpha,  beta, &tt, p_IntpBegin, TotalPathNum, p_Phase, p_CorrMatrix);
	GenCurCHData(alpha,  beta, &tt, p_IntpBegin, TotalPathNum, p_Phase, p_CorrMatrix);

	// front segments
	for(n_Segment=0; n_Segment<SegNum; n_Segment++)
	{
		tt += InterpCoef * Tc;		
		GenCurCHData(alpha,  beta, &tt, p_IntpMiddle, TotalPathNum, p_Phase, p_CorrMatrix);
		
		tt += InterpCoef * Tc;		
		GenCurCHData(alpha,  beta, &tt, p_IntpEnd, TotalPathNum, p_Phase, p_CorrMatrix);
			
		// compute channel datas for one segment
		cmxBiInterp(p_IntpBegin, p_IntpMiddle, p_IntpEnd, TotalPathNum, p_CHData_r, p_CHData_i, BiInterpCoef, BiInterpCoef);
		
		// combine real & image part of channel datas
		for(i=0; i<SegChDataLen; i++)
		{
			p_CHData[i].real  = p_CHData_r[i];    //p_CHData(192*256)
			p_CHData[i].image = p_CHData_i[i];
		}
		p_curCHData = p_CHData;

		// calculate receive signals of one segment
		for(n=0; n<BiInterpCoef; n++)
		{
			//n_Sample = n_Segment*BiInterpCoef + n;
			//p_curCHData = p_CHData+n*TotalPathNum;
			
			for(n_Rx=0; n_Rx<RxAntennaNum; n_Rx++)
			{
				tmp.real = 0;
				tmp.image = 0;
				for(n_Path=0; n_Path<PathNum; n_Path++)   // PathNum=6 
				{
					Index = n_Sample-(int)(p_PathDelay[n_Path]);
					
					if( Index>=0 &&	Index<SigLenPerTx)
					{
						for(n_Tx=0; n_Tx<TxAntennaNum; n_Tx++)
						{
							tmp = add_cmx( tmp, mul_cmx( p_txSignal[Index*TxAntennaNum+n_Tx],
								p_curCHData[(n_Path*TxAntennaNum+n_Tx)*RxAntennaNum+n_Rx] ) );    // p_curCHData : [(0 8 16 24) ... (160 168 176 184)] [(1 9 17 25)...]  .... tmp=tmp+s(n_Tx,Index)*cur(n_Rx+(n_Path*M+n_Tx)*RxAntennaNum) ?
						}    	
					}
				}
				p_rxSignal_r[n_rxSignal] = tmp.real;
				p_rxSignal_i[n_rxSignal] = tmp.image;
				
				// index of next receive signal
				n_rxSignal++;
			}
			// index of sample
			n_Sample++;
			// pointer to current sample
			p_curCHData += TotalPathNum;   //     p_curCHData(192*256)
		}	
		
		// replace pointer of p_IntpBegin & p_IntpEnd
		p_Tmp		= p_IntpBegin;
		p_IntpBegin = p_IntpEnd;
		p_IntpEnd	= p_Tmp;
	}

	// last senment
	if(lastSegLen)
	{
		tt += InterpCoef * Tc;
		GenCurCHData(alpha,  beta, &tt, p_IntpMiddle, TotalPathNum, p_Phase, p_CorrMatrix);
		
		tt += InterpCoef * Tc;
		GenCurCHData(alpha,  beta, &tt, p_IntpEnd, TotalPathNum, p_Phase, p_CorrMatrix);
			
		// compute channel datas for the last segment
		cmxBiInterp(p_IntpBegin, p_IntpMiddle, p_IntpEnd, TotalPathNum, p_CHData_r, p_CHData_i, lastSegLen+1, BiInterpCoef);
		
		// combine real & image part of channel datas
		for(i=0; i<lastSegChDataLen; i++)
		{
			p_CHData[i].real  = p_CHData_r[i];
			p_CHData[i].image = p_CHData_i[i];
		}
		p_curCHData = p_CHData;
		
		// calculate receive signals of one segment
		for(n=0; n<lastSegLen+1; n++)
		{
			//n_Sample = n_Segment*BiInterpCoef + n;
			//p_curCHData = n*TotalPathNum;

			for(n_Rx=0; n_Rx<RxAntennaNum; n_Rx++)
			{
				tmp.real = 0;
				tmp.image = 0;
				for(n_Path=0; n_Path<PathNum; n_Path++)     // PathNum=6 
				{
					Index = n_Sample-(int)(p_PathDelay[n_Path]);
					
					if( Index>=0 &&	Index<SigLenPerTx)
					{
						for(n_Tx=0; n_Tx<TxAntennaNum; n_Tx++)
						{
							tmp = add_cmx( tmp, mul_cmx( p_txSignal[Index*TxAntennaNum+n_Tx],
								p_curCHData[(n_Path*TxAntennaNum+n_Tx)*RxAntennaNum+n_Rx] ) );
						}    	
					}
				}
				p_rxSignal_r[n_rxSignal] = tmp.real;
				p_rxSignal_i[n_rxSignal] = tmp.image;

				// index of next receive signal
				n_rxSignal++;
			}
			// index of sample
			n_Sample++;		
			// pointer to current sample
			if(n<lastSegLen)
			{
				p_curCHData += TotalPathNum;
			}
		}
	}

	// for the last sample
	else
	{	
		// n_Sample = rxSigLen - 1;
		//p_curCHData = p_IntpBegin;	
		for(n_Rx=0; n_Rx<RxAntennaNum; n_Rx++)
		{
			tmp.real = tmp.image = 0;
			for(n_Path=0; n_Path<PathNum; n_Path++)
			{
				Index = n_Sample-(int)(p_PathDelay[n_Path]);
				
				if( Index>=0 &&	Index<SigLenPerTx)
				{
					for(n_Tx=0; n_Tx<TxAntennaNum; n_Tx++)
					{
						tmp = add_cmx( tmp, mul_cmx( p_txSignal[Index*TxAntennaNum+n_Tx],
							p_IntpBegin[(n_Path*TxAntennaNum+n_Tx)*RxAntennaNum+n_Rx] ) );
					}    	
				}
			}
			p_rxSignal_r[n_rxSignal] = tmp.real;
			p_rxSignal_i[n_rxSignal] = tmp.image;

			// index of next receive signal
			n_rxSignal++;
		}
	}

	free(p_CHData);
	free(p_CHData_r);
	free(p_CHData_i);
	free(p_txSignal);

	//free( (char *)alpha );
	//free( (char *)beta );
	free(alpha );
	free(beta );
	free(p_IntpBegin);
	free(p_IntpMiddle);
	free(p_IntpEnd);
}



void main()
{
	int i,j;
	int SigLenPerTx  = 1000;
	int TxAntennaNum = 2;
	int RxAntennaNum = 2;
	double PathDelay[2] = {0, 5};
	int maxDelay = 5;
	int SampleLen = SigLenPerTx+maxDelay;
	int PathNum  = 2;
	double fc = 3.5e9;
	double speed = 500;
	double Tc = 1/(1.28e6);
	double Phase[2] = {1.8, 0.8};
	double TimeBegin = Tc*1e8 + 1e5;
	FILE *fp1 = fopen("rx_r.dat", "w");
	FILE *fp2 = fopen("rx_i.dat", "w");
	
	double *p_CorrMatrix = (double*)malloc(64*sizeof(double));
	double *p_txSignal_r = (double*)malloc(TxAntennaNum*SigLenPerTx*sizeof(double));
	double *p_txSignal_i = (double*)malloc(TxAntennaNum*SigLenPerTx*sizeof(double));
	double *p_rxSignal_r = (double*)malloc(RxAntennaNum*SampleLen*sizeof(double));
	double *p_rxSignal_i = (double*)malloc(RxAntennaNum*SampleLen*sizeof(double));

	for(i=0; i<8; i++)
	{
		for(j=0; j<8; j++)
		{
			if(i==j)
				p_CorrMatrix[i*8+j] = 1.0;
			else 
				p_CorrMatrix[i*8+j] = 0.0;
		}
	}

	for(i=0; i<2000; i++)
	{
		p_txSignal_r[i] = 1.0;
		p_txSignal_i[i] = 1.0;
	}

	MIMO_Channel(p_txSignal_r, p_txSignal_i, SigLenPerTx, TxAntennaNum, RxAntennaNum, PathDelay, maxDelay,
					PathNum, p_CorrMatrix, fc, speed, Tc, Phase, TimeBegin, p_rxSignal_r, p_rxSignal_i);
	for(i=0; i<2010; i++)
	{
		fprintf(fp1, "%f\n", p_rxSignal_r[i]);
		fprintf(fp2, "%f\n", p_rxSignal_i[i]);
	}
	fclose(fp1);
	fclose(fp2);
	free(p_CorrMatrix);
	free(p_txSignal_r);
	free(p_txSignal_i);
	free(p_rxSignal_r);
	free(p_rxSignal_i); 
}