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)
{
double real[256];
double imag[256];
long swap;
long offset;
double tmp_real, tmp_imag, temp;// 临时变量, 记录实部
double tw1, tw2; // 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.
long cur_layer, gr_num, i, k, p;
long step; // 步进
long sample_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<<i;
//每个颗粒的输入样本数N
sample_num = 1<<cur_layer;
//步进. 步进是N/2
step = sample_num/2;
k = 0;
//对颗粒进行循环
for(i=0; i<gr_num; i++)
{
//对样本点进行循环, 注意上限和步进
for(p=0; p<sample_num/2; 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;
}
}
