c语言的256个字符,C语言版 256点FFT算法

从网上找到的一段感觉比较好的FFT函数，可以直接使用，出处由于时间原因已经忘记了@@

//文件名：fft.c

#include "fft.h"

//128等分余弦表

static const double cosTable[] =

{

1.0000,0.9997,0.9988,0.9973,0.9952,0.9925,0.9892,0.9853,

0.9808,0.9757,0.9700,0.9638,0.9569,0.9495,0.9415,0.9330,

0.9239,0.9142,0.9040,0.8932,0.8819,0.8701,0.8577,0.8449,

0.8315,0.8176,0.8032,0.7883,0.7730,0.7572,0.7410,0.7242,

0.7071,0.6895,0.6716,0.6532,0.6344,0.6152,0.5957,0.5758,

0.5556,0.5350,0.5141,0.4929,0.4714,0.4496,0.4276,0.4052,

0.3827,0.3599,0.3369,0.3137,0.2903,0.2667,0.2430,0.2191,

0.1951,0.1710,0.1467,0.1224,0.0980,0.0736,0.0491,0.0245,

-0.0000,-0.0245,-0.0491,-0.0736,-0.0980,-0.1224,-0.1467,-0.1710,

-0.1951,-0.2191,-0.2430,-0.2667,-0.2903,-0.3137,-0.3369,-0.3599,

-0.3827,-0.4052,-0.4276,-0.4496,-0.4714,-0.4929,-0.5141,-0.5350,

-0.5556,-0.5758,-0.5957,-0.6152,-0.6344,-0.6532,-0.6716,-0.6895,

-0.7071,-0.7242,-0.7410,-0.7572,-0.7730,-0.7883,-0.8032,-0.8176,

-0.8315,-0.8449,-0.8577,-0.8701,-0.8819,-0.8932,-0.9040,-0.9142,

-0.9239,-0.9330,-0.9415,-0.9495,-0.9569,-0.9638,-0.9700,-0.9757,

-0.9808,-0.9853,-0.9892,-0.9925,-0.9952,-0.9973,-0.9988,-0.9997,

};

//128等分正弦表

static const double sinTable[] =

{

0.0000,0.0245,0.0491,0.0736,0.0980,0.1224,0.1467,0.1710,

0.1951,0.2191,0.2430,0.2667,0.2903,0.3137,0.3369,0.3599,

0.3827,0.4052,0.4276,0.4496,0.4714,0.4929,0.5141,0.5350,

0.5556,0.5758,0.5957,0.6152,0.6344,0.6532,0.6716,0.6895,

0.7071,0.7242,0.7410,0.7572,0.7730,0.7883,0.8032,0.8176,

0.8315,0.8449,0.8577,0.8701,0.8819,0.8932,0.9040,0.9142,

0.9239,0.9330,0.9415,0.9495,0.9569,0.9638,0.9700,0.9757,

0.9808,0.9853,0.9892,0.9925,0.9952,0.9973,0.9988,0.9997,

1.0000,0.9997,0.9988,0.9973,0.9952,0.9925,0.9892,0.9853,

0.9808,0.9757,0.9700,0.9638,0.9569,0.9495,0.9415,0.9330,

0.9239,0.9142,0.9040,0.8932,0.8819,0.8701,0.8577,0.8449,

0.8315,0.8176,0.8032,0.7883,0.7730,0.7572,0.7410,0.7242,

0.7071,0.6895,0.6716,0.6532,0.6344,0.6152,0.5957,0.5758,

0.5556,0.5350,0.5141,0.4929,0.4714,0.4496,0.4276,0.4052,

0.3827,0.3599,0.3369,0.3137,0.2903,0.2667,0.2430,0.2191,

0.1951,0.1710,0.1467,0.1224,0.0980,0.0736,0.0491,0.0245,

};

//快速求平方根

float sqrt_fast(float x)

{

long lx;

float temp;

float xhalf = x * 0.5f;

lx = *((long*)&x);

lx = 0x5F375A86 - (lx>>1);

temp = *((float*)&lx);

temp = temp*(1.5f - xhalf*temp*temp);

temp = temp*(1.5f - xhalf*temp*temp);

temp = temp*(1.5f - xhalf*temp*temp);

return temp * x;

}

//8bit位反转

unsigned char bitrev(unsigned char i)

{

unsigned char r=0;

r |= (i>>0)&0x01;

r<<= 1;

r |= (i>>1)&0x01;

r<<= 1;

r |= (i>>2)&0x01;

r<<= 1;

r |= (i>>3)&0x01;

r<<= 1;

r |= (i>>4)&0x01;

r<<= 1;

r |= (i>>5)&0x01;

r<<= 1;

r |= (i>>6)&0x01;

r<<= 1;

r |= (i>>7)&0x01;

return r;

}

//256点 fft运算

//输入256个采样数据,只用实数部分，虚数为0

//输出128个频域数据,已将复数转换为幅度(频谱)

void FFT256(long* pInputData, long* pOutputData)

{

doublereal[256];

doubleimag[256];

longswap;

longoffset;

doubletmp_real, tmp_imag, temp;// 临时变量, 记录实部

doubletw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.

longcur_layer, gr_num, i, k, p;

longstep;// 步进

longsample_num;// 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)

//虚部为0

for(i=0; i<256; i++)

{

imag[i] = 0;

real[i] = 0;

}

//倒位序

for(i=0; i<256; i++)

{

swap = bitrev((unsigned char)i);

real[swap] = (float)pInputData[i];

}

// 对层循环

for(cur_layer=1; cur_layer<=8; cur_layer++)//2^8=256,所以循环8次

{

// 求当前层拥有多少个颗粒(gr_num)

i = 8 - cur_layer;

gr_num=1<

//每个颗粒的输入样本数N

sample_num = 1<

//步进. 步进是N/2

step = sample_num/2;

k = 0;

//对颗粒进行循环

for(i=0; i

{

//对样本点进行循环, 注意上限和步进

for(p=0; p

{

offset = p*gr_num;

tw1 = cosTable[offset];

tw2 = sinTable[offset];

tmp_real = real[k+p];

tmp_imag = imag[k+p];

temp = real[k+p+step];

//蝶形算法

real[k+p] = tmp_real + ( tw1*real[k+p+step] - tw2*imag[k+p+step] );

imag[k+p] = tmp_imag + ( tw2*real[k+p+step] + tw1*imag[k+p+step] );

real[k+p+step] = tmp_real - ( tw1*temp - tw2*imag[k+p+step] );

imag[k+p+step] = tmp_imag - ( tw2*temp + tw1*imag[k+p+step] );

}

//开跳

k += 2*step;

}

}

//对复数求平方根

for(i=0; i<128;i++)

{

temp = real[i]*real[i] + imag[i]*imag[i];

temp = sqrt_fast(temp);

pOutputData[i] = (long)temp;

}

}