#include <math.h>
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
#include "mex.h"

#define PI 3.14159265358979

void GenRootRaisedCosine(int , int ,unsigned ,double ,double ,double * ,double * ,unsigned );

void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/**************************************************************************************
* void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
*
* Copyright 2002 The Mobile and Portable Radio Research Group
*
* Gateway function for GenRootRaisedCosine.  Supports that passing of five input parameters
* and two output parameters.  The five input parameters are
*		LBound				Determines the minimum limit of the truncated 
*							impulse response normalized by symbol period
*		UBound				Determines the maximum limit of the truncated 
*							impulse response normalized by symbol period
*		SampleRate			The number of samples per symbol duration
*		RollOff				Rolloff of the Filter.  Must not exceed 1
*		SymbolPeriod		Symbol Duration (in seconds)
*
* The two output parameters are
*		y					Impulse response
*		t					Time index
***************************************************************************************/
{
	int MexLBound,MexUBound,PLength;
	unsigned MexSamplesPerSymbol;
	double MexRollOff,MexSymbolPeriod;
	int mrows,ncols;
	double dbleScale,fraction,integer_portion;
	double *pulse,*TimeIndex;
	const mxArray *tprhs;
	
	
	//Check for the proper number of input and output arguments
	if (nrhs != 5) mexErrMsgTxt("\nExactly FIVE input arguemnts are required!!\n");
	else if (nlhs != 2 ) mexErrMsgTxt("\nExactly TWO output arguements are required!!\n");

	//All Inputs must be scalars
	//Lower Bound
	tprhs = *prhs;
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || !(mrows ==1 && ncols ==1))
		mexErrMsgTxt("\nLowerBoundT must be a noncomplex scalar\n");
	MexLBound = (int) mxGetScalar(tprhs);


	//Upper Bound
	tprhs = *(prhs+1);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || !(mrows ==1 && ncols ==1))
		mexErrMsgTxt("\nUpperBoundT must be a noncomplex scalar\n");
	MexUBound = (int) mxGetScalar(tprhs);

	//Samples per Symbol
	tprhs = *(prhs+2);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || !(mrows ==1 && ncols ==1))
		mexErrMsgTxt("\nSamples Per Symbol must be a noncomplex scalar\n");
	dbleScale=mxGetScalar(tprhs);
	fraction=modf(dbleScale,&integer_portion);
	if (integer_portion <= 0.0) mexErrMsgTxt("\nMexSamplesPerSymbol must be positive\n");
	MexSamplesPerSymbol = (unsigned) integer_portion;
	if (fraction != 0) mexErrMsgTxt("\nMexSamplesPerSymbol must be an integer\n");

	//RollOff
	tprhs = *(prhs+3);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || !(mrows ==1 && ncols ==1))
		mexErrMsgTxt("\nSamples Per Symbol must be a noncomplex scalar\n");
	MexRollOff=mxGetScalar(tprhs);
	if ((MexRollOff <= 0.0) || (MexRollOff > 1.0)) mexErrMsgTxt("\nMexRollOff must be between 0 and 1\n");

	//Symbol Period
	tprhs = *(prhs+4);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || !(mrows ==1 && ncols ==1))
		mexErrMsgTxt("\nSamples Per Symbol must be a noncomplex scalar\n");
	MexSymbolPeriod=mxGetScalar(tprhs);
	if (MexRollOff <= 0.0)  mexErrMsgTxt("\nSample Period must be positive\n");

	PLength = (MexUBound - MexLBound) * (int) MexSamplesPerSymbol + 1;
	*plhs = mxCreateDoubleMatrix( PLength,1,mxREAL);
	*(plhs+1) = mxCreateDoubleMatrix( PLength,1,mxREAL);
	pulse = mxGetPr(*plhs);
	TimeIndex = mxGetPr(*(plhs+1));
	//Call Root Raised Cosine function
	GenRootRaisedCosine(MexLBound, MexUBound,MexSamplesPerSymbol,MexRollOff, MexSymbolPeriod,pulse,TimeIndex,PLength);
	

}





void GenRootRaisedCosine(int LowerBoundT, int UpperBoundT,unsigned SamplesPerSymbol,double RollOff,double SymbolPeriod,double *filter,double *t_index,unsigned NumSamples)
/******************************************************************************************
/ void GenRootRaisedCosine(int LowerBoundT, int UpperBoundT,unsigned SamplesPerSymbol,
/      double RollOff,double SymbolPeriod,double *filter,double *t_index,unsigned NumSamples)
/
/ Copyright 2002 The Mobile and Portable Radio Research Group
/
/
/ GenRootRaisedCosine produces the impulse response of a root raised cosine filter and
/ stores it in an array.  Note that the impulse of a root raised cosine filter ranges
/ from -infinity to infinity.  It is therefore necessary to truncated from LowerBoundT to
/ UperBountT.
/
/ Parameters
/   Input
/      LowerBoundT        Determines the minimum limit of the truncated imulse response
/                         normalized by symbol period
/      UpperBoundT        Determines the maximum limit of the truncated imulse response
/                         normalized by symbol period
/      SamplesPerSymbol   The number of samples per symbol duration
/      RollOff            Rolloff of the Filter.  Must not exceed 1
/      SymbolPeriod       Symbol Duration (in seconds)
/ 
/   Output
/      *filter            Pointer ot an array that contains the filter impulse response
/      *t_index           Pointer ot an array that contains the associated time.
/      NumSamples         Length of the output array
/
/******************************************************************************************/
{
	double *filter_temp,*t_temp;	//Temporary pointers for *filter and *t_index
	double cos_arg;					//Stores arguement to cosine funciton
	double sin_arg;					//Stores arguement to sine function
	double bottom_arg;				//Stores non-trivial term in the denominator
	double coeff;					//Stores the normalizing coefficient
	double SampleDuration;			//Stores the sample period
	int k;							//Counter

	//Determine Array length
	NumSamples = (UpperBoundT-LowerBoundT)*SamplesPerSymbol+1;

	
	//Compute impulse response
	SampleDuration = SymbolPeriod/SamplesPerSymbol;
	coeff = (4.0 * RollOff * sqrt(SampleDuration)) / (PI * sqrt(SymbolPeriod));

    t_temp=t_index;
	filter_temp=filter;
	for (k=0; k< (int) NumSamples;k++)
	{
		//Compute time index
		*t_temp=(LowerBoundT * SymbolPeriod) + (k * SampleDuration);
		if (*t_temp == double (0.0)) *t_temp = 1e-20;
		cos_arg = (1 + RollOff) * (PI * *t_temp / SymbolPeriod);
		sin_arg = (1 - RollOff) * (PI * *t_temp / SymbolPeriod);
		bottom_arg = (double) 4.0 * RollOff * *t_temp / SymbolPeriod;
		t_temp++;
		*filter_temp++ = coeff * (cos(cos_arg) + (sin(sin_arg) / bottom_arg)) / (1 - (bottom_arg * bottom_arg));
	}
}