基于整型运算的FFT计算程序

来源:本站
导读:目前正在解读《基于整型运算的FFT计算程序》的相关信息,《基于整型运算的FFT计算程序》是由用户自行发布的知识型内容!下面请观看由(电工技术网 - www.9ddd.net)用户发布《基于整型运算的FFT计算程序》的详细说明。
简介:FFT计算比较费时,这是由于计算过程中使用浮点数以及需要大量计算sin、cos函数。

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不适用

提醒:《基于整型运算的FFT计算程序》最后刷新时间 2024-03-14 01:02:51,本站为公益型个人网站,仅供个人学习和记录信息,不进行任何商业性质的盈利。如果内容、图片资源失效或内容涉及侵权,请反馈至,我们会及时处理。本站只保证内容的可读性,无法保证真实性,《基于整型运算的FFT计算程序》该内容的真实性请自行鉴别。