fft.c

快速傅里叶变换的C实现

#include	<stdio.h>
#include 	<math.h>
struct _complex
{
	double x;
	double y;	
};
#define N 32
struct _complex F[N];
unsigned int f[N]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31};
unsigned int r;
void conver_bit(void);
int main(void)
{
    struct _complex W[N/2];
    double angle; 
    int DFTn;/*第DFTn级*/
    int k;/*分级有多少个蝶形运算*/
    int d;/*蝶形运算的偏移*/
    int p;/*index */
    int v;
    struct _complex X1,X2;
	unsigned int c;
	unsigned int j;
	/*需要多少级蝶形运算*/
    r = log(N +1) / log(2);	
	/*把时域的转换到频域*/
	for(c = 0; c < N ; c++)
	{
		F[c].x = f[c];
		F[c].y = 0;
	}
	/*倒序*/
	conver_bit();
	
	
	/*计算旋转因子W, 加权系数 */

   
    for( j = 0; j < N / 2; j++) 
    { 
        angle = - j * M_PI * 2 / N; 
        W[j].x = cos(angle);
        W[j].y = sin(angle);
    } 

	/*采用蝶形算法进行快速傅立叶变换*/

    for( DFTn = 0; DFTn < r; DFTn++)
    {/*第几级蝶形运算*/
        for( k= 0; k < (1 << (r - DFTn - 1)); ++k)
        {/*当前列所有的蝶形进行运算*/
            p = 2 * k *  (1 << DFTn);
            for ( d = 0; d < (1 << DFTn); ++d)
            {/*相邻蝶形运算*/
                X1 = F[p+d];
                X2 = F[p+d+(1 << DFTn)];
                if(DFTn==0)/*第0级蝶形运算,第0级蝶形运算有点特别,乘法中没有虚数相乘*/
                {
              	  F[p+d].x = X1.x + X2.x * W[d*(1<<(r - DFTn - 1))].x;
 					  F[p+d].y = X1.y + X2.x * W[d*(1<<(r - DFTn - 1))].y;
					  F[p+d+(1 << DFTn)].x = X1.x - X2.x * W[d*(1<<(r - DFTn - 1))].x;
					  F[p+d+(1 << DFTn)].y = X1.y - X2.x * W[d*(1<<(r - DFTn - 1))].y;
           	}
                else/*乘法中有虚数相乘*/
                {
                   F[p+d].x = X1.x + (X2.x * W[d*(1<<(r - DFTn - 1))].x - X2.y * W[d*(1<<(r - DFTn - 1))].y);
                   F[p+d].y = X1.y + (X2.y * W[d*(1<<(r - DFTn - 1))].x + X2.x * W[d*(1<<(r - DFTn - 1))].y);
                   F[p+d+(1 << DFTn)].x = X1.x - (X2.x * W[d*(1<<(r - DFTn - 1))].x - X2.y * W[d*(1<<(r - DFTn - 1))].y);
				   F[p+d+(1 << DFTn)].y = X1.y - (X2.y * W[d*(1<<(r - DFTn - 1))].x + X2.x * W[d*(1<<(r - DFTn - 1))].y);
                }
            }
        }
    }
    for(v = 0; v < N; v++)
    {
    	F[v].x /= 100;
    	F[v].y /= 100;
    }
    for(v=0;v<N;v++)
    	printf("F[%d] = (%0.4f+(%0.4fi))*100\n",v,F[v].x,F[v].y);
}

/*DIT—FFT序列的倒序算法*/
void conver_bit(void)
{
	int j[N];/*倒序后存放在这个数组里*/
	int k,q,v,g,n,d,e,f,i;
	struct _complex dTemp;
	memset(j,0,sizeof(j));
	j[0]=0;
	j[1]=N/2;
	for(k = 1; k <= r-1; k++)
	{
		f = pow(2,(k-1));
		for(q = 1; q <= f; q++)
		{
			e = ((q-1)*N)/pow(2,(k-1))+1;
			d = ((q-1)*N)/pow(2,(k-1))+2;
			for(v = e; v <= d; v++)
			{
				n = v+N/pow(2,k)-1;/*原先的算法里没有－１，因为作者的数组索引是从１开始的，这里是从０开始的*/
				g = j[v-1] + pow(2,k-1);
				j[n] = g;
			}				
		}
	}
	for(g = 0; g < N; g++)
	{
		i = j[g];
		if(i > g)
		{
			dTemp = F[g];
			F[g] = F[i];
			F[i] = dTemp;
		}
		/*printf("F[%d] = %0.4f\n",g,F[g].x);*/
		/*printf("F[%d] = F[%d]\n",g,i);*/
	}
	/*for(g = 0; g < N ; g++)
		printf("F[%d] = %0.4f\n",g,F[g].x);*/
}