整型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[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不适用。