[原创]基于整型运算的FFT计算程序【恢复】
FFT计算比较费时，这是由于计算过程中使用浮点数以及需要大量计算sin、cos函数。常规方法实现FFT的C代码如下（参见数值计算与信号处理，输入为实数序列）：#include "math.h"
void rfftd(double *x, int n)
{
int i, j, k, m, i1, i2, i3, i4, n1, n2, n4;
double a, e, cc, ss, xt, t1, t2;
for(j = 1, i = 1; i < 16; i ++)
{
m = i;
j = 2 * j;
if(j == n)
break;
}
n1 = n - 1;
for(j = 0, i = 0; i < n1; i ++)
{
if(i < j)
{
xt = x;
x = x;
x = xt;
}
k = n / 2;
while(k < (j + 1))
{
j -= k;
k /= 2;
}
j += k;
}
for(i = 0; i < n; i += 2)
{
xt = x;
x += x;
x = xt - x;
}
n2 = 1;
for(k = 2; k <= m; k ++)
{
n4 = n2;
n2 = 2 * n4;
n1 = 2 * n2;
e = 6.28318530718 / n1;
for(i = 0; i < n; i += n1)
{
xt = x;
x += x;
x = xt - x;
x = -x;
a = e;
for(j = 1; j <= (n4 - 1); j ++)
{
i1 = i + j;
i2 = i - j + n2;
i3 = i + j + n2;
i4 = i - j + n1;
cc = cos(a);
ss = sin(a);
a += e;
t1 = cc * x + ss * x;
t2 = ss * x - cc * x;
x = x - t2;
x = -x - t2;
x = x - t1;
x = x +t1;
}
}
}
}
参数x为要变换的数据的指针，n为数据的个数（必须为2的整数次幂），变换后的结果从x中输出（只存放前n/2+1个值），存储顺序为（Re和Im分别为实部和虚部）。
把这段代码转换成整型运算，可用的方法是：1、把所有浮点数乘以2的N次幂，例如256，取整成整型数；2、sin和cos函数采用查表法实现。修改后的整型FFT运算代码如下：
long SIN_TABLE256 = {0, 4, 9, 13, 18, 22, 27, 31, 36, 40,
44, 49, 53, 58, 62, 66, 71, 75, 79, 83,
88, 92, 96, 100, 104, 108, 112, 116, 120, 124,
128, 132, 136, 139, 143, 147, 150, 154, 158, 161,
165, 168, 171, 175, 178, 181, 184, 187, 190, 193,
196, 199, 202, 204, 207, 210, 212, 215, 217, 219,
222, 224, 226, 228, 230, 232, 234, 236, 237, 239,
241, 242, 243, 245, 246, 247, 248, 249, 250, 251,
252, 253, 254, 254, 255, 255, 255, 256, 256, 256,
256};
long fastsin256(long degree)
{
long result, neg = 0;
if(degree < 0)
degree = -degree + 180;
if(degree>= 360)
degree %= 360;
if(degree>= 180)
{
degree -= 180;
neg = 1;
}
if((degree>= 0) && (degree <= 90))
result = SIN_TABLE256;
else
result = SIN_TABLE256;
if(!neg)
return result;
else
return -result;
}
__inline long fastcos256(long degree)
{
return fastsin256(degree + 90);
}
void rfftl(long *x, int n)
{
int i, j, k, m, i1, i2, i3, i4, n1, n2, n4;
long a, e, cc, ss, xt, t1, t2;
for(j = 1, i = 1; i < 16; i ++)
{
m = i;
j = (j << 1);
if(j == n)
break;
}
n1 = n - 1;
for(j = 0, i = 0; i < n1; i ++)
{
if(i < j)
{
xt = x;
x = x;
x = xt;
}
k = (n>> 1);
while(k < (j + 1))
{
j -= k;
k = (k>> 1);
}
j += k;
}
for(i = 0; i < n; i += 2)
{
xt = x;
x += x;
x = xt - x;
}
n2 = 1;
for(k = 2; k <= m; k ++)
{
n4 = n2;
n2 = (n4 << 1);
n1 = (n2 << 1);
e = 360 / n1;
for(i = 0; i < n; i += n1)
{
xt = x;
x += x;
x = xt - x;
x = -x;
a = e;
for(j = 1; j <= (n4 - 1); j ++)
{
i1 = i + j;
i2 = i - j + n2;
i3 = i + j + n2;
i4 = i - j + n1;
cc = fastcos256(a);
ss = fastsin256(a);
a += e;
t1 = cc * x + ss * x;
t2 = ss * x - cc * x;
t1 = (t1>> 8);
t2 = (t2>> 8);
x = x - t2;
x = -x - t2;
x = x - t1;
x = x + t1;
}
}
}
}
运算过程中所有浮点数都乘以256并取整，输入和输出数据也是乘以256之后的整型数据。sin和cos采用查表实现，精确到1度。程序适用于对精度要求不高的FFT计算，例如音频播放器的频谱显示等。采用更大的取整系数（例如65536）并增加sin、cos表的精度可以提高这个整型FFT计算的精度。
应用转化成整型计算的FFT需要注意的问题是，整型FFT的动态范围会比较小，这是由定点数的性质决定的，因此如果计算对在很大的动态范围内的精度有要求，则整型FFT不适用。 留个脚印 也顶一下，日后有用 顶一下..虽然不知道怎么用..
本贴被 hyz_avr 编辑过,最后修改时间：2008-12-28,17:18:12. 1、把所有浮点数乘以2的N次幂，例如256，取整成整型数；2、sin和cos函数采用查表法实现。
MARK 好的东西，就得支持一下。 MARK 老帖也要mark 看了半天，了解了
去试验下的说
看看是不是真的管用哈~~~ 晕~~~51的DATA不够~~~
想办法缩量ing 看来这东西51说绝对放不下了的说~~~
唉~~~回家跑2148上看看 MARK MARK MARK MARK下，太有用了～ MARK mark 那个正弦余弦表做得不够好，
不应该用角度做单位，应该归一化成256。 顶，已经收集了不少FFT方面的东西，抓紧学习啦，谢谢 呵呵，慎用整型算法。精度不够，出来结果不好说。 mark mark所有FFT FFT 这算出来的结果不对吧 收藏了~~~ Mark 怎么都是乱码？
页:
[1]