www.pudn.com > AudioWave.rar > FFT.cpp


// FFT.cpp: implementation of the CFFT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FFT.h"
#include "math.h"
#include "complex.h"
#include "define.h"
#include "math.h"
#include "malloc.h"
#include "conio.h"
#include "stdio.h"
#include "string.h"
#define pi (double)3.14159265359
#define m 6

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFFT::CFFT()
{
}

CFFT::~CFFT()
{

}

void CFFT::FFTSignal(
int nPower,
double *RealIn,
double *ImagIn,
double *RealOut,
double *ImagOut )
{
int i;
int nSamples = 1<<nPower;
COMPLEX *TV,*FV;
TV = new COMPLEX[nSamples];
FV = new COMPLEX[nSamples];
for(i=0; i<nSamples; i++)
{
TV[i].re = RealIn[i];
TV[i].im = ImagIn[i];
}
FFT(TV, FV, nPower);
for(i=0; i<nSamples; i++)
{
RealOut[i] = FV[i].re/sqrt((double)nSamples);
ImagOut[i] = FV[i].im/sqrt((double)nSamples);
}
delete[] TV;
delete[] FV;
}

void CFFT::IFFTSignal(
int nPower,
double *RealIn,
double *ImagIn,
double *RealOut,
double *ImagOut )
{
int i;
int nSamples = 1<<nPower;
COMPLEX *TV,*FV;
TV = new COMPLEX[nSamples];
FV = new COMPLEX[nSamples];

for(i=0; i<nSamples; i++)
{
FV[i].re = RealOut[i];
FV[i].im = ImagOut[i];
}

IFFT(FV, TV, nPower);

for(i=0; i<nSamples; i++)
{
RealIn[i] = TV[i].re;
ImagIn[i] = TV[i].im;
}
delete[] TV;
delete[] FV;
}



//复数加运算
COMPLEX Add(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re+c2.re;
c.im=c1.im+c2.im;
return c;
}

//复数乘运算
COMPLEX Mul(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re*c2.re-c1.im*c2.im;
c.im=c1.re*c2.im+c2.re*c1.im;
return c;
}

//复数减运算
COMPLEX Sub(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re-c2.re;
c.im=c1.im-c2.im;
return c;
}

//快速傅里叶变换
void CFFT::FFT(COMPLEX *TD,COMPLEX *FD,int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
COMPLEX *W,*X1,*X2,*X;

count=1<<power;

W = new COMPLEX[count];
X1 = new COMPLEX[count];
X2 = new COMPLEX[count];

for (i=0;i<count/2;i++)
{
angle=-i*pi*2/count;
W[i].re=cos(angle);
W[i].im=sin(angle);
}

memcpy(X1,TD,sizeof(COMPLEX)*count);

for (k=0;k<power;k++)
{
for (j=0;j<1<<k;j++)
{
bfsize=1<<(power-k);
for (i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
}
}
X=X1;
X1=X2;
X2=X;
}

for (j=0;j<count;j++)
{
p=0;
for (i=0;i<power;i++)
{
if (j&amt;(1<<i))
p+=1<<(power-i-1);
}
FD[j]=X1[p];
}

delete[] W;
delete[] X1;
delete[] X2;
}

//逆变换
void CFFT::IFFT(COMPLEX *FD,COMPLEX *TD,int power)
{
int i,count;
COMPLEX *x;

count=1<<power;
x = new COMPLEX[count];
memcpy(x,FD,sizeof(COMPLEX)*count);

for (i=0;i<count;i++)
{
x[i].im=-x[i].im;
}
FFT(x,TD,power);

for(i=0;i<count;i++)
{
TD[i].re/=sqrt((double)count);
TD[i].im=-TD[i].im/sqrt((double)count);
}
delete[] x;
}