fft.c

用c语言编写实现FFT算法

#include "main_leak.h"
#include "math.h"

//以输入序列的偶数点和奇数点作为实虚部，然后根据共轭对称求得原序列的变换
//m为对n取以2为底的对数值（2的m次方为n）
FFT(float FFTDATA[],float FFTIMAGE[],int m,int n)
{
	float *leftdata_p;
	float *ibufferfft_p=FFTIMAGE;
	float *rcgene_p,*icgene_p;
	float rpartkb,ipartkb;					//X（K+B）的实部与虚部
	int i,j,l,b,p,k;
	float pi=3.1415926;
	leftdata_p=FFTDATA;


//旋转因子实虚部暂存区
 	rcgene_p=BUFFER;
 	icgene_p=SECONDBUFFER;
 
//倒序
	changeorder(FFTDATA,FFTIMAGE,n);



//计算旋转因子
	for(i=0;i<n;i++)
		{
			*(rcgene_p+i)=cos(i*2.0*pi/n);
			*(icgene_p+i)=sin(i*2.0*pi/n);
		}
		
		
		
//继续进行FFT运算（运算的是以输入序列的偶数点作为实部，奇数点作为虚部）
	for(l=1;l<=m;l++)
		{
			b=(int)(pow(2.0,(l-1))+0.5);					//计算出2的L-1次幂送给b
				
//第二层循环					
			for(j=0;j<b;j++)
				{
					p=j*((int)(pow(2.0,(m-l))+0.5));		//计算出p等于j乘以2的m-l次幂
					
//第三层循环
					for(k=j;k<n;k+=(int)(pow(2.0,l)+0.5))
						{
							rpartkb=(*(leftdata_p+(k+b)))*(*(rcgene_p+p))-(*(ibufferfft_p+(k+b)))*(*(icgene_p+p));
							ipartkb=(*(ibufferfft_p+(k+b)))*(*(rcgene_p+p))+(*(leftdata_p+(k+b)))*(*(icgene_p+p));
							*(leftdata_p+(k+b))=*(leftdata_p+k)-rpartkb;
							*(ibufferfft_p+(k+b))=*(ibufferfft_p+k)-ipartkb;
							*(leftdata_p+k)=*(leftdata_p+k)+rpartkb;
							*(ibufferfft_p+k)=*(ibufferfft_p+k)+ipartkb;
						}
				}
		}

}
						


//倒序子程序    			
changeorder(float COFFTDATA[],float COFFTIMAGE[],int n) 
{   					
    int lh,i,j,n1,k;
    float t;
    lh=n/2;
    j=lh;
    n1=n-1;
    for(i=1;i<=n1;i++)
    	{
    		if(i<j)
    			{
    				t=COFFTDATA[i];
    				COFFTDATA[i]=COFFTDATA[j];
    				COFFTDATA[j]=t;
    				t=COFFTIMAGE[i];
    				COFFTIMAGE[i]=COFFTIMAGE[j];
    				COFFTIMAGE[j]=t;
    			}
    		k=lh;
			while(j>=k)
				{
					j=j-k;
					k=(int)(k/2+0.5);
					if((j==0)&&(k==0))break;
				}
			j=j+k;
		}
}