fftengine.cpp

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

	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;
		
			tllx = -0.5 + kk*BlockSize;
			turx = tllx + BlockSize;
			
			llx = kk*BlockSize;
			urx = llx + BlockSize;

			// 1st vertex
			*(tptr0++) = tllx;
			*(tptr0++) = -0.5;
			*(tptr0++) = Index1Start_-0.5;
			*(tptr0++) = Index2Start_-0.5;
			
			*(tptr1++) = AngSign*(AngStart+j*M_PI);

			*(vptr++) = llx;
			*(vptr++) = 0.0;
			Ns++;
					
			// 2nd vertex
			*(tptr0++) = turx;
			*(tptr0++) = -0.5;
			*(tptr0++) = Index1End_+0.5;
			*(tptr0++) = Index2End_+0.5;
			
			*(tptr1++) = AngSign*(AngEnd+j*M_PI);

			*(vptr++) = urx;
			*(vptr++) = 0.0;
			Ns++;
											
			// 3rd vertex
			*(tptr0++) = turx;
			*(tptr0++) = Height-1+0.5;
			*(tptr0++) = Index1End_+0.5;
			*(tptr0++) = Index2End_+0.5;

			*(tptr1++) = AngSign*(AngEnd+j*M_PI);

			*(vptr++) = urx;
			*(vptr++) = Height;
			Ns++;
			
			// 4th vertex
			*(tptr0++) = tllx;
			*(tptr0++) = Height-1+0.5;
			*(tptr0++) = Index1Start_-0.5;
			*(tptr0++) = Index2Start_-0.5;
			
			*(tptr1++) = AngSign*(AngStart+j*M_PI);

			*(vptr++) = llx;
			*(vptr++) = Height;
			Ns++;
		}	
		Index1Start += Inc;
		Index2Start += Inc;
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::GenDisplayListsX
*
* DESCRIPTION:
*	Generates display lists containing draw arrays for x FFT method 2
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::GenDisplayListsX(void)
{
	int Ns = 0, i;
	int nButterflies = nButterfliesX;
	int nBlocks  = powf(2.0, (float)(nButterflies - 1));
	int SizeV   = nBlocks*2*4*2;
	int SizeT   = nBlocks*2*4*4;

	float *Vs, *Ts0, *Ts1;
	Vs  = new float [SizeV];
	Ts0 = new float [SizeT];
	Ts1 = new float [SizeT];

	ListX = glGenLists(nButterflies);
	for (i = 0; i < nButterflies; i++) {
		memset(Vs, 0, sizeof(float)*SizeV);
		memset(Ts0, 0, sizeof(float)*SizeT);
		memset(Ts1, 0, sizeof(float)*SizeT);

		ComputeVerticesTexCoordsX(i, Vs, Ts0, Ts1, Ns); 

		glEnableClientState(GL_VERTEX_ARRAY); 
		glVertexPointer(2, GL_FLOAT, 0, Vs);

		glClientActiveTextureARB(GL_TEXTURE0_ARB);
		glEnableClientState(GL_TEXTURE_COORD_ARRAY);
		glTexCoordPointer(4 , GL_FLOAT, 0, Ts0); 
		
		glClientActiveTextureARB(GL_TEXTURE1_ARB);
		glEnableClientState(GL_TEXTURE_COORD_ARRAY);
		glTexCoordPointer(1 , GL_FLOAT, 0, Ts1);

		glNewList(ListX+i, GL_COMPILE);
		glDrawArrays(GL_QUADS, 0, Ns);
		glEndList();
		glDisableClientState(GL_VERTEX_ARRAY);
		glClientActiveTextureARB(GL_TEXTURE0_ARB);
		glDisableClientState(GL_TEXTURE_COORD_ARRAY);
		glClientActiveTextureARB(GL_TEXTURE1_ARB);
		glDisableClientState(GL_TEXTURE_COORD_ARRAY);
	}
	delete [] Vs;
	delete [] Ts0;
	delete [] Ts1;
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::ComputeVerticesTexCoordsY
*
* 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::ComputeVerticesTexCoordsY(int Stage, float *Vs, float *Ts0, float *Ts1, int &Ns)
{
	int i = Stage;	
	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);
	}
	float *vptr  = Vs;
	float *tptr0 = Ts0;
	float *tptr1 = Ts1;

	Ns = 0;
	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;
		
			tlly = -0.5 + kk*BlockSize;
			tury = tlly + BlockSize;
			
			lly = kk*BlockSize;
			ury = lly + BlockSize;

			// 1st vertex
			*(tptr0++) = -0.5;
			*(tptr0++) = tlly;
			*(tptr0++) = Index1Start_-0.5;
			*(tptr0++) = Index2Start_-0.5;
			
			*(tptr1++) = AngSign*(AngStart+j*M_PI);

			*(vptr++) = 0.0;
			*(vptr++) = lly;
			Ns++;
									
			// 2nd vertex
			*(tptr0++) = -0.5;
			*(tptr0++) = tury;
			*(tptr0++) = Index1End_+0.5;
			*(tptr0++) = Index2End_+0.5;
			
			*(tptr1++) = AngSign*(AngEnd+j*M_PI);

			*(vptr++) = 0.0;
			*(vptr++) = ury;
			Ns++;
							
			// 3rd vertex				
			*(tptr0++) = Width-1+0.5;
			*(tptr0++) = tury;
			*(tptr0++) = Index1End_+0.5;
			*(tptr0++) = Index2End_+0.5;

			*(tptr1++) = AngSign*(AngEnd+j*M_PI);

			*(vptr++) = Width;
			*(vptr++) = ury;
			Ns++;
					
			// 4th vertex
			*(tptr0++) = Width-1+0.5;
			*(tptr0++) = tlly;
			*(tptr0++) = Index1Start_-0.5;
			*(tptr0++) = Index2Start_-0.5;

			*(tptr1++) = AngSign*(AngStart+j*M_PI);

			*(vptr++) = Width;
			*(vptr++) = lly;
			Ns++;
		}	
		Index1Start += Inc;
		Index2Start += Inc;
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::GenDisplayListsY
*
* DESCRIPTION:
*	Generates display lists containing draw arrays for y FFT method 2
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::GenDisplayListsY(void)
{
	int Ns = 0, i;
	int nButterflies = nButterfliesY;
	int nBlocks  = powf(2.0, (float)(nButterflies - 1));
	int SizeV   = nBlocks*2*4*2;
	int SizeT   = nBlocks*2*4*4;

	float *Vs, *Ts0, *Ts1;
	Vs  = new float [SizeV];
	Ts0 = new float [SizeT];
	Ts1 = new float [SizeT];

	ListY = glGenLists(nButterflies);
	for (i = 0; i < nButterflies; i++) {
		memset(Vs, 0, sizeof(float)*SizeV);
		memset(Ts0, 0, sizeof(float)*SizeT);
		memset(Ts1, 0, sizeof(float)*SizeT);

		ComputeVerticesTexCoordsY(i, Vs, Ts0, Ts1, Ns);

		glEnableClientState(GL_VERTEX_ARRAY); 
		glVertexPointer(2, GL_FLOAT, 0, Vs);

		glClientActiveTextureARB(GL_TEXTURE0_ARB);
		glEnableClientState(GL_TEXTURE_COORD_ARRAY);
		glTexCoordPointer(4 , GL_FLOAT, 0, Ts0); 
		
		glClientActiveTextureARB(GL_TEXTURE1_ARB);
		glEnableClientState(GL_TEXTURE_COORD_ARRAY);
		glTexCoordPointer(1 , GL_FLOAT, 0, Ts1);

		glNewList(ListY+i, GL_COMPILE);
		glDrawArrays(GL_QUADS, 0, Ns);
		glEndList();

		glDisableClientState(GL_VERTEX_ARRAY);
		glClientActiveTextureARB(GL_TEXTURE0_ARB);
		glDisableClientState(GL_TEXTURE_COORD_ARRAY);
		glClientActiveTextureARB(GL_TEXTURE1_ARB);
		glDisableClientState(GL_TEXTURE_COORD_ARRAY);
	}
	delete [] Vs;
	delete [] Ts0;	
	delete [] Ts1;	
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadForFFT
*
* DESCRIPTION:
*	Draws a quad
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DrawQuadForFFT(void)
{
	glEnable(GL_TEXTURE_RECTANGLE_NV);	
	glBegin(GL_QUADS);
	glTexCoord2f(-0.5, -0.5);
	glVertex2f(0.0, 0.0);
	glTexCoord2f(Width-1+0.5, -0.5);
	glVertex2f(Width, 0.0);
	glTexCoord2f(Width-1+0.5, Height-1+0.5);
	glVertex2f(Width, Height);
	glTexCoord2f(-0.5, Height-1+0.5);
	glVertex2f(0.0, Height);
	glEnd();	
	glDisable(GL_TEXTURE_RECTANGLE_NV);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::GenQuadDisplayList
*
* DESCRIPTION:
*	Generates a display list for drawing a quad for method 1
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::GenQuadDisplayList(void)
{
	glEnable(GL_TEXTURE_RECTANGLE_NV);	
	glBegin(GL_QUADS);
	glTexCoord2f(-0.5, -0.5);
	glVertex2f(0.0, 0.0);
	glTexCoord2f(Width-1+0.5, -0.5);
	glVertex2f(Width, 0.0);
	glTexCoord2f(Width-1+0.5, Height-1+0.5);
	glVertex2f(Width, Height);
	glTexCoord2f(-0.5, Height-1+0.5);
	glVertex2f(0.0, Height);
	glEnd();	
	glDisable(GL_TEXTURE_RECTANGLE_NV);

	float *Vs, *Ts;
	int i;
	int Size   = 2*4;
	Vs = new float [Size];
	Ts = new float [Size];

	QList = glGenLists(1);
	memset(Vs, 0, sizeof(float)*Size);
	memset(Ts, 0, sizeof(float)*Size);

	glTexCoord2f(-0.5, -0.5);
	glVertex2f(0.0, 0.0);
	glTexCoord2f(Width-1+0.5, -0.5);
	glVertex2f(Width, 0.0);
	glTexCoord2f(Width-1+0.5, Height-1+0.5);
	glVertex2f(Width, Height);
	glTexCoord2f(-0.5, Height-1+0.5);
	glVertex2f(0.0, Height);

	i = 0;
	Vs[i] = 0.0; i++;
	Vs[i] = 0.0; i++;

	Vs[i] = Width; i++;
	Vs[i] = 0.0; i++;

	Vs[i] = Width; i++;
	Vs[i] = Height; i++;

	Vs[i] = 0.0; i++;
	Vs[i] = Height; i++;

	i = 0;
	Ts[i] = -0.5; i++;
	Ts[i] = -0.5; i++;

	Ts[i] = Width-1+0.5; i++;
	Ts[i] = -0.5; i++;

	Ts[i] = Width-1+0.5; i++;
	Ts[i] = Height-1+0.5; i++;

	Ts[i] = -0.5; i++;
	Ts[i] = Height-1+0.5; i++;


	glEnableClientState(GL_VERTEX_ARRAY); 
	glVertexPointer(2, GL_FLOAT, 0, Vs);

	glClientActiveTextureARB(GL_TEXTURE0_ARB);
	glEnableClientState(GL_TEXTURE_COORD_ARRAY);
	glTexCoordPointer(2 , GL_FLOAT, 0, Ts); 

	glNewList(QList, GL_COMPILE);
	glDrawArrays(GL_QUADS, 0, 4);
	glEndList();

	glDisableClientState(GL_VERTEX_ARRAY);
	glClientActiveTextureARB(GL_TEXTURE0_ARB);
	glDisableClientState(GL_TEXTURE_COORD_ARRAY);

	delete [] Vs;
	delete [] Ts;		
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::DrawQuadForFFT_List
*
* DESCRIPTION:
*	Draws a quad for method 1 in doing FFT 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::DrawQuadForFFT_List(void)
{
	glActiveTextureARB(GL_TEXTURE0_ARB);
	glEnable(GL_TEXTURE_RECTANGLE_NV);
	glCallList(QList);
	glDisable(GL_TEXTURE_RECTANGLE_NV);
}