FFT算法C代码

zouzhichao

基本蝶形运算单元：


注：上图中的运算均为复数运算。
复数表示形式：
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[0].Real = 1;
        Modulus[0].Virtual = 0;
        Modulus[1].Real = cos(Pi);
        Modulus[1].Virtual = sin(Pi);
        for (i = 1, j = 2, N >>= 2; i < N; i++, j += 2)
        {
                Modulus[j].Real = 2.0 * Modulus.Real * Modulus.Real - 1;
                Modulus[j].Virtual = 2.0 * Modulus.Real * Modulus.Virtual;
                Modulus[j + 1].Real = Modulus[1].Real * Modulus[j].Real - Modulus[1].Virtual * Modulus[j].Virtual;
                Modulus[j + 1].Virtual = Modulus[1].Real * Modulus[j].Virtual + Modulus[1].Virtual * Modulus[j].Real;
        }
}

以空间换时间
将内存地址A的1024个数复制到内存地址B处：
（1）
for (i = 0; i < 1024; i++)
        B = A;
（2）
B[0] = A[0];
B[1] = A[1];
B[2] = A[2];
B[3] = A[3];
……
……
B[1022] = A[1022];
B[1023] = A[1023];
显然，（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[1], Data[32]);
        _Change(Data[2], Data[16]);
        _Change(Data[3], Data[48]);
        _Change(Data[4], Data[8]);
        _Change(Data[5], Data[40]);
        _Change(Data[6], Data[24]);
        _Change(Data[7], Data[56]);
        _Change(Data[9], Data[36]);
        _Change(Data[10], Data[20]);
        _Change(Data[11], Data[52]);
        _Change(Data[13], Data[44]);
        _Change(Data[14], Data[28]);
        _Change(Data[15], Data[60]);
        _Change(Data[17], Data[34]);
        _Change(Data[19], Data[50]);
        _Change(Data[21], Data[42]);
        _Change(Data[22], Data[26]);
        _Change(Data[23], Data[58]);
        _Change(Data[25], Data[38]);
        _Change(Data[27], Data[54]);
        _Change(Data[29], Data[46]);
        _Change(Data[31], Data[62]);
        _Change(Data[35], Data[49]);
        _Change(Data[37], Data[41]);
        _Change(Data[39], Data[57]);
        _Change(Data[43], Data[53]);
        _Change(Data[47], Data[61]);
        _Change(Data[55], Data[59]);
}
（2）对于FftK，第0级运算有特殊性，可以特殊处理。
（3）对于整个FftK，FFT，运算顺序，运算次数都是固定的，因此可以取消循环，做成顺序执行的程序，如Change函数一样。对于这点会增加很多运算速度，同事也占用很长的代码，同时需要程序员做写更多的工作，这部分可以参考单片机小精灵自己做一个代码生成器。

zouzhichao

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

cos

好东东,期待中.

zouzhichao

cos 发表于 2013-8-12 23:20
好东东,期待中.

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

leifeng

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

zouzhichao

leifeng 发表于 2013-8-12 23:47
这个真的很历害!楼主NB!

承蒙谬赞，愧不敢当

zouzhichao

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

VC6.0工程：

zouzhichao

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

zouzhichao

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

测试模板：
#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[N], BufferV[N] = {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

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

fiddly

fft，马克下，好东东！

zouzhichao

fiddly 发表于 2013-8-13 06:57
fft，马克下，好东东！

早上好

zouzhichao

上传代码生成器的源代码：

1. #include <stdio.h>
   
2. #include <math.h>
   
3. typedef signed char int8_;
   
4. typedef unsigned char uint8_;
   
5. typedef signed short int16_;
   
6. typedef unsigned short uint16_;
   
7. typedef signed long int32_;
   
8. typedef unsigned long uint32_;
   
9. typedef float flt_;
   
10. typedef double lflt_;
    
11. 
    
12. static int8_ CreateW(FILE* File, uint32_ N)
    
13. {
    
14.         uint32_ i, j;
    
15.         lflt_ Pi = 3.1415926535 * 2.0 / N, Temp = 0.0;
    
16.         for (i = N >> 1, j = 0; i; i--, j++, Temp += Pi)
    
17.         {
    
18.                 if (!fprintf(File, (const char*)"#define _WR_%lu_%lu_ (%.20lf)\n", N, j, cos(Temp)))
    
19.                         break;
    
20.                 if (!fprintf(File, (const char*)"#define _WV_%lu_%lu_ (%.20lf)\n", N, j, sin(Temp)))
    
21.                         break;
    
22.         }
    
23.         return i ? -1 : 0;
    
24. }
    
25. static uint8_ N2K(uint32_ N)
    
26. {
    
27.         uint8_ k;
    
28.         for (k = 0, N--; N ; N >>= 1, k++)
    
29.                 ;
    
30.         return k;
    
31. }
    
32. static int8_ CreateReOrder(FILE* File, uint32_ N)
    
33. {
    
34.         uint32_ u, v, i, j, k = N >> 1;
    
35.         if (!fprintf(File, (const char*)"#define _Exchange(Buffer,A,B,Temp)\\\n{\\\n"))
    
36.                 return -1;
    
37.         if (!fprintf(File, (const char*)"\tTemp = Buffer[A];\\\n"))
    
38.                 return -1;
    
39.         if (!fprintf(File, (const char*)"\tBuffer[A] = Buffer[B];\\\n"))
    
40.                 return -1;
    
41.         if (!fprintf(File, (const char*)"\tBuffer[B] = Temp;\\\n}\n"))
    
42.                 return -1;
    
43.         if (!fprintf(File, (const char*)"#define _ReOrder(Buffer,Temp)\\\n{\\\n"))
    
44.                 return -1;
    
45.         for (u = 0, i = 0, j = 1; j < N; j++)
    
46.         {
    
47.                 for (v = k; i >= v; i-= v, v >>= 1)
    
48.                         ;
    
49.                 i += v;
    
50.                 if (j < i)
    
51.                 {
    
52.                         if (!fprintf(File, (const char*)"\t_Exchange(Buffer, %lu, %lu, Temp);\\\n", j, i))
    
53.                                 break;
    
54.                         u++;
    
55.                 }
    
56.         }
    
57.         if (!fprintf(File, (const char*)"}\n"))
    
58.                 return -1;
    
59.         return (N - (1 << ((N2K(N) + 1) >> 1))) >> 1 == u ? 0 : -1;
    
60. }
    
61. static int8_ CreateButter(FILE* File, uint32_ N, uint8_ K)
    
62. {
    
63.         uint32_ v, i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
    
64.         for (v = 0, i = N >> (K + 1); i; i--)
    
65.         {
    
66.                 for (j = J; j > 0; j--, v++)
    
67.                 {
    
68.                         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))
    
69.                                 return -1;
    
70.                 }
    
71.         }
    
72.         return (N >> 1) == v ? 0 : -1;
    
73. }
    
74. static int8_ CreateButters(FILE* File, uint32_ N)
    
75. {
    
76.         uint8_ k, K = N2K(N);
    
77.         if (!fprintf(File, (const char*)"#define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)\\\n{\\\n"))
    
78.                 return -1;
    
79.         if (!fprintf(File, (const char*)"\tTempR = (MR) * BufferR[B] - (MV) * BufferV[B];\\\n"))
    
80.                 return -1;
    
81.         if (!fprintf(File, (const char*)"\tTempV = (MR) * BufferV[B] + (MV) * BufferR[B];\\\n"))
    
82.                 return -1;
    
83.         if (!fprintf(File, (const char*)"\tBufferR[B] = BufferR[A] - TempR;\\\n"))
    
84.                 return -1;
    
85.         if (!fprintf(File, (const char*)"\tBufferV[B] = BufferV[A] - TempV;\\\n"))
    
86.                 return -1;
    
87.         if (!fprintf(File, (const char*)"\tBufferR[A] += TempR;\\\n"))
    
88.                 return -1;
    
89.         if (!fprintf(File, (const char*)"\tBufferV[A] += TempV;\\\n}\n"))
    
90.                 return -1;
    
91.         if (!fprintf(File, (const char*)"#define _Butters(BufferR,BufferV,TempR,TempV)\\\n{\\\n"))
    
92.                 return -1;
    
93.         for (k = 0; k < K; k++)
    
94.                 if (CreateButter(File, N, k))
    
95.                         return -1;
    
96.         if (!fprintf(File, (const char*)"}\n"))
    
97.                 return -1;
    
98.         return 0 == K - k ? 0: -1;
    
99. }
    
100. 
     
101. int8_ CreateDitFftFile(const char* FileName, uint32_ N)
     
102. {
     
103.         FILE* File = fopen(FileName, "w");
     
104.         if (!File)
     
105.                 return -1;
     
106.         if (CreateW(File, N))
     
107.         {
     
108.                 fclose(File);
     
109.                 return -1;
     
110.         }
     
111.         if (CreateReOrder(File, N))
     
112.         {
     
113.                 fclose(File);
     
114.                 return -1;
     
115.         }
     
116.         if (CreateButters(File, N))
     
117.         {
     
118.                 fclose(File);
     
119.                 return -1;
     
120.         }
     
121.         fclose(File);
     
122.         return 0;
     
123. }

复制代码

zouzhichao

要生成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

//假设数据已准备好
extern  lflt_ Temp1;//临时变量
extern  lflt_ Temp2;//临时变量
extern  lflt_ BufferR[N];//实部
extern  lflt_ BufferV[N];//虚部

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

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

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

zouzhichao

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

zouzhichao

改进后的代码生成工具：

1. #include <stdio.h>
   
2. #include <math.h>
   
3. typedef signed char int8_;
   
4. typedef unsigned char uint8_;
   
5. typedef signed short int16_;
   
6. typedef unsigned short uint16_;
   
7. typedef signed long int32_;
   
8. typedef unsigned long uint32_;
   
9. typedef float flt_;
   
10. typedef double lflt_;
    
11. /**********************************************************************************************
    
12. #define _W_N_n_ Wn
    
13. #define _Exchange(Buffer,A,B,Temp)\
    
14. {\
    
15.         Temp = Buffer[A];\
    
16.         Buffer[A] = Buffer[B];\
    
17.         Buffer[B] = Temp;\
    
18. }
    
19. #define _ReOrder(Buffer,Temp)\
    
20. {\
    
21.         _Exchange(Buffer, %lu, %lu, Temp);\
    
22.         _Exchange(Buffer, %lu, %lu, Temp);\
    
23.         ......
    
24.         ......
    
25.         _Exchange(Buffer, %lu, %lu, Temp);\
    
26.         _Exchange(Buffer, %lu, %lu, Temp);\
    
27. }
    
28. #define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)
    
29. {\
    
30.         TempR = (MR) * BufferR[B] - (MV) * BufferV[B];\
    
31.         TempV = (MR) * BufferV[B] + (MV) * BufferR[B];\
    
32.         BufferR[B] = BufferR[A] - TempR;\
    
33.         BufferV[B] = BufferV[A] - TempV;\
    
34.         BufferR[A] += TempR;\
    
35.         BufferV[A] += TempV;\
    
36. }
    
37. #define _Butter0(BufferR,BufferV,A,B,TempR,TempV)
    
38. {\
    
39.         TempR = BufferR[B];\
    
40.         TempV = BufferV[B];\
    
41.         BufferR[B] = BufferR[A] - TempR;\
    
42.         BufferV[B] = BufferV[A] - TempV;\
    
43.         BufferR[A] += TempR;\
    
44.         BufferV[A] += TempV;\
    
45. }
    
46. #define _Butter1(BufferR,BufferV,A,B,TempR,TempV)
    
47. {\
    
48.         TempR = 0 - BufferV[B];\
    
49.         TempV = BufferR[B];\
    
50.         BufferR[B] = BufferR[A] - TempR;\
    
51.         BufferV[B] = BufferV[A] - TempV;\
    
52.         BufferR[A] += TempR;\
    
53.         BufferV[A] += TempV;\
    
54. }
    
55. #define _Butters(BufferR,BufferV,TempR,TempV)\
    
56. {\
    
57.         _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
    
58.         _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
    
59.         ......
    
60.         ......
    
61.         _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
    
62.         _Butter(BufferR, BufferV, %lu, %lu, %lf, %lf, TempR, TempV);\
    
63. }
    
64. **********************************************************************************************/
    
65. static int8_ CreateW(FILE* File, uint32_ N)
    
66. {
    
67.         uint32_ i, j;
    
68.         lflt_ Pi = 3.1415926535 * 2.0 / N, Temp = 0.0;
    
69.         for (i = N >> 1, j = 0; i; i--, j++, Temp += Pi)
    
70.         {
    
71.                 if (!fprintf(File, (const char*)"#define _WR_%lu_%lu_ (%.20lf)\n", N, j, cos(Temp)))
    
72.                         break;
    
73.                 if (!fprintf(File, (const char*)"#define _WV_%lu_%lu_ (%.20lf)\n", N, j, sin(Temp)))
    
74.                         break;
    
75.         }
    
76.         return i ? -1 : 0;
    
77. }
    
78. static uint8_ N2K(uint32_ N)
    
79. {
    
80.         uint8_ k;
    
81.         for (k = 0, N--; N ; N >>= 1, k++)
    
82.                 ;
    
83.         return k;
    
84. }
    
85. static int8_ CreateReOrder(FILE* File, uint32_ N)
    
86. {
    
87.         uint32_ u, v, i, j, k = N >> 1;
    
88.         if (!fprintf(File, (const char*)"#define _Exchange(Buffer,A,B,Temp)\\\n{\\\n"))
    
89.                 return -1;
    
90.         if (!fprintf(File, (const char*)"\tTemp = Buffer[A];\\\n"))
    
91.                 return -1;
    
92.         if (!fprintf(File, (const char*)"\tBuffer[A] = Buffer[B];\\\n"))
    
93.                 return -1;
    
94.         if (!fprintf(File, (const char*)"\tBuffer[B] = Temp;\\\n}\n"))
    
95.                 return -1;
    
96.         if (!fprintf(File, (const char*)"#define _ReOrder(Buffer,Temp)\\\n{\\\n"))
    
97.                 return -1;
    
98.         for (u = 0, i = 0, j = 1; j < N; j++)
    
99.         {
    
100.                 for (v = k; i >= v; i-= v, v >>= 1)
     
101.                         ;
     
102.                 i += v;
     
103.                 if (j < i)
     
104.                 {
     
105.                         if (!fprintf(File, (const char*)"\t_Exchange(Buffer, %lu, %lu, Temp);\\\n", j, i))
     
106.                                 break;
     
107.                         u++;
     
108.                 }
     
109.         }
     
110.         if (!fprintf(File, (const char*)"}\n"))
     
111.                 return -1;
     
112.         return (N - (1 << ((N2K(N) + 1) >> 1))) >> 1 == u ? 0 : -1;
     
113. }
     
114. static int8_ CreateButter(FILE* File, uint32_ N, uint8_ K)
     
115. {
     
116.         uint32_ v, u, i, j, J = 1 << K, X = 1 << (K + 1), Y = 1 << (N2K(N) - K - 1);
     
117.         for (v = 0, i = N >> (K + 1); i; i--)
     
118.         {
     
119.                 for (j = J; j > 0; j--, v++)
     
120.                 {
     
121.                         u = (j - 1) * Y;
     
122.                         if (0 == u)
     
123.                         {
     
124.                                 if (!fprintf(File, (const char*)"\t_Butter0(BufferR, BufferV, %lu, %lu, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1))
     
125.                                         return -1;
     
126.                         }
     
127.                         else if (N >> 2 == u)
     
128.                         {
     
129.                                 if (!fprintf(File, (const char*)"\t_Butter1(BufferR, BufferV, %lu, %lu, TempR, TempV);\\\n", (i - 1) * X + j - 1, i * X + j - J - 1))
     
130.                                         return -1;
     
131.                         }
     
132.                         else
     
133.                         {
     
134.                                 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))
     
135.                                         return -1;
     
136.                         }       
     
137.                 }
     
138.         }
     
139.         return (N >> 1) == v ? 0 : -1;
     
140. }
     
141. static int8_ CreateButters(FILE* File, uint32_ N)
     
142. {
     
143.         uint8_ k, K = N2K(N);
     
144.         if (!fprintf(File, (const char*)"#define _Butter(BufferR,BufferV,A,B,MR,MV,TempR,TempV)\\\n{\\\n"))
     
145.                 return -1;
     
146.         if (!fprintf(File, (const char*)"\tTempR = (MR) * BufferR[B] - (MV) * BufferV[B];\\\n"))
     
147.                 return -1;
     
148.         if (!fprintf(File, (const char*)"\tTempV = (MR) * BufferV[B] + (MV) * BufferR[B];\\\n"))
     
149.                 return -1;
     
150.         if (!fprintf(File, (const char*)"\tBufferR[B] = BufferR[A] - TempR;\\\n"))
     
151.                 return -1;
     
152.         if (!fprintf(File, (const char*)"\tBufferV[B] = BufferV[A] - TempV;\\\n"))
     
153.                 return -1;
     
154.         if (!fprintf(File, (const char*)"\tBufferR[A] += TempR;\\\n"))
     
155.                 return -1;
     
156.         if (!fprintf(File, (const char*)"\tBufferV[A] += TempV;\\\n}\n"))
     
157.                 return -1;
     
158.         if (!fprintf(File, (const char*)"#define _Butter0(BufferR,BufferV,A,B,TempR,TempV)\\\n{\\\n"))
     
159.                 return -1;
     
160.         if (!fprintf(File, (const char*)"\tTempR = BufferR[B];\\\n"))
     
161.                 return -1;
     
162.         if (!fprintf(File, (const char*)"\tTempV = BufferV[B];\\\n"))
     
163.                 return -1;
     
164.         if (!fprintf(File, (const char*)"\tBufferR[B] = BufferR[A] - TempR;\\\n"))
     
165.                 return -1;
     
166.         if (!fprintf(File, (const char*)"\tBufferV[B] = BufferV[A] - TempV;\\\n"))
     
167.                 return -1;
     
168.         if (!fprintf(File, (const char*)"\tBufferR[A] += TempR;\\\n"))
     
169.                 return -1;
     
170.         if (!fprintf(File, (const char*)"\tBufferV[A] += TempV;\\\n}\n"))
     
171.                 return -1;
     
172.         if (!fprintf(File, (const char*)"#define _Butter1(BufferR,BufferV,A,B,TempR,TempV)\\\n{\\\n"))
     
173.                 return -1;
     
174.         if (!fprintf(File, (const char*)"\tTempR = 0 - BufferV[B];\\\n"))
     
175.                 return -1;
     
176.         if (!fprintf(File, (const char*)"\tTempV = BufferR[B];\\\n"))
     
177.                 return -1;
     
178.         if (!fprintf(File, (const char*)"\tBufferR[B] = BufferR[A] - TempR;\\\n"))
     
179.                 return -1;
     
180.         if (!fprintf(File, (const char*)"\tBufferV[B] = BufferV[A] - TempV;\\\n"))
     
181.                 return -1;
     
182.         if (!fprintf(File, (const char*)"\tBufferR[A] += TempR;\\\n"))
     
183.                 return -1;
     
184.         if (!fprintf(File, (const char*)"\tBufferV[A] += TempV;\\\n}\n"))
     
185.                 return -1;
     
186.         if (!fprintf(File, (const char*)"#define _Butters(BufferR,BufferV,TempR,TempV)\\\n{\\\n"))
     
187.                 return -1;
     
188.         for (k = 0; k < K; k++)
     
189.                 if (CreateButter(File, N, k))
     
190.                         return -1;
     
191.         if (!fprintf(File, (const char*)"}\n"))
     
192.                 return -1;
     
193.         return 0 == K - k ? 0: -1;
     
194. }
     
195. 
     
196. int8_ CreateDitFftFile(const char* FileName, uint32_ N)
     
197. {
     
198.         FILE* File = fopen(FileName, "w");
     
199.         if (!File)
     
200.                 return -1;
     
201.         if (CreateW(File, N))
     
202.         {
     
203.                 fclose(File);
     
204.                 return -1;
     
205.         }
     
206.         if (CreateReOrder(File, N))
     
207.         {
     
208.                 fclose(File);
     
209.                 return -1;
     
210.         }
     
211.         if (CreateButters(File, N))
     
212.         {
     
213.                 fclose(File);
     
214.                 return -1;
     
215.         }
     
216.         fclose(File);
     
217.         return 0;
     
218. }

复制代码

277955973

C fft PC MARK

zouzhichao

277955973 发表于 2013-8-13 12:51
C fft PC MARK

xjl00

赞一个！！

daidaimumu

不错，研究研究

zouzhichao

xjl00 发表于 2013-8-14 14:22
赞一个！！

承蒙谬赞

jxcylxh

非常好，顶一下。。

zouzhichao

jxcylxh 发表于 2013-8-14 19:29
非常好，顶一下。。

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

Jimmyxu

不错~~~~~~~~~

ddds

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

zouzhichao

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

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

ddds

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

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

zouzhichao

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

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

zouzhichao

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

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

zouzhichao

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

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

zouzhichao

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

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

zouzhichao

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

zouzhichao

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

zouzhichao

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

ddds

zouzhichao 发表于 2013-8-14 21:07
32点的DIT-FFT用了0.14S，4M晶振

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

zouzhichao

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

zouzhichao

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

guxingganyue

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

zouzhichao

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

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

zouzhichao

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

哈哈，又遇见你了

guxingganyue

zouzhichao 发表于 2013-8-14 21:48
哈哈，又遇见你了

呵呵，你回复我两次

chinaboy25

make,make,,

cqfeiyu

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

xjmlfm1

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

Pjm2008

不错。有时间再看看

robyn

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

zouzhichao

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

AVR Studio里跑的仿真

zouzhichao

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

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

zouzhichao

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

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

zouzhichao

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

多谢夸奖，不过实在惭愧

zouzhichao

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

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

xjmlfm1

zouzhichao 发表于 2013-8-15 12:47
不知道你们用的是整型运算还是浮点运算？

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

xjmlfm1

zouzhichao 发表于 2013-8-15 12:39
AVR Studio里跑的仿真

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

wmm20031015

好东东，顶起来

wcm_e

mark 学习

zouzhichao

wmm20031015 发表于 2013-8-15 17:07
好东东，顶起来

谢谢帮顶

zouzhichao

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

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

rifjft

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

zouzhichao

rifjft 发表于 2013-8-15 19:48
弄完这个，头发掉了几根？

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

zouzhichao

rifjft 发表于 2013-8-15 19:48
弄完这个，头发掉了几根？

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

rifjft

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

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

已荐精

zouzhichao

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

如若楼主之 ...

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

zouzhichao

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

如若楼主之 ...

还有多谢推荐

pipi516

FFT,先MARK一下!

zhanyanqiang

好东西，，，，

zouzhichao

zhanyanqiang 发表于 2013-8-18 16:43
好东西，，，，

好东西自然要分享

lixianghappy

mark                                      

acai039033

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

siyeb

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

zouzhichao

siyeb 发表于 2013-8-19 15:17
lz回头再把校正加上就更好了。

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

zouzhichao

siyeb 发表于 2013-8-19 15:17
lz回头再把校正加上就更好了。

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

siyeb

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

所以楼主加油啊~

wangfei1956

好东西，收藏了！

zouzhichao

wangfei1956 发表于 2013-8-27 22:24
好东西，收藏了！

欢迎收藏

晓毕8

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

binghe167

不错，感谢楼主分享

404710520

谢谢楼主共享

rainy_T

这个果断得顶起啊

loves6036

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

zouzhichao

loves6036 发表于 2013-8-29 12:29
谢谢楼主分享，真是宝贝啊，收藏了

承蒙夸奖

fangchangqing

这个好。正需要

zouzhichao

fangchangqing 发表于 2013-8-30 12:09
这个好。正需要

尽管拿去

fangchangqing

嘿嘿

zyw19987

zouzhichao 发表于 2013-8-30 12:50
尽管拿去

牛！
录个视频讲解吧

tscyds

很不错啊，学习学习

竹风xu

灰常好！！！感谢

407043680

蝶形算法  厉害

kevin_ldw

非常不错，学习了

touch_mcu

马克一下。

郭本宏

搞得不错，好好

qinkaiabc

收藏备用！

landmuto

fft，mark！

zouzhichao

landmuto 发表于 2013-9-14 19:43
fft，mark！

欢迎MARK

armok

COOL !

kmani

FFT 好东西！

zouzhichao

好意外！！！居然置酷了，Mygod

szmini2006

学习了。。。

moen

好帖留名

TongIC

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

xiaowenshao

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