wcdmalinearinterp.cpp

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

#include <malloc.h>
#include "mex.h"
//#include <math.h>
#include "matrix.h"
#include <stdlib.h>

void Linear(double *,double *,double *,unsigned,
		   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[])
*
* Gateway function for WCDMALinearInterp.  Supports that passing of three input parameters
* and one output parameter.  The five input parameters are
*		Data Abcissa Values			Array containing the abcissa values associated 
*									with the data to be interpolated
*		DataValues					Array containing the data to be interpolated
*		Interpolated Abcissa Values	Array containing the abcissa values of points to 
*									be interpolated
*
* The output parameter is
*		Interpolated Curve			Interpolated data points
***************************************************************************************/
{
	const mxArray *tprhs;
	unsigned mexDataAbcissaLength,mexDataLength,mexAbcissaLength,mrows,ncols;
	double *mexDataAbcissaValues,*mexRealDataValues,*mexImagDataValues;
	double *mexAbcissaValues;
	double *mexRealInterp,*mexImagInterp;
	
	
	//Check for the proper number of input and output arguments
	if (nrhs != 3) mexErrMsgTxt("\nExactly THREE input arguemnts are required!!\n");
	else if (nlhs != 1 ) mexErrMsgTxt("\nExactly ONE output arguement is required!!\n");

	
	//Data Abcissa Values
	//Must be a double vector
	tprhs = *prhs;
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || (mrows ==1 && ncols ==1) || (mrows >1 && ncols >1))
		mexErrMsgTxt("\nX must be a real-valued vector\n");
	if (mrows > ncols) mexDataAbcissaLength = (unsigned ) mrows;
	else mexDataAbcissaLength = (unsigned) ncols;
	mexDataAbcissaValues = mxGetPr(tprhs);
	
	

	//Data Values
	//Must be a double, complex and a vector
	tprhs = *(prhs+1);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || !mxIsComplex(tprhs) || (mrows ==1 && ncols ==1) || (mrows >1 && ncols >1))
		mexErrMsgTxt("\nY must be a complex vector\n");
	if (mrows > ncols) mexDataLength = (unsigned) mrows;
	else mexDataLength = (unsigned) ncols;
	if (mexDataLength != mexDataAbcissaLength)
		mexErrMsgTxt("\nBoth X and Y must have the same length!\n");
	mexRealDataValues = mxGetPr(tprhs);
	mexImagDataValues = mxGetPi(tprhs);

	//Interpolated Abcissa Values
	//Must be a double vector
	tprhs = *(prhs+2);
	mrows = mxGetM(tprhs);
	ncols = mxGetN(tprhs);
	if ( !mxIsDouble(tprhs) || mxIsComplex(tprhs) || (mrows ==1 && ncols ==1) || (mrows >1 && ncols >1))
		mexErrMsgTxt("\nXI must be a real-valued vector\n");
	if (mrows > ncols) mexAbcissaLength = (unsigned) mrows;
	else mexAbcissaLength = (unsigned) ncols;
	mexAbcissaValues = mxGetPr(tprhs);
	
	

	//Allocate output Array for interpolated curve
	*plhs = mxCreateDoubleMatrix(mexAbcissaLength,1,mxCOMPLEX);
	mexRealInterp = mxGetPr(*plhs);
	mexImagInterp = mxGetPi(*plhs);

	//Call Linear interpoaltor
	Linear(mexDataAbcissaValues,mexRealDataValues,mexImagDataValues,mexDataAbcissaLength,
		  mexAbcissaValues,mexRealInterp,mexImagInterp,mexAbcissaLength);

}

void Linear(double *x,double *realy,double *imagy,unsigned Ny,
		   double *xi, double *Realyi,double *Imagyi,unsigned Nyi)
/*********************************************************************************************
*void Linear(double *x,double *realy,double *imagy,unsigned Ny,
*		   double *xi, double *Realyi,double *Imagyi,unsigned Nyi)
*
* Performs linear interpolation.This particular funciton is a scaled down version of the 
* corresponding MATLAB function because many of the error checking features have been removed
*
* Parameters
*	Input
*		*x		Pointer to the vector containing the abcissa values for the data
*		*realy	Pointer to the array contain the real part of the data
*		*imagy	Pointer to the array contain the imaginary part of the data
*		Ny		Lenght of the data array (i.e., *x, *realy, *imagy)
*		*xi		Pointer to the array contianing the abcissa values for the 
*				interpolated sequence
*		Nyi		Lenght of the interpolated array
*	Output
*		*Realyi	Pointer to the array containing the real part of the interpolated sequence
*		*Imagyi	Pointer to the array containing the imaginary part of the interpolated 
*				sequence
*********************************************************************************************/
{
	unsigned long i;
	double *TempAbcissaI;  //Pointers for *xi
	double *TempDataAbcissa;  //Pointers for *x 
	double *TempRealData,*TempImagData;  //Temporary Pointers for *realy, *imagy
	double *TempRealyi,*TempImagyi;  //Temporary Pointers for *Realyi,Imagyi
	double scale,s;
	double *ScaledAbcissa,*ScaledAbcissaTemp;  //Pointers for *ScaledAbcissa
	int *FloorAbcissa,*TempFloorAbcissa;   //Pointer to array that contains the integer portion 
	                                        //of the Scaled Abcissa
	double *DecimalAbcissa,*TempDecimalAbcissa; //Pointer to array that contains the decimal portion 
	                                            //of the Scaled Abcissa

	ScaledAbcissa = (double *) mxCalloc(Nyi,sizeof(double));
	if (ScaledAbcissa == NULL) mexErrMsgTxt("\nMemory Allocation failed for ScaledAbcissa\n");

	FloorAbcissa = (int *) mxCalloc(Nyi,sizeof(int));
	if (FloorAbcissa == NULL) mexErrMsgTxt("\nMemory Allocation failed for ScaledAbcissa\n");

	DecimalAbcissa = (double *) mxCalloc(Nyi,sizeof(double));
	if (DecimalAbcissa == NULL) mexErrMsgTxt("\nMemory Allocation failed for ScaledAbcissa\n");

	//Scale and shift abcissa to index points
	TempAbcissaI = xi;
	TempDataAbcissa = x;
	ScaledAbcissaTemp=ScaledAbcissa;
	
	scale = (double) (Ny-1) / (*(x+Ny-1) - *x);
	for (i=0 ; i<Nyi; i++)
	{
		*ScaledAbcissaTemp++ = (*TempAbcissaI++ - *TempDataAbcissa) * scale;
	}
    TempFloorAbcissa=FloorAbcissa;
	TempDecimalAbcissa=DecimalAbcissa;
	ScaledAbcissaTemp=ScaledAbcissa;
	for (i=0;i<Nyi;i++)
	{
		*TempFloorAbcissa = (int) *ScaledAbcissaTemp;
		*TempDecimalAbcissa++ = *ScaledAbcissaTemp++ - (double) *TempFloorAbcissa++;
	}

    TempFloorAbcissa=FloorAbcissa;
	TempRealyi = Realyi;
	TempImagyi = Imagyi;
	TempDecimalAbcissa=DecimalAbcissa;
	for (i=0; i<Nyi; i++)
	{
		s = *TempDecimalAbcissa++;
		TempRealData = realy + *TempFloorAbcissa;
		TempImagData = imagy + *TempFloorAbcissa++;

		*TempRealyi = *TempRealData++ * ( (double) 1.0 -s); 
		*TempRealyi = *TempRealyi + *TempRealData * s; 
		*TempRealyi++;
		*TempImagyi = *TempImagData++ * ( (double) 1.0 -s); 
		*TempImagyi = *TempImagyi + *TempImagData * s; 
		*TempImagyi++;

	}
	mxFree(ScaledAbcissa);
	mxFree(FloorAbcissa);
	mxFree(DecimalAbcissa);

}