fftengine.cpp

选自<gpu gemes 2>,用gpu实现快速傅立叶变换

		cgGLSetTextureParameter(yBaseTexR2_1, Texs[2]);
		CheckCgError();
		cgGLSetTextureParameter(yBaseTexI2_1, Texs[3]);
		CheckCgError();
	}
	else {
		cgGLSetTextureParameter(yBaseTexR1_2, Texs[0]);
		CheckCgError();
		cgGLSetTextureParameter(yBaseTexI1_2, Texs[1]);
		CheckCgError();
		cgGLSetTextureParameter(yBaseTexR2_2, Texs[2]);
		CheckCgError();
		cgGLSetTextureParameter(yBaseTexI2_2, Texs[3]);
		CheckCgError();
	}
	EnableFFTFragmentProgramY();

	if (CurrentButterflyStage < yCutOff) {
		// Debugging using immediate mode
		// DrawQuadForFFT();
		// Using display lists
		DrawQuadForFFT_List();
	}
	else {
		// Debugging using immediate mode
		// DrawQuadTilesForFFT_Y();
		// Using display lists
		DrawQuadTilesForFFT_Y_List();
	}
	DisableFFTFragmentProgramY();
		
	Buffers->Pong(FirstTime_p);
	Buffers->Swap();
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadTilesForFFT_X
*
* DESCRIPTION:
*	Draws quads for method 2 in doing FFT in x. This function used only
*	for debugging. It has been replaced by DrawQuadTilesForFFT_X_List
*	function, which uses display lists and draw arrays.
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DrawQuadTilesForFFT_X(void)
{
	glActiveTextureARB(GL_TEXTURE0_ARB);
	int i = CurrentButterflyStage;	
	int j, k, kk;
	float tllx, turx, llx, urx;
	int nBlocks  = powf(2.0, (float)(nButterfliesX - 1 - i));
	int Inc      = powf(2.0, (float)(i+1));
	int IncB     = powf(2.0, (float)(i)) - 1.0;
	int nHInputs = powf(2.0, (float)(i));
	float BlockSize = (float)(Width)/(float)nBlocks/2.0;
	float Index1Start, Index1End, Index2Start, Index2End;
	float Index1Start_, Index1End_, Index2Start_, Index2End_;
	float AngStart, AngEnd, AngScale;
	float v1, v2, Grad;
	float AngSign = (ForwardFFT_p) ? 1.0 : -1.0;

	if (i == 0) {
		AngStart = 0.0;
		AngEnd   = 0.0;
	}
	else {
		v1 = 0.0;
		v2 = nBlocks*IncB;
		Grad = (v2 - v1)/(float)(BlockSize-1);
		AngScale = -2.0*M_PI/(float)Width;
		AngStart = AngScale*(v1 - 0.5*Grad);
		AngEnd   = AngScale*(v2 + 0.5*Grad);
	}

	glEnable(GL_TEXTURE_RECTANGLE_NV);
	Index1Start = 0;
	Index2Start = powf(2.0, (float)(i));

	for (k = 0; k < nBlocks; k++) {
		Index1End   = Index1Start + IncB;
		Index2End   = Index2Start + IncB;

		if (i == 0) {
			Index1Start_ = BitReverse((int)Index1Start, Width);
			Index2Start_ = BitReverse((int)Index2Start, Width);
			Index1End_   = BitReverse((int)Index1End, Width);
			Index2End_   = BitReverse((int)Index2End, Width);
		}
		else {
			Index1Start_ = Index1Start;
			Index2Start_ = Index2Start;
			Index1End_   = Index1End;
			Index2End_   = Index2End;
		}

		for (j = 0; j < 2; j++) {
			kk = 2*k + j;

			glBegin(GL_QUADS);

			tllx = -0.5 + kk*BlockSize;
			turx = tllx + BlockSize;

			llx = kk*BlockSize;
			urx = llx + BlockSize;

			glTexCoord4f(tllx, -0.5, Index1Start_-0.5, Index2Start_-0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngStart+j*M_PI));
			glVertex2f(llx, 0.0);

			glTexCoord4f(tllx, Height-1+0.5, Index1Start_-0.5, Index2Start_-0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngStart+j*M_PI));
			glVertex2f(llx, Height);

			glTexCoord4f(turx, Height-1+0.5,  Index1End_+0.5, Index2End_+0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngEnd+j*M_PI));
			glVertex2f(urx, Height);

			glTexCoord4f(turx, -0.5, Index1End_+0.5, Index2End_+0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngEnd+j*M_PI));
			glVertex2f(urx, 0.0);

			glEnd();
		}	
		Index1Start += Inc;
		Index2Start += Inc;
	}
	glDisable(GL_TEXTURE_RECTANGLE_NV);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadTilesForFFT_X_List
*
* DESCRIPTION:
*	Draws quads for method 2 in doing FFT in x using display lists and
*	draw arrays.
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DrawQuadTilesForFFT_X_List(void)
{
	glActiveTextureARB(GL_TEXTURE0_ARB);
	glEnable(GL_TEXTURE_RECTANGLE_NV);
	glCallList(ListX+CurrentButterflyStage);
	glDisable(GL_TEXTURE_RECTANGLE_NV);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadTilesForFFT_Y
*
* DESCRIPTION:
*	Draws quads for method 2 in doing FFT in y. This function used only
*	for debugging. It has been replaced by DrawQuadTilesForFFT_Y_List
*	function, which uses display lists and draw arrays.
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DrawQuadTilesForFFT_Y(void)
{
	glActiveTextureARB(GL_TEXTURE0_ARB);
	int i = CurrentButterflyStage;	
	int j, k, kk;
	float tlly, tury, lly, ury;
	int nBlocks  = powf(2.0, (float)(nButterfliesY - 1 - i));
	int Inc      = powf(2.0, (float)(i+1));
	int IncB     = powf(2.0, (float)(i)) - 1.0;
	int nHInputs = powf(2.0, (float)(i));
	float BlockSize = (float)(Height)/(float)nBlocks/2.0;
	float Index1Start, Index1End, Index2Start, Index2End;
	float Index1Start_, Index1End_, Index2Start_, Index2End_;	
	float AngStart, AngEnd, AngScale;
	float v1, v2, Grad;
	float AngSign = (ForwardFFT_p) ? 1.0 : -1.0;

	if (i == 0) {
		AngStart = 0.0;
		AngEnd   = 0.0;
	}
	else {
		v1 = 0.0;
		v2 = nBlocks*IncB;
		Grad = (v2 - v1)/(float)(BlockSize-1);
		AngScale = -2.0*M_PI/(float)Height;
		AngStart = AngScale*(v1 - 0.5*Grad);
		AngEnd   = AngScale*(v2 + 0.5*Grad);
	}

	glEnable(GL_TEXTURE_RECTANGLE_NV);
	Index1Start = 0;
	Index2Start = powf(2.0, (float)(i));

	for (k = 0; k < nBlocks; k++) {
		Index1End   = Index1Start + IncB;
		Index2End   = Index2Start + IncB;

		if (i == 0) {
			Index1Start_ = BitReverse((int)Index1Start, Height);
			Index2Start_ = BitReverse((int)Index2Start, Height);
			Index1End_   = BitReverse((int)Index1End, Height);
			Index2End_   = BitReverse((int)Index2End, Height);
		}
		else {
			Index1Start_ = Index1Start;
			Index2Start_ = Index2Start;
			Index1End_   = Index1End;
			Index2End_   = Index2End;
		}

		for (j = 0; j < 2; j++) {
			kk = 2*k + j;

			glBegin(GL_QUADS);

			tlly = -0.5 + kk*BlockSize;
			tury = tlly + BlockSize;

			lly = kk*BlockSize;
			ury = lly + BlockSize;

			glTexCoord4f(-0.5, tlly, Index1Start_-0.5, Index2Start_-0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngStart+j*M_PI));
			glVertex2f(0.0, lly);

			glTexCoord4f(Width-1+0.5, tlly, Index1Start_-0.5, Index2Start_-0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngStart+j*M_PI));
			glVertex2f(Width, lly);

			glTexCoord4f(Width-1+0.5, tury,  Index1End_+0.5, Index2End_+0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngEnd+j*M_PI));
			glVertex2f(Width, ury);

			glTexCoord4f(-0.5, tury, Index1End_+0.5, Index2End_+0.5);
			glMultiTexCoord1f(GL_TEXTURE1_ARB, AngSign*(AngEnd+j*M_PI));
			glVertex2f(0.0, ury);

			glEnd();

		}	
		Index1Start += Inc;
		Index2Start += Inc;
	}
	glDisable(GL_TEXTURE_RECTANGLE_NV);

}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadTilesForFFT_Y_List
*
* DESCRIPTION:
*	Draws quads for method 2 in doing FFT in y using display lists and
*	draw arrays.
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DrawQuadTilesForFFT_Y_List(void)
{
	glActiveTextureARB(GL_TEXTURE0_ARB);
	glEnable(GL_TEXTURE_RECTANGLE_NV);
	glCallList(ListY+CurrentButterflyStage);
	glDisable(GL_TEXTURE_RECTANGLE_NV);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::ComputeWeight
*
* DESCRIPTION:
*	Computes the weights.
*
* FORMAL PARAMETERS:
*	N:	Number of samples
*   k:	Current sample
*
* RETURNS:
*	Wr: real part of the weight
*	Wi: imaginary part of the weight
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::ComputeWeight(int N, int k, float &Wr, float &Wi)
{
	Wr =  cosl(2.0*M_PI*k/(float)N);
	Wi = -sinl(2.0*M_PI*k/(float)N);

	Wi = (ForwardFFT_p == true) ? Wi : -Wi;
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::BitReverse
*
* DESCRIPTION:
*	Reverses bits in index
*
* FORMAL PARAMETERS:
*	i:	input index
*   N:	Number of bits in the index
*
* RETURNS:
*	the bit-reversed index
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
int FFT::BitReverse(int i, int N)
{
	int j = i;

	int M   = N;
	int Sum = 0;
	int W   = 1;
	M = M / 2;
	while (M != 0) {
		j = (i & M) > M-1;
		Sum += j*W;
		W *= 2;
		M = M/2;
	}
	return(Sum);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::CreateButterflyLookups
*
* DESCRIPTION:
*	Creates scrambling indices and weights for each butterfly stage
*
* FORMAL PARAMETERS:
*	NButterflies:	number of butterfly stages
*   N:				number of samples
*
* RETURNS:
*	butterflylookupI:	an array containing scrambling lookup table
*	butterflylookupWR:	real part of the weights for each stage
*	butterflylookupWI:	imaginary part of the weights for each stage
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::CreateButterflyLookups(float *butterflylookupI,
								 float *butterflylookupWR, 
								 float *butterflylookupWI, 
								 int NButterflies, int N)
{
	float *ptr0 = butterflylookupI;
	float *ptr1 = butterflylookupWR;
	float *ptr2 = butterflylookupWI;

	int i, j, k, i1, i2, j1, j2;
	int nBlocks, nHInputs;
	float wr, wi;
	float *qtr0, *qtr1, *qtr2;
	float scale = 1.0/((float)(N-1));

	for (i = 0; i < NButterflies; i++) {
		nBlocks  = powf(2.0, (float)(NButterflies - 1 - i));
		nHInputs = powf(2.0, (float)(i));
		qtr0 = ptr0;
		qtr1 = ptr1;
		qtr2 = ptr2;
		for (j = 0; j < nBlocks; j++) {

			for (k = 0; k < nHInputs; k++) {

				if (i == 0) {
					i1 = j*nHInputs*2 + k;
					i2 = j*nHInputs*2 + nHInputs + k;
					j1 = BitReverse(i1, N);
					j2 = BitReverse(i2, N);
				}
				else {
					i1 = j*nHInputs*2 + k;
					i2 = j*nHInputs*2 + nHInputs + k;
					j1 = i1;
					j2 = i2;
				}

				ComputeWeight(N, k*nBlocks, wr, wi);

				*(qtr0 + 2*i1)   = j1;
				*(qtr0 + 2*i1+1) = j2;
				*(qtr1 + i1) = wr;
				*(qtr2 + i1) = wi;

				*(qtr0 + 2*i2)   = j1;
				*(qtr0 + 2*i2+1) = j2;
				*(qtr1 + i2) = -wr;
				*(qtr2 + i2) = -wi;

			}
		}
		ptr0 += 2*N;
		ptr1 += N;
		ptr2 += N;
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::ComputeVerticesTexCoordsX
*
* DESCRIPTION:
*	Computes scrambling indices and weights for each sub-quad for method 2
*
* FORMAL PARAMETERS:
*	Stage:	buttefly stage
*
* RETURNS:
*	Vs:		vertices
*	Ts0:	texture coordinates for texture unit 0
*	Ts1:	texture coordinates for texture unit 1
*	Ns:		number of vertices
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::ComputeVerticesTexCoordsX(int Stage, float *Vs, float *Ts0, float *Ts1, int &Ns)
{
	int i = Stage;	
	int j, k, kk;
	float tllx, turx, llx, urx;
	int nBlocks  = powf(2.0, (float)(nButterfliesX - 1 - i));
	int Inc      = powf(2.0, (float)(i+1));
	int IncB     = powf(2.0, (float)(i)) - 1.0;
	int nHInputs = powf(2.0, (float)(i));
	float BlockSize = (float)(Width)/(float)nBlocks/2.0;
	float Index1Start, Index1End, Index2Start, Index2End;
	float Index1Start_, Index1End_, Index2Start_, Index2End_;
	float AngStart, AngEnd, AngScale;
	float v1, v2, Grad;
	float AngSign = (ForwardFFT_p) ? 1.0 : -1.0;

	if (i == 0) {
		AngStart = 0.0;
		AngEnd   = 0.0;
	}
	else {
		v1 = 0.0;
		v2 = nBlocks*IncB;
		Grad = (v2 - v1)/(float)(BlockSize-1);
		AngScale = -2.0*M_PI/(float)Width;
		AngStart = AngScale*(v1 - 0.5*Grad);
		AngEnd   = AngScale*(v2 + 0.5*Grad);
	}
	float *vptr  = Vs;
	float *tptr0 = Ts0;
	float *tptr1 = Ts1;
	
	Ns = 0;