hannelclass.cpp

用matlab程序实现WCDMA系统的仿真

	TempDopplerPSDPtr = DopplerPSDPtr;
	for (k=0; k<kd; k++) 
		// We are computing the square root of the DopplerPSD
		*TempDopplerPSDPtr++ = sqrt( 1.5 / (PI * FDopplerAdjusted * sqrt( 1.000 - ((deltaF * deltaF * (double) k * k)/(FDopplerAdjusted * FDopplerAdjusted)) ) ));

	//Use a polynomial fit to estimate the last point
	OrdinatePtr = TempDopplerPSDPtr - PolyFitLength;
	for (k=0; k<=PolyFitLength; k++) abcissa[k] = deltaF * (double) k;
	*TempDopplerPSDPtr = PolynomialFit( (double *) abcissa,OrdinatePtr,PolyFitLength,abcissa[4],&error);

	//Now create a complex vector of length N.
	//The first kd+1 elements will contain Gaussian random data that is shaped by the Doppler PSD.
	//The last kd elements will contain the conjugate of the 2nd element to the kd+1st element in reverse order

	FadingSpectrumPtr = (ComplexNumber *) calloc(N,sizeof(ComplexNumber));
	if (FadingSpectrumPtr == NULL)
	{
		printf("\nMemory allocation for FadingSpectrumPtr failed--exiting\n");
		return(NULL);
	}

	//Create the first kd+1 elements by generating Gaussian random data and scaling it by the doppler PSD
	TempFadingSpectrumPtr = FadingSpectrumPtr;
	TempDopplerPSDPtr = DopplerPSDPtr;
	for (k=0; k<=kd; k++)
		*TempFadingSpectrumPtr++ = ComplexRealMultiply( *TempDopplerPSDPtr++,GaussianRandomNumberGenerator(&InitFlag));

	//Create the last kd elements by by generating Gaussian random data and scaling it by the doppler PSD.
	//In this case the Doppler PSD applied to the Fading Spectrum array in a reverse order
	//i.e., FadingSpectrumArray[N-i] = DopplerPSD[i]*Noise
	TempFadingSpectrumPtr = FadingSpectrumPtr + N-1;
	TempDopplerPSDPtr = DopplerPSDPtr + 1;
	for (k=0; k<kd; k++) 
		*TempFadingSpectrumPtr-- = ComplexRealMultiply ( *TempDopplerPSDPtr++, GaussianRandomNumberGenerator(&InitFlag));

	//Compute the IFFT
	IFFToutputPtr = Radix2FFT(FadingSpectrumPtr,N,1);


	// IFFToutputPtr now contains the fading signal/
	// We must now scale it so that it has unit average energy
	Energy = 0.0;
	TempIFFToutputPtr = IFFToutputPtr;
	TempFadingSigPtr = FadingSigPtr;
	for (k=0; k<NSamples; k++)	
		Energy += ComplexEnergy( *TempIFFToutputPtr++ );

	Energy /= (double) NSamples;
	InvAmplitude = 1.00 / sqrt(Energy);

	TempIFFToutputPtr = IFFToutputPtr;
	for (k=0; k<NSamples; k++) *TempFadingSigPtr++ = ComplexRealMultiply(InvAmplitude, *TempIFFToutputPtr++);

	// Clean up arrays
	free(FadingSpectrumPtr);
	free(IFFToutputPtr);
	free(DopplerPSDPtr);

	//return signal
	return (FadingSigPtr);

}


ComplexNumber * ChannelClass::MultiPathChannel(ComplexNumber *SigLastptr,ComplexNumber *SigCurrentptr,
											   ComplexNumber *SigNextptr, unsigned long SignalLength)
/*****************************************************************************************************
/ComplexNumber * ChannelClass::MultiPathChannel(ComplexNumber *SigLastptr,ComplexNumber *SigCurrentptr,
/											   ComplexNumber *SigNextptr, unsigned long SignalLength)
/
/ Copyright 2002 The Mobile and Portable Radio Research Group
/
/ The function simulates the behavior of a specular multipath channel on the signal associated with 
/ the current frame (SigCurrentptr ).  The signal associated with the previous frame (SigLastptr) and 
/ the next frame (SigNextptr) are included so that their signal points can be incorporated as required 
/ when the multipath delays are incorporated.  These signals are placed into one large array called 
/ *CombinedSigPtr. We call this signal the "combined signal".  The combined signal contains three 
/ frames of signal points
/ 
/ For each multipath component, we must (1) determine the fading coefficients associated with that 
/ component, (2) interpolate those coefficients so that we have one coefficient per signal sample 
/ point, (3) extract the portion of the combined signal associated with the specific multipath 
/ component, (4) apply the fading to that signal, and (5) combine the result with the result from the 
/ other multipath components.
/
/ As an example, consider the case where we have two multipath components with at delay of 0 and 10 
/ respectively.  Then the output of this function will be a vector of length 
/ "FrameSampleLength + MaxDelay" (note that MaxDelay=10).  To construct the first multipath component, 
/ we must use that portion of the combined signal that corresponds to the entire current frame and 
/ 10 samples from the next frame.  We then must create the fading signal.  This is accomplished by 
/ extracting the appropriate fading coefficients and interpolating them so that we have one coefficient
/ per sample point.  A nearest neighbor interpolation technique is used.
/ 
/ To construct the second multipath comontent, which has a delay of 10, we must use 10 samples for the 
/ last portion of the previous frame and then the entire current frame.  This will account for the 
/ delay and the inter-symbol interference associated with the delay.  The fading signal is constructed
/ in a similar manner.  The fading signal is then multiplied to the multipath component.  This 
/ component is then combined with the first multipath component to create the channel output.
/
/ Paramters
/	Input
/		*SigLastptr			Pointer to the signal points associated with the previous frame
/		*SigCurrentptr		Pointer to the signal points associated with the current frame
/		*SigNextptr			Pointer to the signal points associated with the next frame
/		SignalLength		Length of signal associated with previous frame
/
/	Returns Pointer to the ComplexNumber array that contains the multipath channel response
/
******************************************************************************************************/

{
	unsigned AmplitudeLength = CHIPS_PER_FRAME/CHIPS_PER_FADE_SAMPLES;

	ComplexNumber *MpathSignalptr;					//Pointer to Multipath Signal array -- function returns this value
	unsigned NMpathSignal;							//length of Multipath signal array
	unsigned long i,k,j,LoopStop;					//Loop counters and parameters
	ComplexNumber *CombinedSigPtr;					//Pointers to the array that stores the combined signal
	ComplexNumber *TempCombinedSigPtr;				//Temporary Pointers to combined signal array
	unsigned long NCombinedSig;						//Length of Combined signal array

	//Temporary pointers for signal frames
	ComplexNumber *TempSigLastptr;					//Temporary pointers for SigLastptr 
	ComplexNumber *TempSigCurrentptr;				//Temporary pointers for SigCurrentptr
	ComplexNumber *TempSigNextptr;					//Temporary pointers for SigNextptr

	//Temporary pointers for multipath channel parameters
	unsigned *TempDelaysPtr;							//Temporary pointers for DelaysPtr
	unsigned MaxDelay;									//Stores the largest values in DelaysPtr
	
	//Temporary pointers for output array
	ComplexNumber *TempMpathSignalptr;//Temporary pointers for RealMpathSignalptr and ImagMpathSignalptr

	unsigned  FrameBegin;		//Determines the location of the beginning of the first frame in the Combined Signal
	unsigned  FrameEnd;			//Determines the location of the end of the first frame in the Combined Signal
	unsigned  FrameSampleLength;//Number of signal samples per frame

	double increment;				//Number of abcissa values between fading signal data points.  Note that the fading 
									//Signal data are stored in AmplitudesPtr 
	unsigned Increment;				//Number of abcissa values between fading signal data points.  Note that the fading 
									//Signal data are stored in AmplitudesPtr 

	unsigned PreSigFadeRatio;		//The next two variables are used to determine the number of fading coefficients 
	unsigned PreRemainingSamples;	//associated with the previous frame the number of samples from the previous frame
									//that must be included in the interpolation
	unsigned PostSigFadeRatio;		//The next two variables are used to determine the number of fading coefficients 
	unsigned PostRemainingSamples;	//associated with the next frame the number of samples from the next frame
									//that must be included in the interpolation
	unsigned faa;
	ComplexNumber *TempFadeptr; //Temporary pointers to the interpolatedfading coefficients
	ComplexNumber *AmpBeginPtr;	//Temporary pointer that marks the begining of the sequence of values that
								//must be included in the interpolation
	unsigned FadeSamplesPerFrame;	//The number of fade samples allocated for each frame

	MaxDelay=0;
	TempDelaysPtr=Delays;
	FrameSampleLength = SamplesPerChip * CHIPS_PER_FRAME;
	FrameBegin = (PulseLength + 1) >> 1;
	FrameEnd = FrameBegin - 1 + FrameSampleLength ;
	FadeSamplesPerFrame = CHIPS_PER_FRAME / CHIPS_PER_FADE_SAMPLES;


	/******************************************************************
	/* Create a combined Signal which contains the previous frame the 
	/* current frame and the future frame
	*/
	NCombinedSig = SignalLength + SignalLength + SignalLength - 2*(PulseLength - SamplesPerChip);
	CombinedSigPtr = (ComplexNumber*) calloc(NCombinedSig,sizeof(ComplexNumber));
	if (CombinedSigPtr  == NULL)
	{
		printf("\nMemory Allocation failed for CombinedSigPtr\n");
		return(NULL);
	}

	TempCombinedSigPtr = CombinedSigPtr;
	TempSigLastptr = SigLastptr;
	LoopStop = SignalLength - PulseLength + SamplesPerChip;
	for (i=0; i<LoopStop; i++) *TempCombinedSigPtr++ = *TempSigLastptr++;

	TempSigCurrentptr = SigCurrentptr;
	for (i=0; i<SignalLength; i++) *TempCombinedSigPtr++ = *TempSigCurrentptr++;

	TempSigNextptr = SigNextptr + PulseLength - SamplesPerChip;
	for (i=0; i<LoopStop; i++) *TempCombinedSigPtr++ = *TempSigNextptr++;

	/* Combined Signal Created
	/******************************************************************/
	//Find the maximum delay

	for (i=0;i<MultiPathComponents;i++)
	{
		if (*TempDelaysPtr > MaxDelay) MaxDelay = *TempDelaysPtr++;
		else TempDelaysPtr++;
	}
	//Maximum delay found
	/******************************************************************/
	//Now allocate space for the multipath signal
	NMpathSignal = SignalLength + MaxDelay;
	MultipathSignalLength = NMpathSignal;
	MpathSignalptr = (ComplexNumber *) calloc(NMpathSignal,sizeof(ComplexNumber));
	if (MpathSignalptr  == NULL)
	{
		printf("\nMemory Allocation failed for MpathSignalptr\n");
		return(NULL);
	}

	increment = (double) FrameSampleLength / (double) AmplitudeLength; //Abcissa increment between data points
	Increment = (unsigned long) increment;
//	RemainingSamples = NMpathSignal % Increment;
															
	TempDelaysPtr = Delays;
	for (i=0;i<MultiPathComponents;i++)
	{
		//Determine the number of fading coefficients need from the previous frame
//		PriorInterpPoints = (unsigned long) ceil((double) *TempDelaysPtr / increment);
//		PriorInterpPoints = (unsigned long) floor((double) *TempDelaysPtr / increment);

		PreSigFadeRatio = (FrameBegin - 1 + *TempDelaysPtr) / Increment;
		PostSigFadeRatio = (NMpathSignal - FrameEnd - *TempDelaysPtr) / Increment;

		PreRemainingSamples = (FrameBegin - 1 + *TempDelaysPtr) - (PreSigFadeRatio * Increment);
		PostRemainingSamples = (NMpathSignal - FrameEnd - *TempDelaysPtr) - (PostSigFadeRatio * Increment);

		AmpBeginPtr = *(FadingSignals + i) + FadeSamplesPerFrame - PreSigFadeRatio-1;

		//Determine the number of fading coefficients need from the next frame

		TempFadeptr = AmpBeginPtr ;
		TempCombinedSigPtr = CombinedSigPtr + FrameSampleLength - *TempDelaysPtr;
		TempMpathSignalptr = MpathSignalptr;

		for (j = 0; j < PreRemainingSamples; j++) 
			*TempMpathSignalptr++ = ComplexAdd(*TempMpathSignalptr, ComplexMultiply(*TempCombinedSigPtr++,*TempFadeptr));
		TempFadeptr++;

		faa = FadeSamplesPerFrame + PreSigFadeRatio + PostSigFadeRatio;
		for (j=0; j < faa; j++)
		{
			for (k=0; k < Increment; k++) 
				*TempMpathSignalptr++ = ComplexAdd(*TempMpathSignalptr, ComplexMultiply(*TempCombinedSigPtr++,*TempFadeptr));
			TempFadeptr++;
		}

		for (j = 0; j < PostRemainingSamples; j++)
		{
			*TempMpathSignalptr++ = ComplexAdd(*TempMpathSignalptr, ComplexMultiply(*TempCombinedSigPtr++,*TempFadeptr));
		}

		TempDelaysPtr++;

		/*******************************************************
		//Debug Code

		FILE *fp;

		fp = fopen("MpathSig.txt","w");
		TempMpathSignalptr = MpathSignalptr;
		for (k=0; k<MultipathSignalLength; k++) 
		{
			fprintf(fp,"%g %g\n",TempMpathSignalptr->real,TempMpathSignalptr->imaginary);
			TempMpathSignalptr++;
		}
		fclose(fp);

		fp = fopen("FadeSig.txt","w");
		faa = TempFadeptr - AmpBeginPtr + 1;
		for (k=0; k<faa; k++) 
			fprintf(fp,"%g %g\n",(AmpBeginPtr+k)->real,(AmpBeginPtr+k)->imaginary);
		fclose(fp);
		/*******************************************************/

	}
	
	free(CombinedSigPtr);
	return(MpathSignalptr);
}



double ChannelClass::PolynomialFit(double *Xdata,double *Ydata,unsigned N,double x,double *error)
/******************************************************************************************
* double PolynomialFit(double *Xdata,double *Ydata,unsigned N,double x,double *error)
*
*
* Copyright 2002 The Mobile and Portable Radio Research Group
*
* Given the abcissa data Xdata and the ordinate data Ydata (each of length N), PolynomialFit
* estimates the ordinate associcated with the abcissa "x" by fitting the data to a 
* polynomial of degree N-1
*
* Returns a double which is an estimate of the ordinate associated with the abcissa, x.
*
* Parameters
*	Input
*		Xdata	double *	Pointer to the abcissa data
*		Ydata	double *	Pointer to the ordinate data
*		N		unsigned	Length of Xdata and Ydata.  N-1 is also the degree of the polynomial
*		x		double		abcissa value of interest
*	Output
*		dy		double *	Pointer to the variable that stores the error indication
*
* NOTE:  Algorithm was found in Numerical Recipies in C
*******************************************************************************************/
{
	unsigned i,m,ns=1;
	double den,dif,dift,ho,hp,w;
	double *c,*d;
	double *cTemp,*dTemp;
	double *TempXdata,*TempXdata1,*TempYdata;
	double y;

	TempXdata = Xdata;
	TempYdata = Ydata;

	dif = fabs(x - *TempXdata);

	c = (double *) calloc(N,sizeof(double));
	if (c == NULL)
	{
		printf("Memory allocation for pointer c failed--exiting\n");
		return(-1);
	}

	d = (double *) calloc(N,sizeof(double));
	if (d == NULL)