dft.cpp

主要是应用VC进行傅立叶变换和反变换的程序

#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "dft.h"

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

#define PIE 3.14159265

//////////////////////////////////////////////////

Complex::Complex()
{
	Re=0.0;
	Im=0.0;
}

Complex::~Complex()
{
}

Complex::Complex(double x,double y)
{
	Re=x;Im=y;
}

double Complex::abs()
{
	return sqrt(Re*Re+Im*Im);
}

double Complex::GetRe()
{
	return Re;
}

double Complex::GetIm()
{
	return Im;
}

void Complex::operator=(wComplex &cm)
{
	Re=cm.GetRe();
	Im=cm.GetIm();
}

void Complex::operator+=(wComplex &cm)
{	
	Re+=cm.GetRe();
	Im+=cm.GetIm();
}

void Complex::operator-=(wComplex &cm)
{
	Re-=cm.GetRe();
	Im-=cm.GetIm();
}

void Complex::operator*=(wComplex &cm)
{
	double x=Re*cm.GetRe()-Im*cm.GetIm();
	double y=Re*cm.GetIm()+Im*cm.GetRe();
	Re=x;
	Im=y;
}

Complex Complex::operator*(wComplex &cm)
{
	Complex temp;
	temp.Re=Re*cm.GetRe()-Im*cm.GetIm();
	temp.Im=Re*cm.GetIm()+Im*cm.GetRe();
	return temp;
}

void Complex::operator*=(float var)
{
	Re*=var;
	Im*=var;
}

wComplex Complex::operator+(wComplex &cm)
{
	double x=Re+cm.GetRe();
	double y=Im+cm.GetIm();
	return wComplex(x,y);
}

void Complex::operator/=(double x)
{
	Re/=x;Im/=x;
}

Complex Complex::operator/(double x)
{
	Complex temp;
	temp.Re=Re/x;
	temp.Im=Im/x;
	return temp;
}

wComplex Complex::operator-(wComplex &cm)
{
	double x=Re-cm.GetRe();
	double y=Im-cm.GetIm();
	return wComplex(x,y);
}

//////////////////////////////////////////////////

Fourior::Fourior()
{
	flag=FALSE;
}

Fourior::Fourior(int N)
{
	flag=FALSE;
	ByteNum=N;
	bWn=new BYTE[ByteNum*sizeof(wComplex)];
	Wn=(wComplex *)bWn;
	if(!Wn)return;
	for(int i=0;i<ByteNum;i++)
	{
		double w=-2.0*PIE*i/ByteNum;
		*(Wn+i)=wComplex(cos(w),sin(w));
	}
	flag=TRUE;
}

void Fourior::SetInverse()
{
	if(flag==TRUE)
	{
		for(int i=0;i<ByteNum;i++)
		{
			double x=(Wn+i)->GetRe();
			double y=-(Wn+i)->GetIm();
			*(Wn+i)=wComplex(x,y);
		}
	}
}

Fourior::~Fourior()
{
	if(flag==TRUE)delete bWn;
}

void Fourior::DFT(wComplex *Input)
{
	int i,j;
	BYTE *btemp=new BYTE[ByteNum*sizeof(wComplex)];
	wComplex *temp=(wComplex *)btemp;
	if(!temp)
	{
		AfxMessageBox("wrf");
		return;
	}

	wComplex var0((double)0,(double)0);
	for(i=0;i<ByteNum;i++)
	{
		*(temp+i)=var0;
		for(j=0;j<ByteNum;j++)
		{
			wComplex temp1=*(Input+j);
			temp1*=*(Wn+(j*i)%ByteNum);
			*(temp+i)+=temp1;
		}
	}
	for(i=0;i<ByteNum;i++)
	{
		*(Input+i)=*(temp+i);
	}
	delete btemp;
}

BOOL Fourior::FFT(wComplex *Input)
{
	int i,j,k,l;
	BYTE *btemp=new BYTE[ByteNum*sizeof(wComplex)];
	wComplex *output=(wComplex*)btemp;
	int M;
	int N1,N2,N3,N4;
	
	int N=ByteNum;
	int temp=N;i=0;
	while(1)
	{
		temp>>=1;i++;
		if(temp&1==1)break;
	}
	M=i;
	BitsNum=i;
	if(ByteNum!=1<<BitsNum)
	{
		AfxMessageBox("FFT要求数据长度为2的n幂次方可");
		return FALSE;
	}
	for(i=0;i<N;i++)   //数据的逆序
	{
		k=0; l=i;
		for(j=0;j<M;j++)
		{
			k<<=1;
			k|=(l&1);
			l>>=1;
		}
		*(output+i)=*(Input+k);
	}
	for(i=1;i<=M;i++)     //蝶形运算
	{
		N1=N;
		for(j=0;j<i;j++)N1>>=1;
		N2=1;
		for(j=0;j<i-1;j++)N2<<=1;
		N3=N2*2;
		N4=1;
		for(j=0;j<M-i;j++)N4<<=1;
		for(j=0;j<N1;j++)
		{
			for(k=0;k<N2;k++)
			{
				wComplex temp=*(output+j*N3+k+N2);
				temp*=Wn[(k*N4)%N];
				*(Input+j*N3+k)=*(output+j*N3+k)+temp;
				*(Input+j*N3+k+N2)=*(output+j*N3+k)-temp;
			}
		}
		for(j=0;j<N;j++)*(output+j)=*(Input+j);
	}
	delete btemp;
	return TRUE;
}