#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "C:\MATLABR11\extern\include\mex.h"


#define CHIPS_PER_FRAME			38400
#define FRAME_DURATION			10e-3
#define ABCISSA_DATA_POINTS		21
#define CHIPS_PER_FADE_SAMPLES	400
#define MIN_FADE_RECORD_LENGTH	512





void MAI_Mpath_Channel(double *RealSigPtr, double *ImagSigPtr,unsigned long NSig,
					   unsigned long *DelayPtr,unsigned long NDelays,double *AmpPtr,
					   double fDoppler,unsigned long SamplesPerChip,
					   double *RealMpathPtr,double *ImagMpathPtr)
/********************************************************************************************************
/void MAI_Mpath_Channel(double *RealSigPtr, double *ImagSigPtr,unsigned long NSig,
/					   unsigned long *DelayPtr,unsigned long NDelays,double *AmpPtr,
/					   double fDoppler,unsigned long SamplesPerChip,
/					   double *RealMpathPtr,double *ImagMpathPtr)
/
/ Copyright 2002 The Mobile and Portable Radio Research Group
/
/This function applies the multipath channel to each interferer.  The function first generates the fading
/signal, through the use of a MATLAB based Rayleigh fading generator.  The fading signal is then applied
/to the interfering signal on a piecewise constant basis.  One fading signal is applied to 400 chiping 
/intervals (see the CHIPS_PER_FADE_SAMPLES constant above).  This is done for each multipath component.  
/Once the fading is applied the multipath component is then offset by its associated multipath delay.  
/
/After all of the multipath compents are scaled by their fading signal and offset by their multipath delay
/then are combined to create the multipath signal for that interferer.
/
/Parameters
/	Input
/		RealSigPtr		*double			Pointer to the real part of the interfering signal
/		ImagSigPtr		*double			Pointer to the imaginary part of the interfering signal
/		NSig			unsigned long	Legnth of the interfereing signal
/		DelayPtr		*unsigned long	Delay of each multipath component in terms of signal samples
/		NDelays			unsigned long	Number of multipath components
/		AmpPtr			*double			Average amplute of each multipath component
/		fDopper			double			Doppler spread in Hz
/		SamplesPerChip	double			Number of samples per chip
/	Output
/		RealMpathPtr	*double			Pointer to the real part of the resulting multipath signal
/		ImagMpathPtr	*double			Pointer to the imaginary part of the resulting multipath signas
**********************************************************************************************************/
{
	double *TempRealSigPtr,*TempImagSigPtr;
	unsigned long *TempDelayPtr;
	double *TempAmpPtr;
	double *TempRealMpathPtr,*TempImagMpathPtr;
	unsigned long i,j,k;
	double SampleRate,SigDuration,FadeSampleRate;
	double increment;
	unsigned long Increment,RemainingSamples;
	mxArray *mxLeftHandSide[1],*mxRightHandSide[3];
	mxArray *mxfDopplerPtr,*mxNAbcissaDataPointsPtr;
	mxArray *mxFadeSampleRatePtr;
	double *FadeSampleRatePtr;
	double *RealFadeSamplesPtr,*ImagFadeSamplesPtr;
	double *TempRealFadeSamplesPtr, *TempImagFadeSamplesPtr;
	double *fDopplerPtr,*NAbcissaDataPointsPtr;
	int ReturnFlag;
	unsigned long FadeSamplesPerFrame,faa,faa1;
	double FadeDelayRatio;
	unsigned long INTFadeDelayRatio;

////////////////////////////////////////////////////////
// Debug Code
//	FILE *fp;
//	double *RealMPathComponentPtr, *TempRealMPathComponentPtr;
//	double *ImagMPathComponentPtr, *TempImagMPathComponentPtr;
// End Debug Code
////////////////////////////////////////////////////////


	SampleRate = (double) CHIPS_PER_FRAME * (double) SamplesPerChip / (double) FRAME_DURATION ;
	SigDuration = (double) NSig / SampleRate;
	FadeSampleRate = (double) CHIPS_PER_FRAME / ((double) FRAME_DURATION * (double) CHIPS_PER_FADE_SAMPLES);
	FadeSamplesPerFrame = CHIPS_PER_FRAME / CHIPS_PER_FADE_SAMPLES;

	//Determine the total number of fade samples available for the signal
	//may not be an integer multiple of the number of samples in the signal 
	//record.  Therefore, it will be necessary  to make the increment (which 
	//is the defined as the number of signal samples per fade sample) take on 
	//next lower integer.  Then we must determine the number of signals samples
	//that are remaining.  These signals will be scaled by the last fading sample

	increment = (double) NSig / ( (double) CHIPS_PER_FRAME / (double) CHIPS_PER_FADE_SAMPLES );
	Increment = (unsigned long) increment;
	RemainingSamples = NSig % Increment;

	if ((mxFadeSampleRatePtr = mxCreateDoubleMatrix(1,1,mxREAL))==NULL)
		mexErrMsgTxt("\nmxFadeSampleRatePtr vector not allocated!--exiting\n");
	FadeSampleRatePtr = mxGetPr(mxFadeSampleRatePtr);
	*FadeSampleRatePtr = FadeSampleRate;

	if ((mxfDopplerPtr = mxCreateDoubleMatrix(1,1,mxREAL))==NULL)
		mexErrMsgTxt("\nmxfDopplerPtr vector not allocated!--exiting\n");
	fDopplerPtr = mxGetPr(mxfDopplerPtr);
	*fDopplerPtr = fDoppler;
	
	if ((mxNAbcissaDataPointsPtr = mxCreateDoubleMatrix(1,1,mxREAL))==NULL)
		mexErrMsgTxt("\nmxNAbcissaDataPointsPtr vector not allocated!--exiting\n");
	NAbcissaDataPointsPtr = mxGetPr(mxNAbcissaDataPointsPtr);
	*NAbcissaDataPointsPtr = (double) CHIPS_PER_FRAME / (double) CHIPS_PER_FADE_SAMPLES;
	if ( *NAbcissaDataPointsPtr < (double) MIN_FADE_RECORD_LENGTH )
		*NAbcissaDataPointsPtr = (double) MIN_FADE_RECORD_LENGTH;

	*mxRightHandSide = mxfDopplerPtr;
	*(mxRightHandSide + 1) = mxFadeSampleRatePtr;
	*(mxRightHandSide + 2) = mxNAbcissaDataPointsPtr;

	TempAmpPtr = AmpPtr;
	TempDelayPtr = DelayPtr;

	for (i = 0; i < NDelays; i++)
	{
		if ( (ReturnFlag = mexCallMATLAB(1,mxLeftHandSide,3,mxRightHandSide,"rayleigh")) != NULL)
			mexErrMsgTxt("\nCall to Matlab Function rayleigh failed\n");

		RealFadeSamplesPtr = mxGetPr(*mxLeftHandSide);
		ImagFadeSamplesPtr = mxGetPi(*mxLeftHandSide);

		TempRealFadeSamplesPtr = RealFadeSamplesPtr;
		TempImagFadeSamplesPtr = ImagFadeSamplesPtr;

		TempRealSigPtr = RealSigPtr + NSig - (int) *TempDelayPtr;
		TempImagSigPtr = ImagSigPtr + NSig - (int) *TempDelayPtr;

		FadeDelayRatio = (double) *TempDelayPtr / ((double) CHIPS_PER_FADE_SAMPLES * (double) SamplesPerChip);
		INTFadeDelayRatio = (unsigned long) floor(FadeDelayRatio);

////////////////////////////////////////////////////
// Debug Code
//		RealMPathComponentPtr = (double *) calloc(NSig,sizeof(double));
//		ImagMPathComponentPtr = (double *) calloc(NSig,sizeof(double));
//		TempRealMPathComponentPtr = RealMPathComponentPtr;
//		TempImagMPathComponentPtr = ImagMPathComponentPtr;
// End Debug Code
////////////////////////////////////////////////////////

		TempRealMpathPtr=RealMpathPtr;
		TempImagMpathPtr=ImagMpathPtr;

		for (j = 0; j < INTFadeDelayRatio; j++)
		{
			for (k = 0; k < Increment; k++)
			{
////////////////////////////////////////////////////
// Debug Code
//				*TempRealMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempRealSigPtr)
//											- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
//				*TempImagMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempImagSigPtr)
//											+ (*TempImagFadeSamplesPtr * *TempRealSigPtr);
// End Debug Code
////////////////////////////////////////////////////////

				*TempRealMpathPtr++ += (*TempRealFadeSamplesPtr * *TempRealSigPtr)
										- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
				*TempImagMpathPtr++ += (*TempRealFadeSamplesPtr * *TempImagSigPtr++)
											+ (*TempImagFadeSamplesPtr * *TempRealSigPtr++);
			}
			TempRealFadeSamplesPtr++;
			TempImagFadeSamplesPtr++;
		}

		faa = *TempDelayPtr - (Increment * INTFadeDelayRatio);
		for (j=0; j < faa; j++)
		{
////////////////////////////////////////////////////
// Debug Code
//			*TempRealMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempRealSigPtr)
//										- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
//			*TempImagMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempImagSigPtr)
//										+ (*TempImagFadeSamplesPtr * *TempRealSigPtr);
// End Debug Code
////////////////////////////////////////////////////////

			*TempRealMpathPtr++ += (*TempRealFadeSamplesPtr * *TempRealSigPtr)
									- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
			*TempImagMpathPtr++ += (*TempRealFadeSamplesPtr * *TempImagSigPtr++)
									+ (*TempImagFadeSamplesPtr * *TempRealSigPtr++);
		}

		TempRealSigPtr = RealSigPtr;
		TempImagSigPtr = ImagSigPtr;
		for (j = faa; j < Increment; j++)
		{
////////////////////////////////////////////////////
// Debug Code
//			*TempRealMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempRealSigPtr)
//										- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
//			*TempImagMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempImagSigPtr)
//										+ (*TempImagFadeSamplesPtr * *TempRealSigPtr);
// End Debug Code
////////////////////////////////////////////////////////

			*TempRealMpathPtr++ += (*TempRealFadeSamplesPtr * *TempRealSigPtr)
									- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
			*TempImagMpathPtr++ += (*TempRealFadeSamplesPtr * *TempImagSigPtr++)
									+ (*TempImagFadeSamplesPtr * *TempRealSigPtr++);
		}
		TempRealFadeSamplesPtr++;
		TempImagFadeSamplesPtr++;

		faa1 = INTFadeDelayRatio +1;
		faa=(FadeSamplesPerFrame-1);
		for (j = faa1; j < faa; j++)
		{
			for (k = 0; k < Increment; k++)
			{
////////////////////////////////////////////////////
// Debug Code
//				*TempRealMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempRealSigPtr)
//											- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
//				*TempImagMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempImagSigPtr)
//											+ (*TempImagFadeSamplesPtr * *TempRealSigPtr);
// End Debug Code
////////////////////////////////////////////////////////

				*TempRealMpathPtr++ += (*TempRealFadeSamplesPtr * *TempRealSigPtr)
										- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
				*TempImagMpathPtr++ += (*TempRealFadeSamplesPtr * *TempImagSigPtr++)
										+ (*TempImagFadeSamplesPtr * *TempRealSigPtr++);
			}
			TempRealFadeSamplesPtr++;
			TempImagFadeSamplesPtr++;
		}

		for (j = 0; j < (Increment + RemainingSamples); j++)
		{
////////////////////////////////////////////////////
// Debug Code
//			*TempRealMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempRealSigPtr)
//										- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
//			*TempImagMPathComponentPtr++ = (*TempRealFadeSamplesPtr * *TempImagSigPtr)
//										+ (*TempImagFadeSamplesPtr * *TempRealSigPtr);
// End Debug Code
////////////////////////////////////////////////////////

			*TempRealMpathPtr++ += (*TempRealFadeSamplesPtr * *TempRealSigPtr)
									- (*TempImagFadeSamplesPtr * *TempImagSigPtr);
			*TempImagMpathPtr++ += (*TempRealFadeSamplesPtr * *TempImagSigPtr++)
									+ (*TempImagFadeSamplesPtr * *TempRealSigPtr++);
		}
////////////////////////////////////////////////////
// Debug Code
//		fp=fopen("fade.txt","w");
//		for (k = 0; k < *NAbcissaDataPointsPtr; k++) 
//			fprintf(fp,"%g %g \n",*(RealFadeSamplesPtr+k),*(ImagFadeSamplesPtr+k));
//		fclose(fp);
//	
//		fp=fopen("BeforeFade.txt","w");
//		for (k=0; k < NSig; k++)
//			fprintf(fp,"%g %g \n",*(RealSigPtr+k),*(ImagSigPtr+k));
//		fclose(fp);
//
//		fp=fopen("AfterFade.txt","w");
//		for (k=0; k < NSig; k++)
//			fprintf(fp,"%g %g\n",*(RealMPathComponentPtr+k),*(ImagMPathComponentPtr+k));
//		fclose(fp);
//		free(RealMPathComponentPtr);
//		free(ImagMPathComponentPtr);
// End Debug Code
////////////////////////////////////////////////////////

		TempAmpPtr++;
		TempDelayPtr++;


	}

////////////////////////////////////////////////////
// Debug Code
//	fp=fopen("MPathSig.txt","w");
//	for (k=0; k < NSig; k++)
//		fprintf(fp,"%g %g\n",*(RealMpathPtr + k),*(ImagMpathPtr + k));
//	fclose(fp);
// End Debug Code
////////////////////////////////////////////////////////

	mxDestroyArray(mxFadeSampleRatePtr);
	mxDestroyArray(mxfDopplerPtr);
	mxDestroyArray(mxNAbcissaDataPointsPtr);
}