论坛首页 › DSP › FFT算法C代码

zouzhichao 发表于 2013-8-11 23:26:38

FFT算法C代码

本帖最后由 zouzhichao 于 2013-8-11 23:30 编辑

基本蝶形运算单元：


注：上图中的运算均为复数运算。
复数表示形式：
Z = R +V * i
加法运算：
    Z1 + Z2 = (R1 + V1 * i) + (R1 + V1 * i) = (R1 + R2) + (V1 + V2) * i
减法运算：
    Z1 - Z2 = (R1 + V1 * i) - (R1 + V1 * i) = (R1 - R2) + (V1 - V2) * i
乘法运算：
Z1 * Z2 = (R1 + V1 * i) * (R2 + V2 * i)
       = R1 * R2 + V1 * V2 * i * i + R1 * V2 * i + R2 * V1 * i
       = (R1 * R2 - V1 * V2) + (R1 * V2 + R2 * V1) * i
除法运算：
Z1 / Z2 = (R1 + V1 * i) / (R2 + V2 * i)
      = ((R1 + V1 * i) * (R2 – V2 * i)) / ( (R2 + V2 * i) * (R2 – V2 * i))
      = (R1 * R2 + V1 * V2) / (R2 * R2 + V2 * V2) +
((R2 * V1 - R1 * V2) /(R2 * R2 + V2 * V2) * i
C 代码实现：
void ButterFly(COMPLEX* SrcDstA, COMPLEX* SrcDstB, const COMPLEX* Modulus)
{
        COMPLEX Temp;
        //Temp = C * B
        Temp.Real = (Modulus->Real) * (SrcDstB->Real) - (Modulus->Virtual) * (SrcDstB->Virtual);
        Temp.Virtual = (Modulus->Real) * (SrcDstB->Virtual) + (Modulus->Virtual) * (SrcDstB->Real);
//B = A - C * B = A - Temp
        SrcDstB->Real = SrcDstA->Real - Temp.Real;
        SrcDstB->Virtual = SrcDstA->Virtual - Temp.Virtual;
//A = A + C * B = A + Temp
        SrcDstA->Real += Temp.Real;
        SrcDstA->Virtual += Temp.Virtual;
}

8点DIT-FFT的蝶形运算流图：



由上图可知一个8点DIT-FFT包含3级的蝶形运算，同理16点DIT-FFT包含4级的蝶形运算，32点DIT-FFT包含5级的蝶形运算。
对于8点DIT-FFT，输入序列的顺序为x0、x4、x2、x6、x1、x5、x3、x7, 输出序列的顺序为X0、X4、X2、X6、X1、X5、X3、X7，因此需要对输入序列进行倒序。
对于8点DIT-FFT，总共3级蝶形运算，即第0级，第1级，第2级.
第0级：分为4组，每组1个运算单元，（0,1）（2,3）（4,4）（6,7），W系数都为W80；
第1级：分为2组，每组2个运算单元，（0,2）（1,3）（4,6）（5,7），W系数分别为W80，W82，W80，W82；
第2级：分为1组，每组4个运算单元，（0,4）（1,5）（2,6）（3,7），W系数分别为W80，W81，W82，W83；
对于N = 2^K点，第i级：分为2^(K-i-1) 组，每组2^i个运算单元。
C代码实现
void FftK(COMPLEX* SrcDst, const COMPLEX* Modulus, uint16_ N, uint8_ K)
{
        uint16_ i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
        for (i = N >> (K + 1); i; i--)
                for (j = J; j > 0; j--)
                        ButterFly(SrcDst + (i - 1) * X + j - 1, SrcDst + i * X + j - J - 1, Modulus + (j - 1) * Y);
}

void FFT(COMPLEX* SrcDst, const COMPLEX* Modulus, uint16_ N)
{
        uint16_ k, K = N2K(N);
        for (k = 0; k < K; k++)
                FftK(SrcDst, Modulus, N, k);
}

W系数：
利用欧拉公式：
Wk = e^(2 * k * Pi * i / N)
    = cos(2 * k * Pi / N) + sin(2 * k * Pi / N)

void InitW(COMPLEX* Modulus, uint16_ N)
{
        double Pi = 3.1415926535 * 2 / N, Temp = 0;
        for (N >>= 1; N; N--, Modulus++, Temp += Pi)
        {
                Modulus->Real = cos(Temp);
                Modulus->Virtual = sin(Temp);
        }
}

DIT-FFT运算的倒序规律：


由上图可知，倒序就是将二进制数倒序即可，例如256点DIT-FFT，35 = 00100011，倒序之后为11000100 = 196，即x35和x196交换即可。
对于8点DIT-FFT，输入有8个数据，位宽为3，根据对称的原则，那么会有2^((3 + 1) / 2) = 4个数与自身对称，不需要倒序，如：
000 = 0，101 = 5，010 = 2，111 = 7
因此可知对于N = 2^K点DIT-FFT，只有2^K - 2^((K + 1) / 2)个数据需要交换，因此需要交换的次数为(2^K - 2^((K + 1) / 2) )/ 2。
倒序交换可以用硬件很容易实现，软件的话需要构造一个倒序的“加法”运算。以8点DIT-FFT为例分析：
（1）        第0个数为000
（2）        第1个数在000基础上“加1”，即最高位加1,变成100 = 4，将x1与x4交换。
（3）        第2个数在100基础上“加1”，即最高位加1,溢出，往次高位进1，变成010 = 2，将x2与x2交换，即不用交换。
（4）        第3个数在010基础上“加1”，即最高位加1,变成110 = 6，将x3与x6交换。
（5）        第4个数在110基础上“加1”，即最高位加1, 溢出，往次高位进1，溢出，再往次次高位进1，变成001= 6，将x4与x1交换，由于前面已经将x1与x4交换，故此步骤省略。
（6）        第5个数在001基础上“加1”，即最高位加1,变成101 = 5，将x5与x5交换，即不用交换。
（7）        第6个数在101基础上“加1”，即最高位加1, 溢出，往次高位进1，变成011 = 3，将x6与x3交换，由于前面已经将x3与x6交换，故此步骤省略。
（8）        第7个数在101基础上“加1”，即最高位加1, 溢出，往次高位进1，变成111 = 7，将x7与x7交换，即不用交换。
换一种思路：
（1）        0
（2）        0 < 4，0 + 4 = 4
（3）        4 >= 4，4 - 4 = 0，0 + 2 = 2
（4）        2 < 4，2 + 4 = 6
（5）        6 >= 4，6 - 4 = 2，2 >= 2，2 – 2 = 0，0 + 1 = 1
（6）        1 < 4，1 + 4 = 5
（7）        5 >= 4，5 - 1 = 1，1 + 2 = 3
（8）        3 < 4，3 + 4 = 7
C代码实现：
void Change(COMPLEX* SrcDst, uint16_ N)
{
        uint16_ v, i, j, k = N >> 1;                       
        COMPLEX Temp;
        for (i = 0, j = 1; j < N; j++)
        {
                for (v = k; i >= v; i-= v, v >>= 1)
                        ;
                i += v;
                if (j < i)
                {
                        Temp = *(SrcDst + i);
                        *(SrcDst + i) = *(SrcDst + j);
                        *(SrcDst + j) = Temp;
                }
        }
}

基本代码如下：

typedef signed char int8_;
typedef unsigned char uint8_;
typedef signed short int16_;
typedef unsigned short uint16_;
typedef signed long int32_;
typedef float flt_;
typedef double lflt_;

typedef struct tagCOMPLEX
{
        /* data */
        lflt_ Real;
        lflt_ Virtual;
}COMPLEX;

void ButterFly(COMPLEX* SrcDstA, COMPLEX* SrcDstB, const COMPLEX* Modulus)
{
        COMPLEX Temp;
        Temp.Real = (Modulus->Real) * (SrcDstB->Real) - (Modulus->Virtual) * (SrcDstB->Virtual);
        Temp.Virtual = (Modulus->Real) * (SrcDstB->Virtual) + (Modulus->Virtual) * (SrcDstB->Real);
        SrcDstB->Real = SrcDstA->Real - Temp.Real;
        SrcDstB->Virtual = SrcDstA->Virtual - Temp.Virtual;
        SrcDstA->Real += Temp.Real;
        SrcDstA->Virtual += Temp.Virtual;
}

void InitW(COMPLEX* Modulus, uint16_ N)
{
        double Pi = 3.1415926535 * 2 / N, Temp = 0;
        for (N >>= 1; N; N--, Modulus++, Temp += Pi)
        {
                Modulus->Real = cos(Temp);
                Modulus->Virtual = sin(Temp);
        }
}

void Change(COMPLEX* SrcDst, uint16_ N)
{
        uint16_ v, i, j, k = N >> 1;                       
        COMPLEX Temp;
        for (i = 0, j = 1; j < N; j++)
        {
                for (v = k; i >= v; i-= v, v >>= 1)
                        ;
                i += v;
                if (j < i)
                {
                        Temp = *(SrcDst + i);
                        *(SrcDst + i) = *(SrcDst + j);
                        *(SrcDst + j) = Temp;
                }
        }
}

uint8_ N2K(uint16_ N)
{
        uint8_ k;
        for (k = 0, N--; N ; N >>= 1, k++)
                ;
        return k;
}

void FftK(COMPLEX* SrcDst, const COMPLEX* Modulus, uint16_ N, uint8_ K)
{
        uint16_ i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
        for (i = N >> (K + 1); i; i--)
                for (j = J; j > 0; j--)
                        ButterFly(SrcDst + (i - 1) * X + j - 1, SrcDst + i * X + j - J - 1, Modulus + (j - 1) * Y);
}

void FFT(COMPLEX* SrcDst, const COMPLEX* Modulus, uint16_ N)
{
        uint16_ k, K = N2K(N);
        for (k = 0; k < K; k++)
                FftK(SrcDst, Modulus, N, k);
}

FFT算法优化
注：这阶段的优化都还是基于浮点型的，暂时还没有把浮点型的运算改到定点运算。
W系数获取的优化
原C函数：
void InitW(COMPLEX* Modulus, uint16_ N)
{
        double Pi = 3.1415926535 * 2 / N, Temp = 0;
        for (N >>= 1; N; N--, Modulus++, Temp += Pi)
        {
                Modulus->Real = cos(Temp);
                Modulus->Virtual = sin(Temp);
        }
}
对于一个256点的DIT-FFT，需要128个复数系数，即 256 个cos/ sin运算，这个有不小的优化空间。
sin(2 * n * a) = 2 * sin(n * a) * cos(n * a)
cos(2 * n * a) = 2 * cos(n * a) * cos(n * a) – 1
cos(2 * n * a) = 1 - 2 * sin(n * a) * sin(n * a)
sin((2 * n + 1) * a) = sin(a) * cos(2 * n * a) + cos(a) * sin(2 * n * a)
cos((2 * n + 1) * a) = cos(a) * cos(2 * n * a) - sin(a) * sin(2 * n * a)
优化后只需要1 个cos/ sin运算，剩下的只是乘法运算。
优化后的C函数：
void InitW1(COMPLEX* Modulus, uint16_ N)
{
        uint16_ i, j;
        double Pi = 3.141592653589793 * 2.0 / N;
        Modulus.Real = 1;
        Modulus.Virtual = 0;
        Modulus.Real = cos(Pi);
        Modulus.Virtual = sin(Pi);
        for (i = 1, j = 2, N >>= 2; i < N; i++, j += 2)
        {
                Modulus.Real = 2.0 * Modulus.Real * Modulus.Real - 1;
                Modulus.Virtual = 2.0 * Modulus.Real * Modulus.Virtual;
                Modulus.Real = Modulus.Real * Modulus.Real - Modulus.Virtual * Modulus.Virtual;
                Modulus.Virtual = Modulus.Real * Modulus.Virtual + Modulus.Virtual * Modulus.Real;
        }
}

以空间换时间
将内存地址A的1024个数复制到内存地址B处：
（1）
for (i = 0; i < 1024; i++)
        B = A;
（2）
B = A;
B = A;
B = A;
B = A;
……
……
B = A;
B = A;
显然，（1）代码长度小于（2），（1）运算时间大于（2），即（1）为代码大小优先，（2）为运算速度优先。
如果你有听说过达夫设备那就很容易理解了。
在FFT运算中通常都是需要运算速度优先，对于一个确定点数的DIT-FFT，它的很多“参数”都是编程时就可以确定的，不需要程序一边运行以便计算，这里的“参数”不单单只是变量，还有运算顺序。

（1）        对于InitW函数，由于在多个运算中，它也只运行一次，因此它速度优先的需求并不高。
（2）        Change函数每次DIT-FFT运算都会要运行一次，因此它有很大的改进空间。
（3）        同理对于FftK，FFT函数，也需要进行优化。
以64点DIT-FFT为例：
（1）Change需要进行(2^6 - 2^((6 + 1) / 2) ) / 2 = (64 - 8 ) / 2 = 28次交换运算。
void Change(Data)
{
        #define _Change(A,B) {Temp = A; A = B; B = Temp;}
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
        _Change(Data, Data);
}
（2）对于FftK，第0级运算有特殊性，可以特殊处理。
（3）对于整个FftK，FFT，运算顺序，运算次数都是固定的，因此可以取消循环，做成顺序执行的程序，如Change函数一样。对于这点会增加很多运算速度，同事也占用很长的代码，同时需要程序员做写更多的工作，这部分可以参考单片机小精灵自己做一个代码生成器。

zouzhichao 发表于 2013-8-11 23:31:56

工程文件还没整理好，迟点上传
还写了一个代码生成工具，也迟点整理再上传

cos 发表于 2013-8-12 23:20:55

好东东,期待中.

zouzhichao 发表于 2013-8-12 23:23:36

cos 发表于 2013-8-12 23:20 static/image/common/back.gif
好东东,期待中.

今天还上传不了，没带优盘

leifeng 发表于 2013-8-12 23:47:54

这个真的很历害!楼主NB!

zouzhichao 发表于 2013-8-13 00:03:24

leifeng 发表于 2013-8-12 23:47 static/image/common/back.gif
这个真的很历害!楼主NB!

承蒙谬赞，愧不敢当

zouzhichao 发表于 2013-8-13 00:06:48

上传dit-fft完整源代码，在TCC，VC6.0，GCC测试过：

VC6.0工程：

zouzhichao 发表于 2013-8-13 00:09:10

更快的DIT-FFT代码需要用DIT-FFT代码生成器生成，这个小软件很快会上传的

zouzhichao 发表于 2013-8-13 02:32:34

代码生成器生成的头文件：

测试模板：
#include "stdio.h"
#include "dit-fft-cal-1024.h"

typedef signed char int8_;
typedef unsigned char uint8_;
typedef signed short int16_;
typedef unsigned short uint16_;
typedef signed long int32_;
typedef unsigned long uint32_;
typedef float flt_;
typedef double lflt_;

int8_ GetData(const char* FileName, lflt_* Buffer, uint32_ N)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for ( ; N; N--, Buffer++)
                if (!fscanf(File, "%lf", Buffer))
                        break;
        fclose(File);
        return N ? -1 : 0;
}

lflt_ GetAverage(const lflt_* Buffer, uint32_ N)
{
        uint32_ i;
        lflt_ Temp = 0.0;
        if (!N)
                return Temp;
        for (i = N; i; i--)
                Temp += *Buffer++;
        return Temp / N;
}

int8_ SubAverage(lflt_* Buffer, uint32_ N)
{
        lflt_ Temp = GetAverage((const lflt_*)Buffer, N);
        for ( ; N; N--)
                *Buffer++ -= Temp;
        return 0;
}

int main(int argc, char* argv[])
{
        #define N 1024
        uint32_ i;
        lflt_ Temp1, Temp2, BufferR, BufferV = {0.0};
        GetData((const char*)"data.txt", (lflt_*)BufferR, N);
        SubAverage((lflt_*)BufferR, N);
//
        _ReOrder(BufferR, Temp1);//头文件里有定义，Temp1为临时变量，在宏调用前定义
        _Butters(BufferR, BufferV, Temp1, Temp2);//头文件里有定义，Temp1、Temp2为临时变量，在宏调用前定义
//
        for (i = 0; i < N; i++)
                printf("%06lu%.20lf + %.20lfi\n", i, BufferR, BufferV);
        return 0;
}

wzavr 发表于 2013-8-13 03:37:51

fft c记号一下电脑好好看看

fiddly 发表于 2013-8-13 06:57:38

fft，马克下，好东东！

zouzhichao 发表于 2013-8-13 06:59:18

fiddly 发表于 2013-8-13 06:57 static/image/common/back.gif
fft，马克下，好东东！

早上好{:lol:}

zouzhichao 发表于 2013-8-13 07:01:30

上传代码生成器的源代码：#include <stdio.h>
#include <math.h>
typedef signed char int8_;
typedef unsigned char uint8_;
typedef signed short int16_;
typedef unsigned short uint16_;
typedef signed long int32_;
typedef unsigned long uint32_;
typedef float flt_;
typedef double lflt_;

static int8_ CreateW(FILE* File, uint32_ N)
{
        uint32_ i, j;
        lflt_ Pi = 3.1415926535 * 2.0 / N, Temp = 0.0;
        for (i = N >> 1, j = 0; i; i--, j++, Temp += Pi)
        {
                if (!fprintf(File, (const char*)"#define _WR_%lu_%lu_ (%.20lf)\n", N, j, cos(Temp)))
                        break;
                if (!fprintf(File, (const char*)"#define _WV_%lu_%lu_ (%.20lf)\n", N, j, sin(Temp)))
                        break;
        }
        return i ? -1 : 0;
}
static uint8_ N2K(uint32_ N)
{
        uint8_ k;
        for (k = 0, N--; N ; N >>= 1, k++)
                ;
        return k;
}
static int8_ CreateReOrder(FILE* File, uint32_ N)
{
        uint32_ u, v, i, j, k = N >> 1;
        if (!fprintf(File, (const char*)"#define _Exchange(Buffer,A,B,Temp)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTemp = Buffer;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBuffer = Buffer;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBuffer = Temp;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _ReOrder(Buffer,Temp)\\\n{\\\n"))
                return -1;
        for (u = 0, i = 0, j = 1; j < N; j++)
        {
                for (v = k; i >= v; i-= v, v >>= 1)
                        ;
                i += v;
                if (j < i)
                {
                        if (!fprintf(File, (const char*)"\t_Exchange(Buffer, %lu, %lu, Temp);\\\n", j, i))
                                break;
                        u++;
                }
        }
        if (!fprintf(File, (const char*)"}\n"))
                return -1;
        return (N - (1 << ((N2K(N) + 1) >> 1))) >> 1 == u ? 0 : -1;
}
static int8_ CreateButter(FILE* File, uint32_ N, uint8_ K)
{
        uint32_ v, i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
        for (v = 0, i = N >> (K + 1); i; i--)
        {
                for (j = J; j > 0; j--, v++)
                {
                        if (!fprintf(File, (const char*)"\t_Butter(BufferR, BufferV, %lu, %lu, _WR_%lu_%lu_, _WV_%lu_%lu_, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1, N, (j - 1) * Y, N, (j - 1) * Y))
                                return -1;
                }
        }
        return (N >> 1) == v ? 0 : -1;
}
static int8_ CreateButters(FILE* File, uint32_ N)
{
        uint8_ k, K = N2K(N);
        if (!fprintf(File, (const char*)"#define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempR = (MR) * BufferR - (MV) * BufferV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempV = (MR) * BufferV + (MV) * BufferR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR = BufferR - TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV = BufferV - TempV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR += TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV += TempV;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _Butters(BufferR,BufferV,TempR,TempV)\\\n{\\\n"))
                return -1;
        for (k = 0; k < K; k++)
                if (CreateButter(File, N, k))
                        return -1;
        if (!fprintf(File, (const char*)"}\n"))
                return -1;
        return 0 == K - k ? 0: -1;
}

int8_ CreateDitFftFile(const char* FileName, uint32_ N)
{
        FILE* File = fopen(FileName, "w");
        if (!File)
                return -1;
        if (CreateW(File, N))
        {
                fclose(File);
                return -1;
        }
        if (CreateReOrder(File, N))
        {
                fclose(File);
                return -1;
        }
        if (CreateButters(File, N))
        {
                fclose(File);
                return -1;
        }
        fclose(File);
        return 0;
}

zouzhichao 发表于 2013-8-13 07:10:59

本帖最后由 zouzhichao 于 2013-8-13 07:13 编辑

要生成128点DIT-FFT的代码直接调用CreateDitFftFile((const char*)"dit-fft-128.h", 128);就行了，会生成一个“dit-fft-128.h”的文件，然后新建一个工程就可以做128点DIT-FFT运算了
工程模板：
//引用头文件
#include "dit-fft-128.h"

#define N 128

//假设数据已准备好
externlflt_ Temp1;//临时变量
externlflt_ Temp2;//临时变量
externlflt_ BufferR;//实部
externlflt_ BufferV;//虚部

//倒序操作，在"dit-fft-cal-1024.h"已经定义
_ReOrder(BufferR, Temp1);
//蝶形运算操作，在"dit-fft-cal-1024.h"已经定义
_Butters(BufferR, BufferV, Temp1, Temp2);

至此，DIT-FFT运算已经结束

这个代码仍然不够完善不够快，还有一些改进空间，算法上也没有加窗这些，自行摸索

zouzhichao 发表于 2013-8-13 07:18:23

此外，还有一点，不知道大家有什么好办法没？
1024点的模板很大，有近500K，7000行，用VC6.0在Debug模式很容易变异出来，在Release模式一编译就死机，GCC测试也没问题
根据VC6.0Debug和GCC测试结果，运算没问题，只是代码非常大，GCC下原来代码编译后只有7K，用模板编译的有800K左右

zouzhichao 发表于 2013-8-13 12:47:57

改进后的代码生成工具：#include <stdio.h>
#include <math.h>
typedef signed char int8_;
typedef unsigned char uint8_;
typedef signed short int16_;
typedef unsigned short uint16_;
typedef signed long int32_;
typedef unsigned long uint32_;
typedef float flt_;
typedef double lflt_;
/**********************************************************************************************
#define _W_N_n_ Wn
#define _Exchange(Buffer,A,B,Temp)\
{\
        Temp = Buffer;\
        Buffer = Buffer;\
        Buffer = Temp;\
}
#define _ReOrder(Buffer,Temp)\
{\
        _Exchange(Buffer, %lu, %lu, Temp);\
        _Exchange(Buffer, %lu, %lu, Temp);\
        ......
        ......
        _Exchange(Buffer, %lu, %lu, Temp);\
        _Exchange(Buffer, %lu, %lu, Temp);\
}
#define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)
{\
        TempR = (MR) * BufferR - (MV) * BufferV;\
        TempV = (MR) * BufferV + (MV) * BufferR;\
        BufferR = BufferR - TempR;\
        BufferV = BufferV - TempV;\
        BufferR += TempR;\
        BufferV += TempV;\
}
#define _Butter0(BufferR,BufferV,A,B,TempR,TempV)
{\
        TempR = BufferR;\
        TempV = BufferV;\
        BufferR = BufferR - TempR;\
        BufferV = BufferV - TempV;\
        BufferR += TempR;\
        BufferV += TempV;\
}
#define _Butter1(BufferR,BufferV,A,B,TempR,TempV)
{\
        TempR = 0 - BufferV;\
        TempV = BufferR;\
        BufferR = BufferR - TempR;\
        BufferV = BufferV - TempV;\
        BufferR += TempR;\
        BufferV += TempV;\
}
#define _Butters(BufferR,BufferV,TempR,TempV)\
{\
        _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
        _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
        ......
        ......
        _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
        _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
}
**********************************************************************************************/
static int8_ CreateW(FILE* File, uint32_ N)
{
        uint32_ i, j;
        lflt_ Pi = 3.1415926535 * 2.0 / N, Temp = 0.0;
        for (i = N >> 1, j = 0; i; i--, j++, Temp += Pi)
        {
                if (!fprintf(File, (const char*)"#define _WR_%lu_%lu_ (%.20lf)\n", N, j, cos(Temp)))
                        break;
                if (!fprintf(File, (const char*)"#define _WV_%lu_%lu_ (%.20lf)\n", N, j, sin(Temp)))
                        break;
        }
        return i ? -1 : 0;
}
static uint8_ N2K(uint32_ N)
{
        uint8_ k;
        for (k = 0, N--; N ; N >>= 1, k++)
                ;
        return k;
}
static int8_ CreateReOrder(FILE* File, uint32_ N)
{
        uint32_ u, v, i, j, k = N >> 1;
        if (!fprintf(File, (const char*)"#define _Exchange(Buffer,A,B,Temp)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTemp = Buffer;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBuffer = Buffer;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBuffer = Temp;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _ReOrder(Buffer,Temp)\\\n{\\\n"))
                return -1;
        for (u = 0, i = 0, j = 1; j < N; j++)
        {
                for (v = k; i >= v; i-= v, v >>= 1)
                        ;
                i += v;
                if (j < i)
                {
                        if (!fprintf(File, (const char*)"\t_Exchange(Buffer, %lu, %lu, Temp);\\\n", j, i))
                                break;
                        u++;
                }
        }
        if (!fprintf(File, (const char*)"}\n"))
                return -1;
        return (N - (1 << ((N2K(N) + 1) >> 1))) >> 1 == u ? 0 : -1;
}
static int8_ CreateButter(FILE* File, uint32_ N, uint8_ K)
{
        uint32_ v, u, i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
        for (v = 0, i = N >> (K + 1); i; i--)
        {
                for (j = J; j > 0; j--, v++)
                {
                        u = (j - 1) * Y;
                        if (0 == u)
                        {
                                if (!fprintf(File, (const char*)"\t_Butter0(BufferR, BufferV, %lu, %lu, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1))
                                        return -1;
                        }
                        else if (N >> 2 == u)
                        {
                                if (!fprintf(File, (const char*)"\t_Butter1(BufferR, BufferV, %lu, %lu, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1))
                                        return -1;
                        }
                        else
                        {
                                if (!fprintf(File, (const char*)"\t_Butter(BufferR, BufferV, %lu, %lu, _WR_%lu_%lu_, _WV_%lu_%lu_, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1, N, u, N, u))
                                        return -1;
                        }       
                }
        }
        return (N >> 1) == v ? 0 : -1;
}
static int8_ CreateButters(FILE* File, uint32_ N)
{
        uint8_ k, K = N2K(N);
        if (!fprintf(File, (const char*)"#define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempR = (MR) * BufferR - (MV) * BufferV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempV = (MR) * BufferV + (MV) * BufferR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR = BufferR - TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV = BufferV - TempV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR += TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV += TempV;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _Butter0(BufferR,BufferV,A,B,TempR,TempV)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempR = BufferR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempV = BufferV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR = BufferR - TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV = BufferV - TempV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR += TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV += TempV;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _Butter1(BufferR,BufferV,A,B,TempR,TempV)\\\n{\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempR = 0 - BufferV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tTempV = BufferR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR = BufferR - TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV = BufferV - TempV;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferR += TempR;\\\n"))
                return -1;
        if (!fprintf(File, (const char*)"\tBufferV += TempV;\\\n}\n"))
                return -1;
        if (!fprintf(File, (const char*)"#define _Butters(BufferR,BufferV,TempR,TempV)\\\n{\\\n"))
                return -1;
        for (k = 0; k < K; k++)
                if (CreateButter(File, N, k))
                        return -1;
        if (!fprintf(File, (const char*)"}\n"))
                return -1;
        return 0 == K - k ? 0: -1;
}

int8_ CreateDitFftFile(const char* FileName, uint32_ N)
{
        FILE* File = fopen(FileName, "w");
        if (!File)
                return -1;
        if (CreateW(File, N))
        {
                fclose(File);
                return -1;
        }
        if (CreateReOrder(File, N))
        {
                fclose(File);
                return -1;
        }
        if (CreateButters(File, N))
        {
                fclose(File);
                return -1;
        }
        fclose(File);
        return 0;
}

277955973 发表于 2013-8-13 12:51:07

C fft PC MARK

zouzhichao 发表于 2013-8-13 12:52:33

277955973 发表于 2013-8-13 12:51 static/image/common/back.gif
C fft PC MARK

{:smile:} {:smile:}

xjl00 发表于 2013-8-14 14:22:18

赞一个！！

daidaimumu 发表于 2013-8-14 17:30:42

{:smile:}不错，研究研究

zouzhichao 发表于 2013-8-14 18:05:50

xjl00 发表于 2013-8-14 14:22 static/image/common/back.gif
赞一个！！

承蒙谬赞

jxcylxh 发表于 2013-8-14 19:29:21

非常好，顶一下。。

zouzhichao 发表于 2013-8-14 19:40:50

jxcylxh 发表于 2013-8-14 19:29 static/image/common/back.gif
非常好，顶一下。。

多谢帮顶，花了不少功夫，希望对大家有帮助

Jimmyxu 发表于 2013-8-14 20:17:52

不错~~~~~~~~~

ddds 发表于 2013-8-14 20:30:20

这个很不错的资料啊，最近也在研究用单片机实现FFT啊，不知道楼主的代码单片机跑的过来吗

zouzhichao 发表于 2013-8-14 20:40:25

ddds 发表于 2013-8-14 20:30 static/image/common/back.gif
这个很不错的资料啊，最近也在研究用单片机实现FFT啊，不知道楼主的代码单片机跑的过来吗 ...

我是照着《数字信号处理》这本书一步步写的代码，用的浮点型的，估计单片机跑要慢的很，
代码经过GCC，TCC，VC都测试过的，不挑平台，移植工作几乎没有
要在单片机里运行，如果对速度有要求，估计要考虑改成定点运算
PS：我手上有AVR的DIT-FFT的汇编源码（网上down的），不过我没有弄明白，不保证可用

ddds 发表于 2013-8-14 20:44:54

zouzhichao 发表于 2013-8-14 20:40 static/image/common/back.gif
我是照着《数字信号处理》这本书一步步写的代码，用的浮点型的，估计单片机跑要慢的很，
代码经过GCC，TC ...

呵呵，汇编的从学到现在都快忘记的差不多了，不过慢慢研究楼主的代码，这个可能要DSP才能跑的过来

zouzhichao 发表于 2013-8-14 20:45:09

ddds 发表于 2013-8-14 20:30 static/image/common/back.gif
这个很不错的资料啊，最近也在研究用单片机实现FFT啊，不知道楼主的代码单片机跑的过来吗 ...

帖子里提出了两种DIT-FFT的风格，一种是空间换速度，一种是速度换空间，单片机如果要跑点数多的话，前一种会消耗不少的Flash空间，但是速度会快一些，RAM消耗少一些

zouzhichao 发表于 2013-8-14 20:45:44

ddds 发表于 2013-8-14 20:30 static/image/common/back.gif
这个很不错的资料啊，最近也在研究用单片机实现FFT啊，不知道楼主的代码单片机跑的过来吗 ...

帖子里提出了两种DIT-FFT的风格，一种是空间换速度，一种是速度换空间，单片机如果要跑点数多的话，前一种会消耗不少的Flash空间，但是速度会快一些，RAM消耗少一些

zouzhichao 发表于 2013-8-14 20:46:17

ddds 发表于 2013-8-14 20:44 static/image/common/back.gif
呵呵，汇编的从学到现在都快忘记的差不多了，不过慢慢研究楼主的代码，这个可能要DSP才能跑的过来 ...

你要运算多少点的FFT？一秒需要运算多少次？

zouzhichao 发表于 2013-8-14 20:51:59

ddds 发表于 2013-8-14 20:44 static/image/common/back.gif
呵呵，汇编的从学到现在都快忘记的差不多了，不过慢慢研究楼主的代码，这个可能要DSP才能跑的过来 ...

我的那点汇编现在也不灵光了，主要是用汇编写对平台限制太死，又没有特别的动力，总感觉不值当，现在主要还是喜欢用C实现一些算法（几乎100%可移植），目前搞定的有曲线拟合，DIT-FFT，下一步准备攻占DIF-FFT

zouzhichao 发表于 2013-8-14 21:05:12

用WinAVR GCC + AVR Studio跑了一下仿真，256点的DIT-FFT用了1.5S，4M晶振

zouzhichao 发表于 2013-8-14 21:07:21

32点的DIT-FFT用了0.14S，4M晶振

zouzhichao 发表于 2013-8-14 21:19:17

用代码生成向导生成的32点的DIT-FFT代码比较快，用了0.07S，4M晶振，是前一个速度的两倍

ddds 发表于 2013-8-14 21:20:20

zouzhichao 发表于 2013-8-14 21:07 static/image/common/back.gif
32点的DIT-FFT用了0.14S，4M晶振

这个算很不错啦，呵呵，改天有空移植跑一下

zouzhichao 发表于 2013-8-14 21:21:52

不过后者编译后有21kb，前者编译后5kb多

zouzhichao 发表于 2013-8-14 21:29:10

用代码生成向导生成的128点的DIT-FFT代码比较快，用了0.4S，4M晶振

guxingganyue 发表于 2013-8-14 21:38:35

写的不错，顶一个，一直没弄懂那个蝴蝶图是啥意思

zouzhichao 发表于 2013-8-14 21:45:29

guxingganyue 发表于 2013-8-14 21:38 static/image/common/back.gif
写的不错，顶一个，一直没弄懂那个蝴蝶图是啥意思

原理我也一知半解，就是照着书上的流程一步步转换成C语言代码实现

zouzhichao 发表于 2013-8-14 21:48:23

guxingganyue 发表于 2013-8-14 21:38 static/image/common/back.gif
写的不错，顶一个，一直没弄懂那个蝴蝶图是啥意思

哈哈，又遇见你了

guxingganyue 发表于 2013-8-15 08:50:51

zouzhichao 发表于 2013-8-14 21:48 static/image/common/back.gif
哈哈，又遇见你了

呵呵，你回复我两次

chinaboy25 发表于 2013-8-15 08:58:58

make,make,,

cqfeiyu 发表于 2013-8-15 09:05:39

这个不错,老师都这么教的话,中国赶超日本还是有望的.

xjmlfm1 发表于 2013-8-15 09:25:57

楼主说的0.07S是在哪测的？单片机里面吗？
我们实际应用时，32点的FFT，在22M晶振的51单片机里面，跑一次用时6ms左右。
在STM32里的72M晶振时，64点FFT耗时300us。

Pjm2008 发表于 2013-8-15 09:26:07

不错。有时间再看看

robyn 发表于 2013-8-15 09:47:10

顶楼主，不过楼主能告诉我一下，是看的哪本书吗？我也想好好学习一下这个

zouzhichao 发表于 2013-8-15 12:39:30

xjmlfm1 发表于 2013-8-15 09:25 static/image/common/back.gif
楼主说的0.07S是在哪测的？单片机里面吗？
我们实际应用时，32点的FFT，在22M晶振的51单片机里面，跑一次用 ...

AVR Studio里跑的仿真

zouzhichao 发表于 2013-8-15 12:42:03

xjmlfm1 发表于 2013-8-15 09:25 static/image/common/back.gif
楼主说的0.07S是在哪测的？单片机里面吗？
我们实际应用时，32点的FFT，在22M晶振的51单片机里面，跑一次用 ...

你手上的的代码应该比我写的质量要好很多

zouzhichao 发表于 2013-8-15 12:44:39

robyn 发表于 2013-8-15 09:47 static/image/common/back.gif
顶楼主，不过楼主能告诉我一下，是看的哪本书吗？我也想好好学习一下这个 ...

《数字信号处理》第三版 高西全 丁玉美 编著 西安电子科技大学出版社
这个是们学校集成电路专业用的教科书，不过我是自动化专业的，没有修这个课，书是从我室友那拿的

zouzhichao 发表于 2013-8-15 12:45:35

cqfeiyu 发表于 2013-8-15 09:05 static/image/common/back.gif
这个不错,老师都这么教的话,中国赶超日本还是有望的.

多谢夸奖，不过实在惭愧

zouzhichao 发表于 2013-8-15 12:47:48

xjmlfm1 发表于 2013-8-15 09:25 static/image/common/back.gif
楼主说的0.07S是在哪测的？单片机里面吗？
我们实际应用时，32点的FFT，在22M晶振的51单片机里面，跑一次用 ...

不知道你们用的是整型运算还是浮点运算？

xjmlfm1 发表于 2013-8-15 15:43:37

zouzhichao 发表于 2013-8-15 12:47 static/image/common/back.gif
不知道你们用的是整型运算还是浮点运算？

我们用的是浮点的。
我一直想把它转成定点的，但是51里面没法定义long long型变量，转不过来。

xjmlfm1 发表于 2013-8-15 15:44:22

zouzhichao 发表于 2013-8-15 12:39 static/image/common/back.gif
AVR Studio里跑的仿真

这种虚拟的仿真可能会有一定的误差，也不一定是代码的事。

wmm20031015 发表于 2013-8-15 17:07:47

好东东，顶起来

wcm_e 发表于 2013-8-15 17:12:40

mark 学习

zouzhichao 发表于 2013-8-15 18:48:09

wmm20031015 发表于 2013-8-15 17:07 static/image/common/back.gif
好东东，顶起来

谢谢帮顶

zouzhichao 发表于 2013-8-15 18:55:48

xjmlfm1 发表于 2013-8-15 15:44 static/image/common/back.gif
这种虚拟的仿真可能会有一定的误差，也不一定是代码的事。

刚工作，开发板还在包裹里没打开，等心情好了再拿出来跑跑试试

rifjft 发表于 2013-8-15 19:48:10

弄完这个，头发掉了几根？

zouzhichao 发表于 2013-8-15 22:19:51

rifjft 发表于 2013-8-15 19:48 static/image/common/back.gif
弄完这个，头发掉了几根？

还好啦，照着书的流程写的，没啥动脑子的地方

zouzhichao 发表于 2013-8-15 22:22:34

rifjft 发表于 2013-8-15 19:48 static/image/common/back.gif
弄完这个，头发掉了几根？

这个里面的最小二乘法拟合指数曲线才伤了很久的脑筋
http://www.amobbs.com/thread-5535356-1-1.html

rifjft 发表于 2013-8-16 00:49:02

嗯，看得出都是花了不少心思折腾。应届生中能下这份功夫的少见，难得的是还喜欢逛论坛{:lol:}

如若楼主之流能在应届生中都占半数，未来可谓人才济济

已荐精

zouzhichao 发表于 2013-8-16 07:14:35

rifjft 发表于 2013-8-16 00:49 static/image/common/back.gif
嗯，看得出都是花了不少心思折腾。应届生中能下这份功夫的少见，难得的是还喜欢逛论坛

如若楼主之 ...

逛论坛有四年了，不过以前是潜水党，在论坛学到很多东西，也希望自己能带给别人收获吧，以后有什么好的东西也会分享出来的

zouzhichao 发表于 2013-8-16 07:15:12

rifjft 发表于 2013-8-16 00:49 static/image/common/back.gif
嗯，看得出都是花了不少心思折腾。应届生中能下这份功夫的少见，难得的是还喜欢逛论坛

如若楼主之 ...

还有多谢推荐

pipi516 发表于 2013-8-16 08:40:16

FFT,先MARK一下!

zhanyanqiang 发表于 2013-8-18 16:43:56

好东西，，，，

zouzhichao 发表于 2013-8-18 16:57:11

zhanyanqiang 发表于 2013-8-18 16:43 static/image/common/back.gif
好东西，，，，

好东西自然要分享

lixianghappy 发表于 2013-8-18 20:01:15

mark                                    

acai039033 发表于 2013-8-19 14:53:44

这么好的东西 不顶说不过去了

siyeb 发表于 2013-8-19 15:17:53

lz回头再把校正加上就更好了。

zouzhichao 发表于 2013-8-19 22:37:40

siyeb 发表于 2013-8-19 15:17 static/image/common/back.gif
lz回头再把校正加上就更好了。

工作了时间上不如大学时随意了，想抽点时间高点东西都得从牙缝里挤，忙完这两周的任务再考虑弄这个吧

zouzhichao 发表于 2013-8-19 22:39:31

siyeb 发表于 2013-8-19 15:17 static/image/common/back.gif
lz回头再把校正加上就更好了。

BOSS派任务下来了，得加紧了

siyeb 发表于 2013-8-20 10:55:08

zouzhichao 发表于 2013-8-19 22:37 static/image/common/back.gif
工作了时间上不如大学时随意了，想抽点时间高点东西都得从牙缝里挤，忙完这两周的任务再考虑弄这个吧 ...

所以楼主加油啊~

wangfei1956 发表于 2013-8-27 22:24:23

好东西，收藏了！

zouzhichao 发表于 2013-8-27 22:26:56

wangfei1956 发表于 2013-8-27 22:24 static/image/common/back.gif
好东西，收藏了！

欢迎收藏

晓毕8 发表于 2013-8-28 08:46:07

赞一个，最近也想搞单片机的FFT，试试下看看能不能成功{:lol:}

binghe167 发表于 2013-8-28 08:58:46

不错，感谢楼主分享

404710520 发表于 2013-8-28 09:11:31

谢谢楼主共享

rainy_T 发表于 2013-8-29 01:04:19

这个果断得顶起啊

loves6036 发表于 2013-8-29 12:29:00

谢谢楼主分享，真是宝贝啊，收藏了

zouzhichao 发表于 2013-8-29 12:34:07

loves6036 发表于 2013-8-29 12:29 static/image/common/back.gif
谢谢楼主分享，真是宝贝啊，收藏了

承蒙夸奖

fangchangqing 发表于 2013-8-30 12:09:50

这个好。正需要

zouzhichao 发表于 2013-8-30 12:50:18

fangchangqing 发表于 2013-8-30 12:09 static/image/common/back.gif
这个好。正需要

尽管拿去

fangchangqing 发表于 2013-8-30 13:39:53

嘿嘿{:smile:}{:smile:}{:smile:}{:smile:}

zyw19987 发表于 2013-8-30 13:51:55

zouzhichao 发表于 2013-8-30 12:50 static/image/common/back.gif
尽管拿去

牛！
录个视频讲解吧{:lol:}

tscyds 发表于 2013-8-30 22:06:32

很不错啊，学习学习

竹风xu 发表于 2013-8-31 19:41:44

灰常好！！！感谢

407043680 发表于 2013-9-8 18:53:15

蝶形算法厉害

kevin_ldw 发表于 2013-9-8 22:38:30

非常不错，学习了

touch_mcu 发表于 2013-9-9 01:11:32

马克一下。

郭本宏 发表于 2013-9-9 07:58:59

搞得不错，好好

qinkaiabc 发表于 2013-9-10 16:38:26

收藏备用！

landmuto 发表于 2013-9-14 19:43:22

fft，mark！

zouzhichao 发表于 2013-9-15 11:09:39

landmuto 发表于 2013-9-14 19:43 static/image/common/back.gif
fft，mark！

欢迎MARK

armok 发表于 2013-9-19 12:05:44

COOL !

kmani 发表于 2013-9-19 12:39:19

FFT 好东西！

zouzhichao 发表于 2013-9-23 12:40:43

{:biggrin:}好意外！！！居然置酷了，Mygod

szmini2006 发表于 2013-9-23 13:58:27

学习了。。。

moen 发表于 2013-9-23 14:07:26

好帖留名

TongIC 发表于 2013-11-22 22:14:08

楼主强大，非常好的代码分享！

xiaowenshao 发表于 2013-11-23 08:58:06

niubility，大家的fft都应用在哪里了呢

查看完整版本: FFT算法C代码