FFT计算比较费时,这是由于计算过程中使用浮点数以及需要大量计算sin、cos函数。常规方法实现FFT的C代码如下(参见数值计算与信号处理,输入为实数序列):
#i nclude "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[i]; x[i] = xt; } k = n / 2; while(k < (j + 1)) { j -= k; k /= 2; } j += k; } for(i = 0; i < n; i += 2) { xt = x[i]; x[i] += 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[i]; x[i] += 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[i]; x[i] = xt; } k = (n >> 1); while(k < (j + 1)) { j -= k; k = (k >> 1); } j += k; } for(i = 0; i < n; i += 2) { xt = x[i]; x[i] += 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[i]; x[i] += 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不适用