/*------------------------------------------------------------------------------*
 * File Name:				 													*
 * Creation: 																	*
 * Purpose: OriginC Source C file												*
 * Copyright (c) ABCD Corp.	2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010		*
 * All Rights Reserved															*
 * 						
 *
 * Modification Log:															*
 *		Soapy 12/30/03 CHANGE_RETURN_ON_TIME_RESOLUTION							*
 *		Soapy 12/30/03 ADD_THRESHOLD_FILTERING_TO_DLL							*
 *		Soapy 12/30/03 MODIFICATION_ON_STFT										*
 *		Annie 3/1/04 SAMPLING_RESOLUTION_PARA_MODIFICATION						*
 *		Annie 3/4/04 STFT_MODIFICATION											*
 *		Annie 3/18/04 STFT_INTERFACE_MODIFICATION								*
 *		Annie 3/19/04 WINDOW_FUNCTION_CREATE									*
 *		Iris 04/01/04 ADD_PREFIX_FOR_METHOD_NAME								*
 *		Annie 4/15/04 STFT_INTERFACE_UPDATE										*
 *		Annie 4/15/04 CHECK_VECTOR												*
 *		Annie 4/15/04 STFT_OVERLOAD_ROUTINE										*
 *		Annie 4/16/04 FFT_ROUTINE_ADDING										*
 *		Annie 4/22/04 FFT_RELATED_OVERLOAD										*
 *		Annie 4/23/04 WINDOW_SIZE1_MODIFICATION									*
 *		Annie 4/23/04 CORRELATE_RETURN_VALUE_ADD								*
 *		Annie 4/23/04 HIGH_LEVEL_ROUTINE_INTERFACE_MODIFICATION					*
 *		Annie 4/23/04 NORMALIZE_MODIFICATION									*
 *		Annie 4/23/04 CLEAR_ADD_PREFIX_FOR_METHOD_NAME							*
 *		Annie 4/23/04 ROBUSTNESS_CHECK											*
 *		Annie 4/26/04 RETURN_VALUE_UNIFICATION									*
 *		Annie 4/26/04 FFT_SAMPLING_RESOLUTION_INTERFACE_UPDATE					*
 *		AW 05/21/04 QA70-6387 v7.0876 CENTRALIZE_CODE_FOR_STFT					*
 *		AW 05/21/04 QA70-6387 v7.0876 WILL_LEAD_CRASH_HERE						*
 *  	AW Soapy 05/28/04 QA75-6457 v7.5847 CLEANUP_STFT						*
 *		AW 05/28/04 QA75-6456 v7.05848 CLEANUP_STFT_HANDLE_ERROR_CODES			*
 *		AW 06/29/04 MORE_ON_STFT_SIZE											*
 *		Forest 08/16/04 RETURN_VALUE_UNIFICATION                                *
 *      Forest 08/16/04 STFT_INTERFACE_UPDATE									*
 *		Forest 08/27/04 MODIFY_ERROR_HANDING									*
 *		Forest 09/03/04 DEAL_WITH_SAME_VECTOR                                   *
 *		Forest 09/21/04 QA70-6227 CHANGE_PROTOTYPE								*
 *		ER 05/27/05 CHANGE_FUNCTION_NAME										*
 *		Leo 06/04/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH							*
 *		Soapy 07/03/05 FIX_NORMALIZATION_FUNCTION								*
 *		Soapy 07/03/05 ADD_CORRELATION_OPTION_LINEAR_CIRCULAR					*
 *		Soapy 07/10/05 QA70-7865 ADD_2D_CORRELATION								*
 * 		Soapy 07/10/05 QA70-5060 ADD_COHERENCE_FUNCTION							*
 *		Leo 09/19/05 QA70-8103 REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT				*
 *		Austin/Leo 09/20/05	QA70-8103 REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT		*
 *		Leo 2005-9-22 v8.0340 REWRITE_CHECK_CURVE 								*
 *		Leo/Raine 2006-2-16 LINEAR_OR_CIRCULAR_FOR_CONVOLUTION					*
 *		Raine/Leo 09/03/06 ADD_WELCH_WIN										*
 *		Raine/Leo 2006-3-10 STFT_RESULT_IS_DIFFERENT_FROM_AUTOSIGNAL			*
 *		Raine/Leo 2006-3-16 REWRITE_FFT_FFT_COMPLEX								*
 *		Leo 3006-3-21 FIX_COMPLEX_VECTOR_POINTER_BUG							*
 *		Leo 2006-3-31 FFTW_FFT_2D_COMPLEX										*
 *		Raine 2006-04-06 ADD_STFT_COMPLEX										*
 *		Raine 2006-04-13 CHANGE_STFT_PARAMETER									*
 *		Raine 2006-04-20 CORRECT_CORR2_DIMENSION								*
 *      Sandy 2006-5-22  CLEAN_UP_STFT_RELATED_FUNCTION							*
 * 		Arvin 2006-5-23 ADD_GET_WINDOW_DATA										*
 * 		Arvin 2006-5-22 ADD_GET_FREQUENCE_FROM_INTERNAL							*
 * 		Arvin 2006-5-22 ADD_FFT_NORM_AMP										*
 * 		Arvin 2006-5-22 ADD_FFT_FACTOR											*
 * 		ARVIN 2006-5-22 ADD_FFT_SHIFT											*
 * 		Arvin 2006-5-22 ADD_FFT_UNWRAP_PHASE									*
 * 		Arvin 2006-5-10 ADD_APPLY_WINDOW_DATA									*
 *------------------------------------------------------------------------------*/
 
////////////////////////////////////////////////////////////////////////////////////
// Including the system header file Origin.h should be sufficient for most Origin
// applications and is recommended. Origin.h includes many of the most common system
// header files and is automatically pre-compiled when Origin runs the first time.
// Programs including Origin.h subsequently compile much more quickly as long as
// the size and number of other included header files is minimized. All NAG header
// files are now included in Origin.h and no longer need be separately included.
//
// Right-click on the line below and select 'Open "Origin.h"' to open the Origin.h
// system header file.
#include <Origin.h>
#include <..\originlab\FFT.h>
#include <fft_utils.h> /// AW 05/21/04 QA70-6387 v7.0876 CENTRALIZE_CODE_FOR_STFT
#include <ReportTree.h>
#include <fftEx_utils.h>
#include <oErrMsg.h>
//#include <matrix.h>

#pragma labtalk(0) //--- CPY 7/26/05 hide all functions in this file from LT access, as XF should be used instead

enum FFT_FILTER_TYPE{
	FFT_LOWPASS = 0, FFT_HIGHPASS, FFT_BANDPASS, FFT_BANDBLOCK
};


/*
enum {
	FFT_REPORT_TABLE = 101,
};

enum {
	IDST_FFT_FREQ = 0x00010010,
	IDST_FFT_REAL,
	IDST_FFT_IMAG,
	IDST_FFT_AMP,
	IDST_FFT_PHASE,
	IDST_FFT_POWER,
};
*/


/*
/// Annie 4/23/04 CLEAR_ADD_PREFIX_FOR_METHOD_NAME
/// Annie 3/19/04 WINDOW_FUNCTION_CREATE
int RectWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return fft_error_handle(INVALID_SIZE);
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION
	
	int nRet = stft_RectWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int TriWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return fft_error_handle(INVALID_SIZE);
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_TriWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int BartlettWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return fft_error_handle(INVALID_SIZE);
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_BartlettWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int HanningWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return INVALID_SIZE;
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_HanningWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int HammingWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return INVALID_SIZE;
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_HammingWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int BlackmanWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return INVALID_SIZE;
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_BlackmanWin(vWin, iSize);
	if (nRet)
		return nRet;
	return 0;
}

int GaussWin(vector &vWin, int iSize, double dAlpha)
{
	if (iSize<=0)
		return INVALID_SIZE;
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_GaussWin(vWin, iSize, dAlpha);
	if (nRet)
		return nRet;
	return 0;
}

int KaiserWin(vector &vWin, int iSize, double dBeta)
{
	if (iSize<=0)
		return fft_error_handle(INVALID_SIZE);
	vWin.SetSize(iSize);
	///Annie 4/23/04 WINDOW_SIZE1_MODIFICATION
	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	///End WINDOW_SIZE1_MODIFICATION

	int nRet = stft_KaiserWin(vWin, iSize, dBeta);
	if (nRet)
		return nRet;
	return 0;
}
/// End WINDOW_FUNCTION_CREATE
/// End CLEAR_ADD_PREFIX_FOR_METHOD_NAME

/// Raine/Leo 09/03/06 ADD_WELCH_WIN
int WelchWin(vector &vWin, int iSize)
{
	if (iSize<=0)
		return fft_error_handle(INVALID_SIZE);
	vWin.SetSize(iSize);

	if (iSize==1)
	{
		if(SetSize1(vWin)) 
			return fft_error_handle(UNKNOWN_ERROR);
		return 0;
	}
	
	int nRet = stft_WelchWin(vWin, iSize);
	
	if (nRet)
	{
		return nRet;
	}

	return 0;
}
*/

int fft_2d_complex(int nFFTRow, int nFFTCol, matrix<complex>& mSig, FFT_SIGN iSign)
{
	///Forest 08/16/04 RETURN_VALUE_UNIFICATION 
	if(nFFTRow < 1 || nFFTCol < 1)
		//return fft_error_handle(INVALID_SIZE);
		return CER_INVALID_SIZE;
	///End RETURN_VALUE_UNIFICATION 
	
	//int nRet = check_vector_2d(mSig); ///Leo 06/04/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH	
	//if (nRet)	return nRet;
	
	int nRow = mSig.GetNumRows();
	int nCol = mSig.GetNumCols();
	if(nFFTRow != nRow || nFFTCol != nCol)
		mSig.SetSize(nFFTRow,nFFTCol, true);
	
	matrix mReal, mImag;
	///Cheney 2006-7-14 CENTRALIZE_CODE
	//mSig.GetReal(mReal);
	//mSig.GetImaginary(mImag);
	
	///Sandy 2006-11-29 REMOVE_MISS_CHECK
	//int nRet = _check_matrix_real_and_imag(mSig, mReal, mImag, nFFTRow, nFFTCol);
	//if(nRet != 0)
		//return nRet;
	///end REMOVE_MISS_CHECK
	
	///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	///Leo 06/04/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH
	//int nRet = ocmath_count_missing_values(nFFTRow*nFFTCol, mReal);
	//int nRet1 = ocmath_count_missing_values(nFFTRow*nFFTCol, mImag);
	//int nRet = ocmath_count(NANUM, nFFTRow*nFFTCol, mReal);
	//int nRet1 = ocmath_count(NANUM, nFFTRow*nFFTCol, mImag);
	//if (nRet || nRet1)
	//	return INVALID_DATA;
	///End REPLACE_CHECK_VECTOR_2D_BY_OCMATH	
	///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	///end CENTRALIZE_CODE
	
	///Leo 2006-3-31 FFTW_FFT_2D_COMPLEX
	//nRet = fft_fft_2d_complex(nFFTRow, nFFTCol, mReal, mImag, iSign);
	//mSig.MakeComplex(mReal,mImag);
	int nRet = fftw_fft_2d_complex(nFFTRow, nFFTCol, mSig, iSign);
	return nRet;
}

// int fft_real_2d(int nFFTRow, int nFFTCol, matrixbase& mSig, FFT_SIGN iSign)
// ER 05/27/05 CHANGE_FUNCTION_NAME
// Changed name to be compatible with fft_2d_complex
// int fft_2d_real(int nFFTRow, int nFFTCol, matrixbase& mSig, FFT_SIGN iSign)
// Leo 6/28/05 CHANGE_FUNCTION_PROTOTYPE
int fft_2d_real(int nFFTRow, int nFFTCol, matrix<complex>& mSig, FFT_SIGN iSign)
{
	//int nRet = check_vector_2d(mSig);	///Leo 06/04/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH
	//if (nRet)	return nRet;
	if(nFFTRow < 1 || nFFTCol < 1)
		return CER_INVALID_SIZE;
		
	int nRet;
	int nRow = mSig.GetNumRows();
	int nCol = mSig.GetNumCols();
	if(nFFTRow!=nRow || nFFTCol!=nCol)
		mSig.SetSize(nFFTRow,nFFTCol,true);
	
	matrix mReal, mImag;
	if(iSign==FFT_FORWARD)
	{
		mSig.GetReal(mReal);

		///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
		///Leo 06/28/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH	
		//nRet = ocmath_count_missing_values(nFFTRow*nFFTCol, mReal);
		
		///Sandy 2006-11-25 NO_NEED_TO_CHECK_MISS_FOR_2D
		//if (nRet = ocmath_count(NANUM, nFFTRow*nFFTCol, mReal))
			//return CER_INVALID_DATA;
		///end NO_NEED_TO_CHECK_MISS_FOR_2D
		///End REPLACE_CHECK_VECTOR_2D_BY_OCMATH	
		///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT

		mImag.SetSize(nFFTRow,nFFTCol);
		nRet = fft_fft_2d_real(nFFTRow, nFFTCol, mReal, mImag, FFT_FORWARD);
		mSig.MakeComplex(mReal,mImag);
	}
	///Forest 08/27/04 MODIFY_ERROR_HANDING
	//else
	else if(iSign == FFT_BACKWARD)
	{
		///Sandy 2007-8-27 RockBack Cheney's remove, it is source data
		//Cheney 2006-7-14 CENTRALIZE_CODE
		mSig.GetReal(mReal);
		mSig.GetImaginary(mImag);
		
		///Sandy 2006-11-25 NO_NEED_TO_CHECK_MISS_FOR_2D
		//nRet = _check_matrix_real_and_imag(mSig, mReal, mImag, nFFTRow, nFFTCol);
		//if(nRet != 0)
			//return nRet;
		///end 2006-11-25 NO_NEED_TO_CHECK_MISS_FOR_2D
		///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
		///Leo 06/28/05 REPLACE_CHECK_VECTOR_2D_BY_OCMATH
		//nRet = ocmath_count_missing_values(nFFTRow*nFFTCol, mReal);
		//int nRet1 = ocmath_count_missing_values(nFFTRow*nFFTCol, mImag);
		//nRet = ocmath_count(NANUM, nFFTRow*nFFTCol, mReal);
		//int nRet1 = ocmath_count(NANUM, nFFTRow*nFFTCol, mImag);
		//if (nRet || nRet1)
		//	return INVALID_DATA;
		///End REPLACE_CHECK_VECTOR_2D_BY_OCMATH	
		///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
		///end CENTRALIZE_CODE

		nRet = fft_fft_2d_real(nFFTRow, nFFTCol, mReal, mImag, FFT_BACKWARD);
		mSig = mReal;
	}
	else
		return	CER_INVALID_ARGUMENT;
	///End MODIFY_ERROR_HANDING
	return nRet;
}
///End FFT_ROUTINE_ADDING

///Cheney 2006-7-14 ADD_STATIC_FUNC_FOR_FFT_2D_COMPLEX_AND_REAL
static int _check_matrix_real_and_imag(matrix<complex>& mSrc, matrix& mReal, matrix& mImag, int nFFTRow, int nFFTCol)
{
	mSrc.GetReal(mReal);
	mSrc.GetImaginary(mImag);
	
	int nRet = ocmath_count(NANUM, nFFTRow*nFFTCol, mReal);
	int nRet1 = ocmath_count(NANUM, nFFTRow*nFFTCol, mImag);
	if (nRet || nRet1)
		return CER_INVALID_DATA;
	
	return 0;
}
///add ADD_STATIC_FUNC_FOR_FFT_2D_COMPLEX_AND_REAL

///Soapy 07/10/05 QA70-7865 ADD_2D_CORRELATION
///Cheney 2006-7-14 CHANGE_FUNCTION_TO_STATIC
///Sandy 2008-1-21 MATRXBASE_SHOULD_BE_USED
//static int _check_2d_correlation_matrices(matrix& matSignal1, matrix& matSignal2, int& nRows, int& nCols)
static int _check_2d_correlation_matrices(matrixbase& matSignal1, matrixbase& matSignal2, int& nRows, int& nCols)
///end CHANGE_FUNCTION_TO_STATIC
{
	int nRows1, nCols1;
	nRows1 = matSignal1.GetNumRows();
	nCols1 = matSignal1.GetNumCols();
	nRows = matSignal2.GetNumRows();
	nCols = matSignal2.GetNumCols();
	
	//Check the validy of the matrix data
	if(nRows1<1 || nRows<1 || nCols1<1 || nCols<1)
		return CER_INVALID_SIZE;
	
	///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//if(ocmath_count_missing_values(nRows1*nCols1, matSignal1))
	//if(ocmath_count_missing_values(nRows*nCols, matSignal2))
	//if ( ocmath_count(NANUM, nRows1*nCols1, matSignal1))
		//return CER_INVALID_DATA;
	//if ( ocmath_count(NANUM, nRows*nCols, matSignal2))
		//return CER_INVALID_DATA;
	///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
		
	// Get length of signal and pad with zeroes for the shorter one
	nRows = max(nRows, nRows1);
	nCols = max(nCols, nCols1);
	
	return 0;
}
	
int fft_2d_correlation(matrixbase& matSignal1, matrixbase& matSignal2, matrix& matCorrelation, bool bNormalize, bool bCircular)
{
	int iRet;	

	int nRows, nCols;	
	int mat1Rows = matSignal1.GetNumRows(), mat1Cols = matSignal1.GetNumCols();	
	
	///Sandy 2006-11-7 ADD_SUPPORT_FOR_MORE_MATRIX_TYPE
	matrix mat1 = matSignal1;
	matrix mat2 = matSignal2;
	//iRet = _check_2d_correlation_matrices(matSignal1, matSignal2, nRows, nCols);
	iRet = _check_2d_correlation_matrices(mat1, mat2, nRows, nCols);
	//iRet = check_2d_correlation_matrices(matSignal1, matSignal2, maxRows, maxCols);
	if(iRet)
		//return fft_error_handle(iRet);
		return iRet;
		
	if(!bCircular)
	{
		///Raine 2006-04-20 CORRECT_CORR2_DIMENSION
		//nRows = nRows * 2 - 1;
		//nCols = nCols * 2 - 1;
		nRows = matSignal1.GetNumRows() + matSignal2.GetNumRows() - 1;
		nCols = matSignal1.GetNumCols() + matSignal2.GetNumCols() - 1;
	}
	
	matrix mSignal1(nRows,nCols), mSignal2(nRows,nCols);
	mSignal1.SetSubMatrix(matSignal1, 0, 0);
	mSignal2.SetSubMatrix(matSignal2, 0, 0);
		
	if(bNormalize) 	//normalize
	{
		_fft_normalize_response(mSignal1);
		_fft_normalize_response(mSignal2);
	}

	matCorrelation = mSignal1;
	if((iRet = fft_fft_2d_correlation(nRows, nCols, matCorrelation, mSignal2)) != 0) 
		//return fft_error_handle(iRet);
		return iRet;
	
	if(!bCircular)
	{
		//matCorrelation.Wrap((nRows+1)/2, (nCols+1)/2);
		///Raine 2006-04-20 CORRECT_CORR2_DIMENSION
		matCorrelation.Wrap(-mat1Rows+1, -mat1Cols+1);
		matCorrelation.Rotate(180);
	}
	
	return 0;	
}

//int correlation_2d_shift_accumulate(matrix& matSignal1, matrix& matSignal2, matrix& matCorrelation, bool bNormalize, bool bCircular)
int correlation_2d_shift_accumulate(matrixbase& matSignal1, matrixbase& matSignal2, matrixbase& matCorrelation, bool bNormalize, bool bCircular)
{
	int iRet;	

	int nRows, nCols;
	
	iRet = _check_2d_correlation_matrices(matSignal1, matSignal2, nRows, nCols);
	if(iRet)
		//return fft_error_handle(iRet);
		return iRet;

	///Sandy  2007-8-7 need to resize input matrices 
	//matrix mSignal1, mSignal2;
	//mSignal1 = matSignal1;
	//mSignal2 = matSignal2;	
	matrix mSignal1(nRows,nCols), mSignal2(nRows,nCols);
	mSignal1.SetSubMatrix(matSignal1, 0, 0);
	mSignal2.SetSubMatrix(matSignal2, 0, 0);
	//end 
	
	
	
	if(bNormalize) 	//normalize
	{
		_fft_normalize_response(mSignal1);
		_fft_normalize_response(mSignal2);
	}
		
	if(bCircular)
	{
		matCorrelation.SetSize(nRows, nCols);
		
		///Sandy  2008-1-21 MATRIXBASE_SHOULD_BE_USED
		matrix mCorrelation = matCorrelation;
		///Kevin 11/29/05 CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		//shift_correlation2D_circular(nRows, nCols, mSignal1, mSignal2, matCorrelation);
		//shift_2d_correlation_circular(nRows, nCols, mSignal1, mSignal2, matCorrelation);
		shift_2d_correlation_circular(nRows, nCols, mSignal1, mSignal2, mCorrelation);
		///End CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		
		matCorrelation = mCorrelation;
	}
	else
	{
		matCorrelation.SetSize(nRows*2-1, nCols*2-1);
		///Sandy  2008-1-21 MATRIXBASE_SHOULD_BE_USED
		matrix mCorrelation = matCorrelation;
		///Kevin 11/29/05 CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		//shift_correlation2D_linear(nRows, nCols, mSignal1, mSignal2, matCorrelation);
		//shift_2d_correlation_linear(nRows, nCols, mSignal1, mSignal2, matCorrelation);
		shift_2d_correlation_linear(nRows, nCols, mSignal1, mSignal2, mCorrelation);
		///End CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		
		matCorrelation = mCorrelation;
	}
	
	return 0;
}
///END ADD_2D_CORRELATION	

///Cheney 2006-7-17 MOVE_CODE_TO_XF_AND_LLVC
///Soapy 07/10/05 QA70-5060 ADD_COHERENCE_FUNCTION
/*int coherence(vector& vecSignal1, vector& vecSignal2, vector& vecCoherence, 
		int nWinSize, int iWindowType, int nOverlap, int nfft, double dAlpha, double dBeta)
{
	int iRet;

	///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//if(iRet = check_vector(vecSignal1))
	//if(iRet = check_vector(vecSignal2))
	if ( iRet = ocmath_count(NANUM, vecSignal1.GetSize(), vecSignal1) )
		//return fft_error_handle(iRet);
		return iRet;
	if ( iRet = ocmath_count(NANUM, vecSignal2.GetSize(), vecSignal2) )
		//return fft_error_handle(iRet);
		return iRet;
	///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	
	int nSize, nSize1;
	nSize = vecSignal1.GetSize();
	nSize1 = vecSignal2.GetSize();
	nSize = max(nWinSize, max(nSize, nSize1));
	vecSignal1.SetSize(nSize);
	vecSignal2.SetSize(nSize);
	
	
	// Check the arguments
	if(nWinSize == -1)
		nWinSize = nSize;
	else if(nWinSize < 1)
		//return fft_error_handle(INVALID_ARGUMENT);
		return CER_INVALID_ARGUMENT;
	
	if(nfft == -1)
		nfft = nWinSize;
	else if(nfft < 1)
		//return fft_error_handle(INVALID_ARGUMENT);
		return CER_INVALID_ARGUMENT;
	
	if(nOverlap > nWinSize)
		//return fft_error_handle(INVALID_ARGUMENT);
		return CER_INVALID_ARGUMENT;
		
	// Create a window function
	vector vecWindow(nWinSize);
	iRet = get_window_data(iWindowType, nWinSize, vecWindow, dAlpha, dBeta);
	if(iRet)
		//return fft_error_handle(INVALID_ARGUMENT);
		return CER_INVALID_ARGUMENT;

	// Number of segment
	int nStep = nWinSize - nOverlap;
	int nSegment = (nSize - nWinSize) / nStep + 1 ;
	
	// Predefine the vector of Pxx, Pyy, Pxy to sum
	vector<complex> vecPxx(nfft), vecPyy(nfft), vecPxy(nfft);
	
	vector<complex> vecFFT1, vecFFT2;   // For FFT
			
	for(int ii=0; ii < nSegment; ii++)
	{
		// Get a segment to caculate
		vecSignal1.GetSubVector(vecFFT1, ii*nStep, ii*nStep + nWinSize - 1);
		vecSignal2.GetSubVector(vecFFT2, ii*nStep, ii*nStep + nWinSize - 1);
		
		//Aplly the window to the signal	
		vecFFT1 *= vecWindow;
		vecFFT2 *= vecWindow; 
		
		vecFFT1.SetSize(nfft);
		vecFFT2.SetSize(nfft);

		// Perform FFT of both signals
		if( ( fft_complex(nfft, vecFFT1, FFT_FORWARD)) != 0) 
			//return fft_error_handle(UNKNOWN_ERROR);
			return CER_UNKNOWN_ERROR;
		if( ( fft_complex(nfft, vecFFT2, FFT_FORWARD)) != 0) 
			//return fft_error_handle(UNKNOWN_ERROR);
			return CER_UNKNOWN_ERROR;
	
		// Get conjugate of FFT
		vector<complex> vecFFTCONJ1(nfft),vecFFTCONJ2(nfft);
		vecFFTCONJ1 = vecFFT1;
		vecFFTCONJ2 = vecFFT2;
		vecFFTCONJ1.Conjugate();
		vecFFTCONJ2.Conjugate();
	
		// Sum the Pxx, Pyy, Pxy
		vecPxx += vecFFT1 * vecFFTCONJ1;
		vecPyy += vecFFT2 * vecFFTCONJ2;
		vecPxy += vecFFT2 * vecFFTCONJ1;	
	}
	
	// Select first half (positive frequency components) 
	int iselect;
	iselect = nfft / 2 + 1;
	vecPxx.SetSize(iselect);
	vecPyy.SetSize(iselect);
	vecPxy.SetSize(iselect);
	
	// Coherence function estimate
	vector vecReal(iselect);
	vecCoherence.SetSize(iselect);
	vecPxy.GetAmplitude(vecCoherence);
	vecCoherence *= vecCoherence;
	vecPxx *= vecPyy;
	vecPxx.GetAmplitude(vecReal);
	vecCoherence /= vecReal;
	
	return 0;
}*/
///END ADD_COHERENCE_FUNCTION
///end MOVE_CODE_TO_XF_AND_LLVC

/**
///Arvin 2006-5-19 ADD_MATRIX_ROTATE
int matrix_rotate(matrix<complex>& matData)
{
	int nRows = matData.GetNumRows();
	int nCols = matData.GetNumCols();
	
	if(0 == nRows||0 == nCols)
	{
		return INVALID_SIZE;
	}
	
	matrix<complex> matTemp;
	matTemp.SetSize(nRows, nCols);
	
	for(int ii = 0; ii < nRows; ii++)
		for(int jj = 0; jj < nCols; jj++)
		{
			matTemp[ii][jj] = matData[nRows-ii-1][nCols-jj-1];
		}
		
	matData = matTemp;
	return 0;
}
///END ADD_MATRIX_ROTATE

////Arvin 2006-5-19 ADD_FFT_TRANSFER_TO_MATLAB_FORMAT
int fft_transfer_to_matlab_format(matrix<complex>& matData)
{
	if(0 == matData.GetNumRows() || 0 == matData.GetNumCols())
		return INVALID_SIZE;
	
	matrix<complex> matTemp1, matTemp2, matTemp3;
	matData.CopyTo(matTemp1, 1, 1, -1, -1);
	matrix_rotate(matTemp1);
	matData.CopyTo(matTemp2, 0, 1, 0, -1);
	matrix_rotate(matTemp2);
	matData.CopyTo(matTemp3, 1, 0, -1, 0);
	matrix_rotate(matTemp3);
	matData.SetSubMatrix(matTemp1, 1, 1);
	matData.SetSubMatrix(matTemp2, 1, 0);
	matData.SetSubMatrix(matTemp3, 0, 1);
	return 0;
}
///END ADD_FFT_TRANSFER_TO_MATLAB_FORMAT
*/

///Cheney 2006-7-14 ADD_FOR_STFT
enum{
	REAL_TYPE = 0, COMPLEX_TYPE,
};
///end ADD_FOR_STFT

int stft_real(vector& vSignal, vector& vWin, double dSmplFreq, int nOverlap, int nFFTPts , matrix& matSTFT, double& dScaleX, double& dScaleY )
{

	///Cheney 2006-7-14 CENTRALIZE_CODE
	int nSize = vSignal.GetSize();
	int nWinSize, nShift, nStep;
	
	///Sandy 2006-11-21 ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX_AND_MOVE_THE_CODE_BACK
	/*///Sandy 2006-7-19 CHANGE_STFT_WITH_COMPLEX_OUTPUT_SPLIT_FUNCTION
		//if(_prepare_for_stft(vWin, matSTFT, nSize, nWinSize, nShift, nStep, nFFTPts, nOverlap, REAL_TYPE))
		if(_prepare_for_stft_real(vWin, matSTFT, nSize, nWinSize, nShift, nStep, nFFTPts, nOverlap))
			vSignal.SetSize(nWinSize);
	*/
	nWinSize = vWin.GetSize();
	nShift = nWinSize - nOverlap;	
	nStep = floor((double)(nSize - nWinSize) / (double)nShift) + 1 ;
				
	if(nStep < 1)
	{
		vSignal.SetSize(nWinSize);
		nStep = 1;
		nSize = nWinSize;
	}
	
	///Sandy 2006-11-21 Hard code for report the error, temp solution
	if(nStep>=STFT_MAX_OUTPUT_MATRIX_COLS ||(nFFTPts/2+1)>=STFT_MAX_OUTPUT_MATRIX_COLS)
		return CWAR_MATRUX_SIZE_NOT_SUPPORT;
	
	if(!matSTFT.SetSize(nStep, nFFTPts/2+1))
		return CER_FAILED_CREATE_MATRIX;
	
	matSTFT = 0;	
	///end ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX_AND_MOVE_THE_CODE_BACK
		
	
	//int nSize = vSignal.GetSize();
	//int nWinSize = vWin.GetSize();
	//int nPadding =	nFFTPts - nWinSize;
	//int nShift = nWinSize - nOverlap;
	
	//int nStep = floor((double)(nSize - nWinSize) / (double)nShift) + 1;
	//if(nStep < 1)
	//{
	//	vSignal.SetSize(nWinSize);		
	//	nStep = 1;
	//	nSize = nWinSize;
	//}
	
	//matSTFT.SetSize(nStep, nFFTPts/2+1);
	
	//matSTFT = 0;
	///end CENTRALIZE_CODE
	
	int nRet = fft_stft(vSignal, nSize, vWin, nWinSize, dSmplFreq, nShift, nFFTPts, nStep, matSTFT, &dScaleX, &dScaleY );
	if(nRet != 0 )
	{
		return nRet;
	}
	if(!matSTFT.Transpose())
		return XFERR_MATRIX_TRANSPOSE_FAILED;
	
	
	return 0;
}
///Sandy 2008-3-20 CLEAN_UP_bGetAmp
///Sandy 2006-7-19 CHANGE_STFT_WITH_COMPLEX_OUTPUT
//int stft_complex(vector<complex>& vSignal, vector& vWin, double dSmplFreq, int nOverlap, int nFFTPts ,matrix<complex>& matSTFT, matrix& matAmp, double& dScaleX, double& dScaleY, bool bGetAmp)
int stft_complex(vector<complex>& vSignal, vector& vWin, double dSmplIntv, int nOverlap, int nFFTPts ,matrix<complex>& matSTFT, matrix& matAmp, double& dScaleTime, double& dScaleFreq)
{
	
	///Cheney 2006-7-14 CENTRALIZE_CODE
	int nSize = vSignal.GetSize();
	int nWinSize, nShift, nStep;
	
	/*
	///Sandy 2006-7-19 CHANGE_STFT_WITH_COMPLEX_OUTPUT_SPLIT_FUNCTION
	//if(_prepare_for_stft(vWin, matSTFT, nSize, nWinSize, nShift, nStep, nFFTPts, nOverlap, REAL_TYPE))
	if(_prepare_for_stft_complex(vWin, matSTFT, nSize, nWinSize, nShift, nStep, nFFTPts, nOverlap))
		vSignal.SetSize(nWinSize);
	//int nSize = vSignal.GetSize();
	//int nWinSize = vWin.GetSize();
	//int nPadding =	nFFTPts - nWinSize;	
	//int nShift = nWinSize - nOverlap;
	
	//int nStep = floor((double)(nSize - nWinSize) / (double)(nWinSize - nOverlap)) + 1;
	//if(nStep < 1)
	//{
	//	vSignal.SetSize(nWinSize);		
	//	nStep = 1;
	//	nSize = nWinSize;
	//}
	
	//matSTFT.SetSize(nStep, nFFTPts);	
	//matSTFT = 0;
	///end CENTRALIZE_CODE
	*/
	nWinSize = vWin.GetSize();
	nShift = nWinSize - nOverlap;	
	nStep = floor((double)(nSize - nWinSize) / (double)(nWinSize - nOverlap)) + 1;

	bool bFlag = false;
	if(nStep < 1)
	{	
		nStep = 1;
		vSignal.SetSize(nWinSize);
		nSize = nWinSize;
	}
	
	///Sandy 2006-11-21 Hard code for report the error, temp solution
	if(nStep>=STFT_MAX_OUTPUT_MATRIX_COLS ||nFFTPts/2+1>=STFT_MAX_OUTPUT_MATRIX_COLS)
		return CWAR_MATRUX_SIZE_NOT_SUPPORT;
		
	if(!matSTFT.SetSize(nStep, nFFTPts))
		return CER_FAILED_CREATE_MATRIX;
	matSTFT = 0;

	
	int nRet;
	
	//if(bGetAmp)
	if(matAmp != NULL)
	{
		matAmp.SetSize(nStep, nFFTPts);
		nRet = fft_stft_complex(vSignal, nSize, vWin, nWinSize, dSmplIntv, nShift, nFFTPts, nStep, matSTFT, matAmp, &dScaleTime, &dScaleFreq )

	}
	else
	{
		nRet = fft_stft_complex(vSignal, nSize, vWin, nWinSize, dSmplIntv, nShift, nFFTPts, nStep, matSTFT, NULL, &dScaleTime, &dScaleFreq )
		
	}
	
	if(nRet!= 0)
		return nRet;
	
	//vector<complex> vComp;
	//matSTFT.GetAsVector(vComp);
	
	if(!matSTFT.Transpose())
		return XFERR_MATRIX_TRANSPOSE_FAILED;
	
	if(matAmp != NULL)
	{
		if(!matAmp.Transpose())
			return XFERR_MATRIX_TRANSPOSE_FAILED;	
	}

	return 0;
}


///Sandy 2006-11-21 ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX_AND_MOVE_THE_CODE_BACK
/*
///Sandy 2006-7-19 CHANGE_STFT_WITH_COMPLEX_OUTPUT_SPLIT_FUNCTION
/////Cheney 2007-7-14 ADD_FUNCTION_FOR_STFT
//static bool _prepare_for_stft(vector& vWin, matrix& matSTFT, int& nSize, 
			//int& nWinSize, int& nShift, int& nStep, int nFFTPts, int nOverlap, int nType = REAL_TYPE)
static bool _prepare_for_stft_real(vector& vWin, matrix& matSTFT, int& nSize, 
			int& nWinSize, int& nShift, int& nStep, int nFFTPts, int nOverlap)
{
	nWinSize = vWin.GetSize();
	nShift = nWinSize - nOverlap;	
	//nStep = (nType == REAL_TYPE) ? floor((double)(nSize - nWinSize) / (double)nShift) + 1 :
				//floor((double)(nSize - nWinSize) / (double)(nWinSize - nOverlap)) + 1;
	nStep = floor((double)(nSize - nWinSize) / (double)nShift) + 1 ;
				
	
	bool bFlag = false;
	if(nStep < 1)
	{	
		nStep = 1;
		nSize = nWinSize;
		bFlag = true;
	}
	
	
	///Sandy 2006-11-21 ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX
	//matSTFT.SetSize(nStep, nFFTPts/2+1);
	if(!matSTFT.SetSize(nStep, nFFTPts/2+1))
		return false;
	///end ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX
	
	matSTFT = 0;
	return bFlag;
}
///end ADD_FUNCTION_FOR_STFT

static bool _prepare_for_stft_complex(vector& vWin, matrix<complex>& matSTFT, int& nSize, 
			int& nWinSize, int& nShift, int& nStep, int nFFTPts, int nOverlap)
{
	nWinSize = vWin.GetSize();
	nShift = nWinSize - nOverlap;	
	nStep = floor((double)(nSize - nWinSize) / (double)(nWinSize - nOverlap)) + 1;

	bool bFlag = false;
	if(nStep < 1)
	{	
		nStep = 1;
		nSize = nWinSize;
		bFlag = true;
	}
	
	matSTFT.SetSize(nStep, nFFTPts);
	
	matSTFT = 0;
	return bFlag;
}
///end CHANGE_STFT_WITH_COMPLEX_OUTPUT_SPLIT_FUNCTION
*/ ///end ADD_WARNING_FOR_TOO_LARGE_RESULT_MATRIX_AND_MOVE_THE_CODE_BACK
/*
/// AW Soapy 05/28/04 QA75-6457 v7.5847 CLEANUP_STFT 
///Forest 08/16/04 STFT_INTERFACE_UPDATE, vWindow is redundant argument
//int stft(vectorbase& vbSig, vector& vWindow, const STFTOptions& stftOp, matrix& matRet, fpoint& fpScale)
int stft(vectorbase& vbSig, const STFTOptions& stftOp, matrix& matRet, fpoint& fpScale)
{
	if( stftOp.FFTLength < 1 || stftOp.OverlapPt < 0 || stftOp.WinLength < 0 )
		return	fft_error_handle(INVALID_ARGUMENT);
	if( stftOp.OverlapPt >= stftOp.WinLength || stftOp.FFTLength < stftOp.WinLength)
		return	fft_error_handle(INVALID_ARGUMENT);
	
	vector vWindow(stftOp.WinLength);
///End STFT_INTERFACE_UPDATE
	int nRet;
	if ( 0 > ( nRet = get_stft_window_data( stftOp.WindowMethod, stftOp.WinLength, stftOp.Alpha, stftOp.Belta, vWindow ) ))
	{	
		ASSERT(FALSE);
		return nRet;
	}
	
	if ( 0 != ( nRet = stft(vbSig,vWindow,matRet,stftOp.OverlapPt,stftOp.FFTLength, stftOp.SamplingFq, fpScale) ) )
	{	
		ASSERT(FALSE);
		return nRet;
	}
	//2006-3-7 sandy only_comment
	//double d1 = matSTFT[0][2];
	return 0;
}

//
//|============================|
       //Win Size    |============================|
       	          //Overlap            |============================|
//
//|==================================|
       //FFT Size           |=========|
                            //Padding
//

int stft(vectorbase& vbSig, vector& vecWin, matrix& matSTFT, int nOverlap, int nFFTPts, double dFs, fpoint& fpScale )
{
	///Raine/Leo 2006-3-10 STFT_RESULT_IS_DIFFERENT_FROM_AUTOSIGNAL
	//	int nRet;
	//	///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//	//if(nRet = check_vector(vecWin))	return nRet;
	//	if ( nRet = ocmath_count(NANUM, vecWin.GetSize(), vecWin) ) 
	//		return INVALID_DATA;	///Mouqx 2005-9-22
	//	///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//	
	//	int nWinSize = vecWin.GetSize();
	//	int nShiftPts = nWinSize - nOverlap;
	//	int nPadding =	nFFTPts - nWinSize;
	//	
	//	//if((nShiftPts < 1) || (nPadding < 0) || (dFs <= 0))	///Forest 08/16/04 RETURN_VALUE_UNIFICATION, Matlab supports negative dFs
	//	if((nShiftPts < 1) || (nPadding < 0) || (dFs == 0))
	//		return fft_error_handle(INVALID_ARGUMENT);
	//
	//	///Austin/Leo 09/20/05	REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//	//nRet = check_vector(vbSig);
	//	if (nRet = ocmath_count(NANUM, vbSig.GetSize(), (vector)vbSig))
	//		return INVALID_DATA;	///Mouqx 2005-9-22
	//	///End REPLACE_CHECK_VECTOR_BY_OCMATH_COUNT
	//	
	//	int iSigSize = vbSig.GetSize();
	//	
	//	// padding zero to the beginning of the signal to make the maximum point of the window at 
	//	// the first point of the signal when apply the window to the signal.
	//	int iPosWinMax;
	//	double min, max;
	//	vecWin.GetMinMax(min, max, NULL, &iPosWinMax);
	//	/// AW 06/29/04 MORE_ON_STFT_SIZE
	//	// also need add the end padding part
	//	
	//	vector vecPadSig(iPosWinMax);
	//
	//	if(iPosWinMax>0)
	//	{
	//		vecPadSig = 0;
	//		vecPadSig.Append(vbSig);
	//	}
	//	else
	//		vecPadSig=vbSig;
	//	
	//	iSigSize += iPosWinMax;
	//	
	//	vector vecPadSig(iPosWinMax), vecPadEndSig(nWinSize - iPosWinMax);
	//
	//	if(iPosWinMax>0)
	//	{
	//		vecPadSig = 0;
	//		vecPadSig.Append(vbSig);
	//		vecPadEndSig = 0;
	//		vecPadSig.Append(vecPadEndSig);
	//	}
	//	else
	//	{
	//		vecPadSig=vbSig;
	//		vecPadEndSig = 0;
	//		vecPadSig.Append(vecPadEndSig);
	//	}
	//		
	//	iSigSize += nWinSize;
	//	/// END MORE_ON_STFT_SIZE
	//	
	//	
	//	//Set the property of the output Matrix
	//	int iStep = floor((double)(iSigSize-nOverlap) / (double)nShiftPts);
	//	if(iStep<1)
	//	{
	//		vecPadSig.SetSize(nWinSize);
	//		iStep = 1;
	//		iSigSize = nWinSize;
	//	}
	//	
	//	matSTFT.SetSize(iStep, nFFTPts/2+1);
	//	
	//	matSTFT = 0;	// clear all data
	//	
	//	//Perform STFT
	//	if(nRet = fft_stft(vecPadSig,iSigSize,vecWin,nWinSize,nShiftPts,nPadding,iStep,matSTFT))
	//		return nRet;
	//	
	//	matSTFT.Transpose();

	//	fpScale.y = (nFFTPts/2)*dFs/nFFTPts;
	//	fpScale.x = 1/dFs*nShiftPts* (iStep-1);

	int nRet;
	if ( nRet = ocmath_count(NANUM, vecWin.GetSize(), vecWin) ) 
	{
		return INVALID_DATA;
	}

	int nWinSize = vecWin.GetSize();
	int nPadding =	nFFTPts - nWinSize;
	
	if(((nWinSize - nOverlap) < 1) || (nPadding < 0) || (dFs == 0))
	{
		return fft_error_handle(INVALID_ARGUMENT);
	}

	if (nRet = ocmath_count(NANUM, vbSig.GetSize(), (vector)vbSig))
	{
		return INVALID_DATA;
	}
	
	int nDataSize = vbSig.GetSize();
	vector vecSig;
	vecSig.Append(vbSig);	

	int nSegments = floor((double)(nDataSize - nWinSize) / (double)(nWinSize - nOverlap)) + 1;
	if(nSegments<1)
	{
		vecSig.SetSize(nWinSize);		
		nSegments = 1;
		nDataSize = nWinSize;
	}
	
	matSTFT.SetSize(nSegments, nFFTPts/2+1);
	
	matSTFT = 0;

	if(nRet = fft_stft(vecSig, nDataSize, vecWin, nWinSize, (nWinSize - nOverlap), nPadding, nSegments, matSTFT))
	{
		return nRet;
	}
	matSTFT.Transpose();	
	
	fpScale.y = (nFFTPts/2)*dFs/nFFTPts;
	fpScale.x = 1/dFs*(nWinSize - nOverlap)* (nSegments-1);
	
	return 0;
}
///END CLEANUP_STFT
///END STFT_RESULT_IS_DIFFERENT_FROM_AUTOSIGNAL

///Raine 2006-04-06 ADD_STFT_COMPLEX
int stft_complex(vector& vReal, vector& vImag, matrix& matSTFT, const STFTOptions& stftOp, fpoint& fpScale)
{
	if( stftOp.FFTLength < 1 || stftOp.OverlapPt < 0 || stftOp.WinLength < 0 )
		return	fft_error_handle(INVALID_ARGUMENT);
	if( stftOp.OverlapPt >= stftOp.WinLength || stftOp.FFTLength < stftOp.WinLength)
		return	fft_error_handle(INVALID_ARGUMENT);
	
	vector vecWin(stftOp.WinLength);
	int nRet;	
	if ( 0 > ( nRet = get_stft_window_data( stftOp.WindowMethod, stftOp.WinLength, stftOp.Alpha, stftOp.Belta, vecWin ) ))
	{	
		ASSERT(FALSE);
		return nRet;
	}	
	if ( nRet = ocmath_count(NANUM, vecWin.GetSize(), vecWin) ) 
	{
		return INVALID_DATA;
	}		
	if (nRet = ocmath_count(NANUM, vReal.GetSize(), vReal))
	{
		return INVALID_DATA;
	}
	if (nRet = ocmath_count(NANUM, vImag.GetSize(), vImag))
	{
		return INVALID_DATA;
	}
	
	int nDataSize = vReal.GetSize();
	int nPadding =	stftOp.FFTLength - stftOp.WinLength;	
	int nSegments = floor((double)(nDataSize - stftOp.WinLength) / (double)(stftOp.WinLength - stftOp.OverlapPt)) + 1;
	if(nSegments<1)
	{
		vReal.SetSize(stftOp.WinLength);		
		vImag.SetSize(stftOp.WinLength);
		nSegments = 1;
		nDataSize = stftOp.WinLength;
	}
	
	matSTFT.SetSize(nSegments, stftOp.FFTLength);	
	matSTFT = 0;
	
	///Raine 2006-04-13 CHANGE_STFT_PARAMETER
	int nSigSize = stftOp.FFTLength * nSegments;
	vector<complex> vSig(nSigSize);	
	if(nRet = fft_apply_window_complex(nDataSize, vReal, vImag, stftOp.WinLength, vecWin, (stftOp.WinLength - stftOp.OverlapPt), nSigSize, vSig, nPadding, nSegments))
	{
		return nRet;
	}	
	if(nRet = fft_stft_complex(vSig, nDataSize, vecWin, stftOp.WinLength, (stftOp.WinLength - stftOp.OverlapPt), nPadding, nSegments, matSTFT))
	{
		return nRet;
	}
	///END CHANGE_STFT_PARAMETER
	matSTFT.Transpose();
	
	//2006-3-7 sandy only_comment
	//double d1 = matSTFT[0][2];
	
	fpScale.y = (stftOp.FFTLength/2)*stftOp.SamplingFq/stftOp.FFTLength;
	fpScale.x = 1/stftOp.SamplingFq*(stftOp.WinLength - stftOp.OverlapPt)* (nSegments-1);
	
	return 0;
}
///END ADD_STFT_COMPLEX
*/

///Soapy 07/10/05 QA70-7865 ADD_2D_CORRELATION
static void _fft_normalize_response(matrix& matResponse)
{
	matrix msquare;
	msquare = matResponse;
	msquare.DotMultiply(matResponse);
	
	double dSumResp;
	dSumResp = msquare.GetMean() * matResponse.GetNumCols() * matResponse.GetNumRows();
	
	matResponse /= sqrt(dSumResp);
}
///END ADD_2D_CORRELATION


///Sandy  2006-7-17 COMMENT_SEEM_TO_NO_USED_CODE
///// Sandy 7/12/06 MOVE_FFT_RELATIVE_FUNCIONS_BACK_TO_FFT_UTILS
//
//// This function will do the 2d fft with the data form required by nag funtion.
////static int fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImag, matrix<complex>& matFFT2Result);
//static int _fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImag, matrix<complex>& matFFT2Result)
//{
	/////Iris 2/13/04 USING_FFT_REPLACE_NAG
	///*
	//vector vecTrigRows;
	//vecTrigRows.SetSize(2 * nNumRows);
	////cal nage function to prepare the trig value
	//int nRet;
	//if( nRet = nag_fft_init_trig(nNumRows, vecTrigRows))
		//return nRet;
	//
	//vector vecTrigCols;
	//vecTrigCols.SetSize(2 * nNumCols);
	//if( nRet = nag_fft_init_trig(nNumCols, vecTrigCols))
		//return nRet;
	//
	////call the nag function to do the fft2
	//if(nRet = nag_fft_2d_complex(nNumRows, nNumCols, matReal, matImag, vecTrigRows, vecTrigCols))
		//return nRet;
	//*/
	//int nRet;
	//if(nRet = fft_fft_2d_complex(nNumRows,nNumCols,matReal,matImag))
		//return nRet;	
	/////End USING_FFT_REPLACE_NAG
	//
	////set the data back
	//if( nRet = matFFT2Result.MakeComplex(matReal, matImag))
		//return nRet;
	//
//}
//
///**
//
	//int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	//int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
			//These two Functions will do the 2 dimension Fast Fourier Transform.
	//Remarks:
			//Right now, we only support the data matrix in double and complex form. If user have the data
			//matrix in integer form, user should first transform it to double. It is a member function of matrix.
			//Also, if user has a vector data, it should firstly be transformed into a matrix and then do the 
			//FFT.
	//Example:
	//Parameters:
			//matSource: The data source matrix
			//matFFT2Result: The matrix containing the 2 dimension FFT results
			//nRowProcess:
			//nColProcess: The number of rows or columns user provided. If they are -1, there is no padding
	 					//or truncating. If they are more than the sizes of the matrix, the function will pad zero to meet the size
						//If less than, it will truncate the data matrix.
	//Return:
			//If succeed, it will return 0;	Otherwise it returns error indications.
			//Error:
			//-1: Provided column number or row number. Must be -1 or positive integer
			//-2: The data source matrix is empty 
			//-3: When calling member function, it fails.
			//-4: When calling GetReal or GetImage member function of matrix class, it failed.
			//Nag Error:
			 //11: The column or the row size of matrix is less then 1
 			 //73: System failed to allocate memory
//
//*/
//int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
//{
	//matrix matReal;
	//matrix matTemp(matSource);
	//int nRet;
	//if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess))
		//return nRet;
			//
	//int nNumCols = matReal.GetNumCols();
	//int nNumRows = matReal.GetNumRows();
	//matrix matImag;
	//matImag.SetSize(nNumRows, nNumCols);
	//matImag = 0;	
	//return _fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImag, matFFT2Result);
//}
	//
//int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess )
//{
	//int nRet;
	//matrix matReal, matImag;
	//if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImag, nRowProcess, nColProcess))
		//return nRet;
	//
	//int nNumCols = matImag.GetNumCols();
	//int nNumRows = matImag.GetNumRows();
	//return _fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImag, matFFT2Result);	
//}
//
///*
	//int IFFT2(matrix<complex>& matSource, matrix<complex>& matIFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	//Remarks:
			//Right now, we only support the data matrix in complex form. If user have the data
			//matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			//Matrix class has a member function to make the complex matrix. Also, if user has a vector data,
			//it should firstly be transformed into a matrix and then one can do the inverse FFT.
	//Example:
	//Parameters:
			//matSource: The data source matrix
			//matIFFT2Result: The matrix containing the 2 dimension inverse FFT results
			//nRowProcess:
			//nColProcess:The number of rows or columns user provided. if they are -1, there are no padding
	 					//or truncating. if they are larger than the size of the matrix, the function will pad zero to meet
						//the size. If less than, it will truncate the data matrix.
	//Return:
			//If succeed, it will return 0;	Otherwise it returns error indications.
			//Error:
			//-1: Provided column number or row number. Must be -1 or positive integer
			//-2: The data source matrix is empty 
			//-3: When calling member function, it fails.
			//-4: When calling GetReal or GetImage member function of matrix class, it failed.
			//Nag Error:
			 //11: The column or the row size of matrix is less then 1
 			 //73: System failed to allocate memory
//
  //
//*/
//int IFFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
//{
//
	//int nRet;
	//matrix<complex> matSourceCopy(matSource);
	//if(nRet = matSourceCopy.Conjugate())
		//return nRet;
	//
	//if(nRet = FFT2(matSourceCopy, matFFT2Result,nRowProcess, nColProcess ))
		//return nRet;
	//
	//return matFFT2Result.Conjugate();
//}		
///*
	//int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	//int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
//
				//These two Functions will do the 1 dimension Fast Fourier Transform.
  	//Remarks:
			//Right now, we only support the data matrix in double and complex form. If user has the data
			//matrix in integer form, user should firstly transform it to double. It is a member function of matrix.
			//Also, if user has a vector data, he should firstly transform data into a matrix and then one can do the 
			//FFT. Unlike the matlab, this function will do the FFT on each row (Matlab do FFT on each column).
	//Example:
	//Parameters:
			//matSource: The data source matrix
			//matFFTResult: The matrix containing the 1 dimension FFT results
			//nColProcess: The number of columns user provided. If it is -1, function will not do any padding
	 					//or truncating. But if it is larger than the size of the  matrix, it will pad zero to meet the size
						//if less than, it will truncate the data matrix.
	//Return:
			//If succeed, it will return 0;	Otherwise it returns error indications.
			//Error:
			//-1: Provided column number or row number. Must be -1 or positive integer
			//-2: The data source matrix is empty 
			//-3: When calling member function, it fails.
			//-4: When calling GetReal or GetImage member function of matrix class, it failed.
			//Nag Error:
			 //11: The column or the row size of matrix is less then 1
 			 //73: System failed to allocate memory
//
//*/
//
//int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess)
//{
	//
	//matrix matTemp;
	//int nRet;
	//if(nRet = set_matrix_with_padding_truncting(matSource, matTemp, -1, nColProcess))
		//return nRet;
	//int nNumCols = matTemp.GetNumCols();
	//int nNumRows = matTemp.GetNumRows();
	//
	/////Iris 2/13/04 USING_FFT_REPLACE_NAG
	///*
	//vector vecTrig;
	//vecTrig.SetSize(2 * nNumCols); 
	//if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		//return nRet;
	//if(nRet = nag_fft_multiple_real(nNumRows, nNumCols, matTemp, vecTrig))
		//return nRet;
	//matrix matReal, matImag;
	//matReal.SetSize(nNumRows, nNumCols);
	//matImag.SetSize(nNumRows, nNumCols);
	//
	//if(nRet = nag_multiple_hermitian_to_complex(nNumRows, nNumCols, matTemp, matReal, matImag))
		//return nRet;
	//*/
	//matrix matReal, matImag;
	//matReal.SetSize(nNumRows, nNumCols);
	//matImag.SetSize(nNumRows, nNumCols);
	//matReal = matTemp;
	//if(nRet = fft_fft_multiple_complex(nNumRows, nNumCols, matReal, matImag, FFT_FORWARD))
		//return nRet;
	/////End USING_FFT_REPLACE_NAG
	//
	//return matFFTResult.MakeComplex(matReal, matImag);
//}
		//
//int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
//{
	//int nRet;
	//matrix matReal, matImag;
	//if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImag, -1, nColProcess))
		//return nRet;
	//int nNumCols = matImag.GetNumCols();
	//int nNumRows = matImag.GetNumRows();
	/////Iris 2/13/04 USING_FFT_REPLACE_NAG
	///*
	//vector vecTrig;
	//vecTrig.SetSize(2 * nNumCols); 
	//if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		//return nRet;
	//if(nRet = nag_fft_multiple_complex(nNumRows, nNumCols, matReal, matImag, vecTrig ))
		//return nRet;
	//*/
	/////End USING_FFT_REPLACE_NAG
	//if(nRet = fft_fft_multiple_complex(nNumRows, nNumCols, matReal,matImag,FFT_FORWARD))
		//return nRet;
	//
	//return matFFTResult.MakeComplex(matReal, matImag);
//}
//
//// ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
//// The following function performs 1-D FFT on a complex vector.
//// It uses the matrix-based 1-D FFT function.
//
//int FFT(vector<complex>& vecSource, vector<complex>& vecFFTResult)
//{
	//// Get size of source vector
	//int iSize = vecSource.GetSize();
	//if(iSize < 1) 
		//return CER_INVALID_SIZE;
	//
	//// Create a complex matrix of size 1xn and copy source vector into matrix
	//matrix<complex> matSource;
	//matSource.SetSize(1,iSize);
	//matSource.SetByVector(vecSource);
	//
	//// Create complex matrix to hold result of FFT operation
	//matrix<complex> matFFTResult;
	//
	//// Call the matrix-based FFT function to perform FFT on this one row
	//int iRet = FFT(matSource, matFFTResult);
	//if(iRet) 
		//return iRet;
//
	//// Extract FFT result from result matrix, into a complex matrix, and return
	//matFFTResult.GetAsVector(vecFFTResult);
	//
	//// return 0 on success
	//return 0;
//}
//
///*	  
	//int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	//
			//This Function will do the 1 dimension inverse Fast Fourier Transform.
	//Remarks:
			//Right now, we only support the data matrix in complex form. If user has the data
			//matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			//Matrix class has a member function to make the complex matrix. also, if user has a vector data,
			//it should be firstly transformed into a matrix and then one can do the inverse FFT.
//
	//Example:
	//Parameters:
			//matSource:The data source matrix
			//matIFFT2Result: The matrix containing the 1 dimension inverse FFT results
			//nColProcess:The number of columns user provided. if it is -1, there is no padding
	 					//or truncting. if it is more than the size of the  matrix, it will pad zero to meet the size
						//if less than, it will truncate the data matrix.
	//Return:
			//If succeed, it will return 0;	Otherwise it returns error indications.
			//Error:
			//-1: Provided column number or row number must be -1 or positve integer
			//-2: The data source matrix is empty 
			//-3: When calling member function, it fails.
			//-4: When calling GetReal or GetImage member function of matrix class, it failed.
			//Nag Error:
			 //11: The coloumn or the row size of matrix is less then 1;
 			 //73: System failed to allocate memory
//
  //
//*/
//
//
		//
//int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
//{
	//int nRet;
	//matrix<complex> matSourceCopy(matSource);
	//if(nRet = matSourceCopy.Conjugate())
		//return nRet;
//
	//if(nRet = FFT(matSourceCopy, matFFTResult, nColProcess))
		//return nRet;
//
	//return matFFTResult.Conjugate();
//}	
//
//
//// ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
//// The following function performs 1-D inverse FFT on a complex vector.
//// It uses the matrix-based 1-D inverse FFT function.
//
//int IFFT(vector<complex>& vecSource, vector<complex>& vecFFTResult)
//{
	//// Get size of source vector
	//int iSize = vecSource.GetSize();
	//if(iSize < 1) 
		//return CER_INVALID_SIZE;
	//
	//// Create a complex matrix of size 1xn and copy source vector into matrix
	//matrix<complex> matSource;
	//matSource.SetSize(1,iSize);
	//matSource.SetByVector(vecSource);
	//
	//// Create complex matrix to hold result of IFFT operation
	//matrix<complex> matFFTResult;
	//
	//// Call the matrix-based IFFT function to perform IFFT on this one row
	//int iRet = IFFT(matSource, matFFTResult);
	//if(iRet) 
		//return iRet;
//
	//// Extract IFFT result from result matrix, into a complex matrix, and return
	//matFFTResult.GetAsVector(vecFFTResult);
	//
	//// return 0 on success
	//return 0;
//}
//
//
/////end MOVE_FFT_RELATIVE_FUNCIONS_BACK_TO_FFT_UTILS
//
///end  2006-7-17 COMMENT_SEEM_TO_NO_USED_CODE