fft.cpp

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

*   N:	The number of elements
*   Ar:	The array
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::PrintArrray(int N, float *Ar)
{
	int i;

	for (i = 0; i < N; i++) {
		cerr << i << " " << Ar[i] << endl;
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::GetBestCutOff
*
* DESCRIPTION:
*	Computes the best transition point from method 1 to method 2.
*
* FORMAL PARAMETERS:
*   N:	The number of elements
*   Ar:	The array
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
int FFT::GetBestCutOff(int N, float *Ar)
{
	int i;
	int index = -1;
	float maxf = -1.0;

	for (i = 0; i < N; i++) {
		if (maxf < Ar[i]) {
			index = i;
			maxf = Ar[i];
		}
	}
	return(index);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::SetGPUVendor
*
* DESCRIPTION:
*	Set the vendor
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	The GPU vendor
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
Vendor FFT::SetGPUVendor(void)
{
	char *str = strdup((char *)glGetString(GL_VENDOR));
	if (strstr(str, "ATI") != NULL) {
		GPU = GPU_ATI;
	}
	else {
		if (strstr(str, "NVIDIA") != NULL) {
			GPU = GPU_NVIDIA;
		}
		else {
			GPU = GPU_NONE;
		}
	}

	if (GPU == GPU_NVIDIA) {
		cerr << "NVIDIA board found" << endl;
		str = strdup((char *)glGetString(GL_RENDERER));
	}
	else {
		cerr << "Currently only implemented on NVIDIA" << endl;
		exit(0);
	}

	return(GPU);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::DO_FFT_GPU_X
*
* DESCRIPTION:
*	Do the FFT in x only
*
* FORMAL PARAMETERS:
*	FirstTime_p:	true for the very first time only
*   Texs:			an array containig (texR1, texI1, texR2, texI2), where
*	texR1:			texture for the real part of image 1
*	texI1:			texture for the imaginary part of image 1
*	texR2:			texture for the real part of image 2
*	texI2:			texture for the imaginary part of image 2
*
* RETURNS:
*	FirstTime_p
*   Texs
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DO_FFT_GPU_X(bool&FirstTime_p, GLuint *Texs)
{
	int i = 0;

	for (i = 0; i < nButterfliesXWorking; i++) {
		CurrentButterflyStage = i;

		RenderFFTStageX(FirstTime_p,
			texButterflyLookupI_X[i],  
			texButterflyLookupWR_X[i],  
			texButterflyLookupWI_X[i], Texs);

		FirstTime_p = false;
		SetTexs(Texs, texTmpR1, texTmpI1, texTmpR2, texTmpI2);
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::DO_FFT_GPU_Y
*
* DESCRIPTION:
*	Do the FFT in y only
*
* FORMAL PARAMETERS:
*	FirstTime_p:	true for the very first time only
*   Texs:			an array containig (texR1, texI1, texR2, texI2), where
*	texR1:			texture for the real part of image 1
*	texI1:			texture for the imaginary part of image 1
*	texR2:			texture for the real part of image 2
*	texI2:			texture for the imaginary part of image 2
*
* RETURNS:
*	FirstTime_p
*   Texs
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DO_FFT_GPU_Y(bool&FirstTime_p, GLuint *Texs)
{
	int i = 0;

	for (i = 0; i < nButterfliesYWorking; i++) {
		CurrentButterflyStage = i;
		RenderFFTStageY(FirstTime_p,
			texButterflyLookupI_Y[i],  
			texButterflyLookupWR_Y[i],  
			texButterflyLookupWI_Y[i], Texs);

		FirstTime_p = false;
		SetTexs(Texs, texTmpR1, texTmpI1, texTmpR2, texTmpI2);
	}
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::DO_FFT_GPU_Y
*
* DESCRIPTION:
*	Do the FFT in y only
*
* FORMAL PARAMETERS:
*	FirstTime_p:	true for the very first time only
*   Texs:			an array containig (texR1, texI1, texR2, texI2), where
*	texR1:			texture for the real part of image 1
*	texI1:			texture for the imaginary part of image 1
*	texR2:			texture for the real part of image 2
*	texI2:			texture for the imaginary part of image 2
*
* RETURNS:
*	FirstTime_p
*   Texs
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DO_FFT_GPU(bool&FirstTime_p, GLuint *Texs)
{
	DO_FFT_GPU_X(FirstTime_p, Texs);
	DO_FFT_GPU_Y(FirstTime_p, Texs);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::DoFFT
*
* DESCRIPTION:
*	Do the FFT in x, y or xy
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::DoFFT(void)
{
	GLuint Texs[4];
	bool FirstTime_p;

	Buffers->Bind();

	switch (type) {
	case X:
		FirstTime_p = true;
		SetTexs(Texs, texReal1, texImag1, texReal2, texImag2);
		DO_FFT_GPU_X(FirstTime_p, Texs);
		break;
	case Y:
		FirstTime_p = true;
		SetTexs(Texs, texReal1, texImag1, texReal2, texImag2);
		DO_FFT_GPU_Y(FirstTime_p, Texs);
		break;
	case XY:
		FirstTime_p = true;
		SetTexs(Texs, texReal1, texImag1, texReal2, texImag2);
		DO_FFT_GPU(FirstTime_p, Texs);
		break;
	default:
		break;
	}
	Buffers->Unbind();
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::ComputeFrameRate
*
* DESCRIPTION:
*	Computes the frame rate.
*
* FORMAL PARAMETERS:
*	Count:		Number of frames processed
*	thisClock:	Current clock value
*	prevClock:	Previous clock value
*
* RETURNS:
*	FrameRate:	The computed frame rate
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
float FFT::ComputeFrameRate(int Count, long thisClock, long prevClock)
{
	int NumSimultanuesComplex2DImages = 2;

	float duration  = (double)(thisClock - prevClock)/ CLOCKS_PER_SEC;
	float FrameRate = ((float)Count*NumSimultanuesComplex2DImages)/duration;
	return(FrameRate);
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::FindOptimalTransitionPoints
*
* DESCRIPTION:
*	Vary the xCutOff and yCutOff sequentially for all possible combinations
*	and picks the best values for the fastest frame rate. xCutOff is the
*	transition point between doing FFT in x using method 1 and method 2. 
*	Similarly, yCutOff.
*
* FORMAL PARAMETERS:
*	none
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::FindOptimalTransitionPoints(void)
{
	long thisClock, startClock;

	int numTimes        = 0;
	int avgCount        = 100;
	float FrameRate     = 0.0;
	float PrevFrameRate = 0.0;
	bool FirstTimeX_p   = true;
	bool DoneX_p        = false;
	bool DoneAdjX_p     = false;
	bool FirstTimeY_p   = true;
	bool DoneY_p        = false;
	bool DoneAdjY_p     = false;

	startClock = clock();

	while (!DoneAdjX_p || !DoneAdjY_p) {
		if (FirstTimeX_p && (numTimes == 0)) {
			cerr << "*** Optimizing Frame Rate ***" << endl;
			cerr << "Do not touch anything until finished" << endl;
			cerr << "Will be done in a few minutes" << endl;
			cerr << endl;
		}

		if (numTimes == avgCount) {
			thisClock = clock();
			FrameRate = ComputeFrameRate(avgCount, thisClock, startClock);
			startClock = thisClock;
			numTimes = 0;

			if (FirstTimeX_p) {
				cerr << "Detecting best X frame rate...";
				cerr << "do not touch anything..." << endl;
				xCutOff = 0;
			}
			else {
				if (!DoneX_p) {
					cerr << xCutOff << " " << FrameRate << endl;
					xFrameRate[xCutOff] = FrameRate;
					xCutOff++;
					DoneX_p = (xCutOff >= nButterfliesX) || (PrevFrameRate > FrameRate);
					PrevFrameRate = FrameRate;
				}
			}
			if (DoneX_p) {
				if (!DoneAdjX_p) {
					xCutOff = GetBestCutOff(nButterfliesX, xFrameRate);
					cerr << "Best xCutOff: " << xCutOff << endl;
					cerr << "Done X" << endl;
					DoneAdjX_p = true;
				}
			}
			FirstTimeX_p = false;

			if (DoneAdjX_p) {
				if (FirstTimeY_p) {
					cerr << "Detecting best Y frame rate...";
					cerr << "do not touch anything..." << endl;
					yCutOff = 0;
					PrevFrameRate = 0.0;
				}
				else {
					if (!DoneY_p) {
						cerr << yCutOff << " " << FrameRate << endl;
						yFrameRate[yCutOff] = FrameRate;
						yCutOff++;
						DoneY_p = (yCutOff >= nButterfliesY) || (PrevFrameRate > FrameRate);
						PrevFrameRate = FrameRate;
					}
				}
				if (DoneY_p) {
					if (!DoneAdjY_p) {
						yCutOff = GetBestCutOff(nButterfliesY, yFrameRate);
						cerr << "Best yCutOff: " << yCutOff << endl;
						cerr << "Done Y" << endl;
						cerr << "*** DONE Optimizing Frame Rate ***" << endl;
						DoneAdjY_p = true;
					}
				}
				FirstTimeY_p = false;
			}

			if (DoneAdjX_p && DoneAdjY_p) {
				cerr << "Optimal Frame Rate = " << FrameRate << endl;
			}
		}
		DoFFT();
		numTimes++;
	}
	cerr << endl;
}

/*************************[MAN-BEG]*****************************************
*
* NAME:
*   FFT::Print
*
* DESCRIPTION:
*	Prints a string to the display window.
*	Usage: Print(X, Y, "My string", GLUT_BITMAP_HELVETICA_12);
*
* FORMAL PARAMETERS:
*	X:		x coordinate
*	Y:		y coordinate
*	str:	the string to print
*	font:	font to use (e.g.: GLUT_BITMAP_HELVETICA_12)
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::Print(float X, float Y, char *str, void *font)
{
    int i; 
    int fontWidth = glutBitmapWidth(font, 'A'); 
    
	glRasterPos2f(X, Y); 
    for(i = 0; str[i] != '\0'; i++) 
        glutBitmapCharacter(font, str[i]);
} 

/*************************[MAN-BEG]*****************************************
*
* NAME:
*	FFT::SetTexs
*
* DESCRIPTION:
*	Sets Texs with 4 texture names.
*
* FORMAL PARAMETERS:
*	texR1, texI1,texR2, texI2: the four texture names
*	Texs: an array
*
* RETURNS:
*	none
*
* REVISION HISTORY:
* Rev     When      Who         What
* V1      15Dec2004 Thilaka     Created.
**************************[MAN-END]*****************************************/
void FFT::SetTexs(GLuint* Texs, GLuint texR1, GLuint texI1, GLuint texR2, GLuint texI2)
{
	int i = 0;
	Texs[i] = texR1; i++;
	Texs[i] = texI1; i++;
	Texs[i] = texR2; i++;
	Texs[i] = texI2; i++;
}