|
|
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[j];
x[j] = 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[i + 1];
x[i + 1] = xt - x[i + 1];
}
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[i + n2];
x[i + n2] = xt - x[i + n2];
x[i + n2 + n4] = -x[i + n2 + n4];
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[i3] + ss * x[i4];
t2 = ss * x[i3] - cc * x[i4];
x[i4] = x[i2] - t2;
x[i3] = -x[i2] - t2;
x[i2] = x[i1] - t1;
x[i1] = x[i1] +t1;
}
}
}
}
参数x为要变换的数据的指针，n为数据的个数（必须为2的整数次幂），变换后的结果从x中输出（只存放前n/2+1个值），存储顺序为[Re0, Re1, …, Ren/2, Imn/2-1, …, Im1]（Re和Im分别为实部和虚部）。
把这段代码转换成整型运算，可用的方法是：1、把所有浮点数乘以2的N次幂，例如256，取整成整型数；2、sin和cos函数采用查表法实现。修改后的整型FFT运算代码如下：
long SIN_TABLE256[91] = {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[degree];
else
result = SIN_TABLE256[180 - degree];
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[j];
x[j] = 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[i + 1];
x[i + 1] = xt - x[i + 1];
}
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[i + n2];
x[i + n2] = xt - x[i + n2];
x[i + n2 + n4] = -x[i + n2 + n4];
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[i3] + ss * x[i4];
t2 = ss * x[i3] - cc * x[i4];
t1 = (t1>> 8);
t2 = (t2>> 8);
x[i4] = x[i2] - t2;
x[i3] = -x[i2] - t2;
x[i2] = x[i1] - t1;
x[i1] = x[i1] + t1;
}
}
}
}
运算过程中所有浮点数都乘以256并取整，输入和输出数据也是乘以256之后的整型数据。sin和cos采用查表实现，精确到1度。程序适用于对精度要求不高的FFT计算，例如音频播放器的频谱显示等。采用更大的取整系数（例如65536）并增加sin、cos表的精度可以提高这个整型FFT计算的精度。
应用转化成整型计算的FFT需要注意的问题是，整型FFT的动态范围会比较小，这是由定点数的性质决定的，因此如果计算对在很大的动态范围内的精度有要求，则整型FFT不适用。 |
阿莫论坛20周年了！感谢大家的支持与爱护！！
月入3000的是反美的。收入3万是亲美的。收入30万是移民美国的。收入300万是取得绿卡后回国，教唆那些3000来反美的！
|
发表于 2008-12-28 16:30:48
