dsp.cpp

《Visual C／C++图形图像与游戏编程典型实例解析》所有源代码。特别适合初学者。

///////////////////////////////////////////////////////////
//                                                      
//	类说明：				
//							
//		1。DFT_2D_FFT()		二维离散快速傅立叶变换						
//		2。DFT_2D_IFFT()	二维离散快速傅立叶逆变换						
//		3。DFT_FFT()		一维离散快速傅立叶变换						
//		4。DFT_IFFT()		一维离散快速傅立叶逆变换						
//      5。ReverseOrder()	倒序函数							
//
///////////////////////////////////////////////////////////

#include "stdafx.h"
#include "BMPChange.h"
#include "DSP.h"
#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CDSP

CDSP::CDSP()
{
}

CDSP::~CDSP()
{
}


BEGIN_MESSAGE_MAP(CDSP, CWnd)
	//{{AFX_MSG_MAP(CDSP)
		// NOTE - the ClassWizard will add and remove mapping macros here.
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()


/////////////////////////////////////////////////////////////////////////////
// CDSP message handlers

////////////////////////////////////////////////////
//          倒序函数                              //
////////////////////////////////////////////////////
void CDSP::ReverseOrder(Complex *A, int N)
{
	int NV2=N/2;
	int NM1=N-1;
	int I,J,K=0;
	Complex T;

	I=J=1;
	while(I<=NM1)
	{
		if(I<J)
		{
			T=A[J-1];
			A[J-1]=A[I-1];
			A[I-1]=T;
		}
		K=NV2;
		while(K<J)
		{
			J-=K;
			K/=2;
		}
		J+=K;
		I++;
	}
}

//////////////////////////////////////////////////////////////
//          时间抽选奇偶分解快速离散傅立叶变换（一维）      //
//////////////////////////////////////////////////////////////
void CDSP::DFT_FFT(Complex *A, int M)
{
	int N=(int)pow(2,M);
	ReverseOrder(A,N);
	for(int i=1;i<=M;i++)
	{
		int b=(int)pow(2,(i-1));
		for(int j=0;j<=(b-1);j++)
		{
			float p=(float)(pow(2,(M-i))*j*2.0*PI/(float)N);
			for(int k=j;k<=(N-1);)
			{
				float tr=(float)(A[k+b].Re*cos(p)+A[k+b].Im*sin(p));
				float ti=(float)(A[k+b].Im*cos(p)-A[k+b].Re*sin(p));
				A[k+b].Re=A[k].Re-tr;
				A[k+b].Im=A[k].Im-ti;
				A[k].Re+=tr;
				A[k].Im+=ti;
				k+=b*2;
			}
		}
	}
}

//////////////////////////////////////////////////////////////
//          时间抽选奇偶分解快速离散傅立叶变换（二维）      //
//////////////////////////////////////////////////////////////
void CDSP::DFT_2D_FFT(unsigned char *D, int W, int H,Complex* U)
{
	int WM=(int)(log(W)/log(2)+1.0f);
	int HM=(int)(log(H)/log(2)+1.0f);
	WM=HM=max(WM,HM);
	int WN=(int)pow(2,WM);
	int HN=(int)pow(2,HM);

	for(int i=0;i<HN;i++)
	{
		for(int j=0;j<WN*3;j++)
		{
			if(i<H && j<W*3)
			{
				U[i*WN*3+j].Re=D[i*W*3+j];
				U[i*WN*3+j].Im=0.0f;
			}
			else
				U[i*WN*3+j].Re=U[i*WN*3+j].Im=0.0f;
		}
	}

	Complex* UH=new Complex[HN];
	for(i=0;i<WN*3;i++)
	{
		for(int j=0;j<HN;j++)
		{
			UH[j].Re=U[j*WN*3+i].Re;
			UH[j].Im=U[j*WN*3+i].Im;
		}
		DFT_FFT(UH,HM);
		for(j=0;j<HN;j++)
		{
			U[j*WN*3+i].Re=HN*UH[j].Re;
			U[j*WN*3+i].Im=HN*UH[j].Im;
		}
	}
	delete []UH;

	Complex* UW=new Complex[WN];
	for(i=0;i<HN;i++)
	{
		for(int k=0;k<3;k++)
		{
			for(int j=0;j<WN;j++)
			{
				UW[j].Re=U[i*WN*3+j*3+k].Re;
				UW[j].Im=U[i*WN*3+j*3+k].Im;
			}
			DFT_FFT(UW,WM);
			for(j=0;j<WN;j++)
			{
				U[i*WN*3+j*3+k].Re=UW[j].Re;
				U[i*WN*3+j*3+k].Im=UW[j].Im;
			}
		}
	}
	delete []UH;
}

//////////////////////////////////////////////////////////////
//          时间抽选奇偶分解快速离散傅立叶逆变换（一维）    //
//////////////////////////////////////////////////////////////
void CDSP::DFT_IFFT(Complex *A, int M)
{
	int N=(int)pow(2,M);
	
	ReverseOrder(A,N);

	for(int i=1;i<=M;i++)
	{
		int b=(int)pow(2,(i-1));
		for(int j=0;j<=(b-1);j++)
		{
			float p=(float)(pow(2,(M-i))*j*2.0*PI/(float)N);
			for(int k=j;k<=(N-1);)
			{
				float tr=(float)(A[k+b].Re*cos(p)-A[k+b].Im*sin(p));
				float ti=(float)(A[k+b].Im*cos(p)+A[k+b].Re*sin(p));
				A[k+b].Re=A[k].Re-tr;
				A[k+b].Im=A[k].Im-ti;
				A[k].Re+=tr;
				A[k].Im+=ti;
				k+=b*2;
			}
		}
	}

	for(i=0;i<N;i++)
	{
		A[i].Re/=N;
		A[i].Im/=N;
	}
}

//////////////////////////////////////////////////////////////
//          时间抽选奇偶分解快速离散傅立叶逆变换（二维）    //
//////////////////////////////////////////////////////////////
void CDSP::DFT_2D_IFFT(unsigned char *D, int W, int H, Complex *U)
{
	int WM=(int)(log(W)/log(2)+1.0f);
	int HM=(int)(log(H)/log(2)+1.0f);
	WM=HM=max(WM,HM);
	int WN=(int)pow(2,WM);
	int HN=(int)pow(2,HM);

	Complex* UW=new Complex[WN];
	for(int i=0;i<HN;i++)
	{
		for(int k=0;k<3;k++)
		{
			for(int j=0;j<WN;j++)
			{
				UW[j].Re=U[i*WN*3+j*3+k].Re;
				UW[j].Im=U[i*WN*3+j*3+k].Im;
			}
			DFT_IFFT(UW,WM);
			for(j=0;j<WN;j++)
			{
				U[i*WN*3+j*3+k].Re=UW[j].Re;
				U[i*WN*3+j*3+k].Im=UW[j].Im;
			}
		}
	}
	delete []UW;

	Complex* UH=new Complex[HN];
	for(i=0;i<WN*3;i++)
	{
		for(int j=0;j<HN;j++)
		{
			UH[j].Re=U[j*WN*3+i].Re;
			UH[j].Im=U[j*WN*3+i].Im;
		}
		DFT_IFFT(UH,HM);
		for(j=0;j<HN;j++)
		{
			U[j*WN*3+i].Re=UH[j].Re/HN;
			U[j*WN*3+i].Im=UH[j].Im/HN;
		}
	}
	delete []UH;

	for(i=0;i<HN;i++)
	{
		for(int j=0;j<WN*3;j++)
		{
			if(i<H && j<W*3)
				D[i*W*3+j]=U[i*WN*3+j].Re;
		}
	}
}