fft.cpp

适用于MP3音乐中波形处理（FFT算法）

#include <windows.h>
#include <stdio.h>
#include <math.h>

#define SAMPLENUM		8
float PI=(float)3.1415926;

typedef struct// tagComplex
{
	float re; //复数的实部
	float im; //复数的虚部
}ST_Complex;

float RealData[SAMPLENUM];

ST_Complex inData[SAMPLENUM];
ST_Complex outData[SAMPLENUM];

void ChangeToComplex(int n,float *RealData,ST_Complex x[])
{
	int i;
	for(i=0;i<n;i++)
	{
		x[i].re = RealData[i];
		x[i].im = 0;
	}	
}
void init_x(int n,ST_Complex x[])
{
	int i;
	for(i=0;i<n;i++)
	{
		x[i].im = 0;
	}
	x[0].re=1;
	x[1].re=2;
	x[2].re=2;
	x[3].re=2;
	x[4].re=0;
	x[5].re=1;
	x[6].re=1;
	x[7].re=1;	
}

int GetInverse(int n,int preNum)
{
	int K = n/2;
	if(preNum<K)
	{
		return K+preNum;
	}
	else
	{
		return GetInverse(K,preNum-K);
	}
}

void ComplexAdd(ST_Complex *a,ST_Complex *b,ST_Complex *ret)
{
	ret->re = a->re + b->re;
	ret->im = a->im + b->im;
}

void ComplexSub(ST_Complex *a,ST_Complex *b,ST_Complex *ret)
{
	ret->re = a->re - b->re;
	ret->im = a->im - b->im;
}

void ComplexMul(ST_Complex *m,ST_Complex *n,ST_Complex *ret)
{
	float a,b,c;

	a = ( m->re - m->im ) * n->im;
	b = m->re * ( n->re - n->im );
	c = m->im * ( n->re + n->im );
	ret->re = a + b;
	ret->im = a + c;
}

void ComplexCopy(ST_Complex dest[],ST_Complex src[],int n)
{
	for (int i=0;i<n;i++)
	{
		dest[i].re = src[i].re;
		dest[i].im = src[i].im;
	}	
}
void TransToMod(ST_Complex src[],float mod[],int n)
{
	int i;	
	for(i=0;i<n;i++)
	{		
		if(abs(src[i].re) < 0.000001)
		{
			src[i].re = 0.0f;
		}
		if(abs(src[i].im) < 0.000001)
		{
			src[i].im = 0.0f;
		}
		mod[i] = (float)sqrt(src[i].re*src[i].re + src[i].im*src[i].im);
	}
}
/*****************************************
*
*	Name     :FFT
*   Function :Fast Fuliye Transformation
*   Params   :N -- Total count of sampling points
*             X -- Input sequence
*   Return   :XN(k)=sum[x(n)*Pow(e,j2*Pi/N)] 
*                   k,n:0..N-1
*
*
*****************************************/
void FFT(int n,ST_Complex X[],ST_Complex XK[])
{
	int i,j=0;
	int M,L;
	ST_Complex U={0,0};
	ST_Complex W={0,0};
	ST_Complex T={0,0};
	//ST_Complex *A = XK;	
	ST_Complex *A = (ST_Complex *)malloc(sizeof(ST_Complex) * n);	

	//Adjust sequence
	for(i=0;i<n;i++)
	{
		if(i==0)
		{
			A[0].re = X[0].re;
			A[0].im = X[0].im;
		}
		else
		{
			j=GetInverse(n,j);
			A[i].re = X[j].re;
			A[i].im = X[j].im;
		}
	}

	//确定级别数
	for(M=0;M<n;M++)
	{
		if((1<<M)==n)
			break;
	}

	for(L=1;L<=M;L++)//1-M级依次确定
	{
		int LE = (int)pow(2,L);//间隔
		int LE1 = LE/2;//W级数,如W0,W1,W2...

		W.re = (float)cos(PI/LE1);
		W.im = (float)-sin(PI/LE1);
#if 0
		if(abs(W.re)<1.0e-12)
			W.re = 0.0f;
		if(abs(W.im)<1.0e-12)
			W.im = 0.0f;
#endif	
		U.re = 1;
		U.im = 0;		
		for(j=0;j<LE1;j++)//
		{
			i=j;
			while(i<n)
			{
				int IP = i+LE1;
				ComplexMul(&A[IP],&U,&T);			//T=A[IP]*U;
				ComplexSub(&A[i],&T,&A[IP]);		//A[IP]=A[i]-T;//蝶形计算
				ComplexAdd(&A[i],&T,&A[i]);			//A[i]=A[i]+T;				
				i+=LE;
			}
			ComplexMul(&U,&W,&U);			//U=U*W;//不同的W次幂
		}		
	}
	ComplexCopy(X,A,n);	
	free(A);
}
void main()
{
	int i;
	init_x(SAMPLENUM,inData);	
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\ninData[%d]=%f + %f i",i,inData[i].re,inData[i].im);
	}
	printf("\n");
	FFT(SAMPLENUM,inData,NULL);
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\ninData[%d]=%f + %f i",i,inData[i].re,inData[i].im);
	}
	printf("\n");
	TransToMod(inData,RealData,SAMPLENUM);
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\nRealData[%d]=%f",i,RealData[i]);
	}	
	printf("\n");
}
/*
void main()
{
	int i;
	init_x(SAMPLENUM,inData);
	FFT(SAMPLENUM,inData,outData);
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\ninData[%d]=%f + %f i",i,inData[i].re,inData[i].im);
	}
	printf("\n");
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\noutData[%d]=%f + %f i",i,outData[i].re,outData[i].im);
	}	
	printf("\n");
	TransToMod(outData,RealData,SAMPLENUM);
	for(i=0;i<SAMPLENUM;i++)
	{
		printf("\nRealData[%d]=%f",i,RealData[i]);
	}	
	printf("\n");
}*/

/*
结果:

inData[0]=1.000000 + 0.000000 i
inData[1]=2.000000 + 0.000000 i
inData[2]=2.000000 + 0.000000 i
inData[3]=2.000000 + 0.000000 i
inData[4]=0.000000 + 0.000000 i
inData[5]=1.000000 + 0.000000 i
inData[6]=1.000000 + 0.000000 i
inData[7]=1.000000 + 0.000000 i

outData[0]=10.000000 + 0.000000 i
outData[1]=1.000000 + -2.414214 i
outData[2]=-2.000000 + 0.000000 i
outData[3]=1.000000 + -0.414213 i
outData[4]=-2.000000 + 0.000000 i
outData[5]=1.000000 + 0.414214 i
outData[6]=-2.000000 + 0.000000 i
outData[7]=1.000000 + 2.414213 i

*/
/*

inData[0]=1.000000 + 0.000000 i
inData[1]=2.000000 + 0.000000 i
inData[2]=2.000000 + 0.000000 i
inData[3]=2.000000 + 0.000000 i
inData[4]=0.000000 + 0.000000 i
inData[5]=1.000000 + 0.000000 i
inData[6]=1.000000 + 0.000000 i
inData[7]=1.000000 + 0.000000 i

inData[0]=10.000000 + 0.000000 i
inData[1]=1.000000 + -2.414214 i
inData[2]=-2.000000 + 0.000000 i
inData[3]=1.000000 + -0.414213 i
inData[4]=-2.000000 + 0.000000 i
inData[5]=1.000000 + 0.414214 i
inData[6]=-2.000000 + 0.000000 i
inData[7]=1.000000 + 2.414213 i

RealData[0]=10.000000
RealData[1]=2.613126
RealData[2]=2.000000
RealData[3]=1.082392
RealData[4]=2.000000
RealData[5]=1.082392
RealData[6]=2.000000
RealData[7]=2.613126
Press any key to continue
*/