【翻译】最简单的DFT算法的C代码

0xAA55

翻译原文：https://batchloaf.wordpress.com/2013/12/07/simple-dft-in-c/
原文作者：Ted Burke

我最近花了一些时间研究离散傅立叶变换（DFT），以及快速傅里叶变换（FFT）能进行更高效的DFT变换计算。

DFT定义如下。指定离散时间序列x［n］，其中n = 0,1,2，...，N-1


De Moivre的定理表明：


因此，我们可以将DFT求和重写为X［k］的实部和虚部的表达式，如下所示（假设x［n］为实数）。


其中：


且


作为我学习过程的一部分，我一直在尝试编写C和MATLAB / Octave中DFT和FFT的不同实现。这让我感觉很爽，因为我写了个非常简单的DFT暴力求解算法。

1. //
   
2. // dft.c - 简单暴力求解DFT
   
3. // 作者：Ted Burke
   
4. // 最后更新：7-12-2013
   
5. //
   
6. // 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
   
7. //    gcc dft.c -o dft.exe
   
8. //
   
9. // 运行：
   
10. //    dft.exe
    
11. // （Linux上则用“./dft”来运行）
    
12. //
    
13. 
    
14. #include <stdio.h>
    
15. #include <stdlib.h>
    
16. #include <math.h>
    
17. 
    
18. #define N 16
    
19. #define PI 3.1415926535897932384626433832795
    
20. #define PI2 (PI * 2)
    
21. 
    
22. int main()
    
23. {
    
24.         // 时域和频域数据的数组
    
25.         int n, k;             // 时域和频域的索引
    
26.         float x[N];           // 离散时间信号x
    
27.         float Xre[N], Xim[N]; // DFT变换后的x（实部和虚部）
    
28.         float P[N];           // x的功率图
    
29.          
    
30.         // 生成在(-1,+1)区间内的随机离散时间信号x
    
31.         srand(time(0));
    
32.         for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
    
33.          
    
34.         // 用暴力算法计算x的DFT变换
    
35.         for (k=0 ; k<N ; ++k)
    
36.         {
    
37.                 // X[k]的实部
    
38.                 Xre[k] = 0;
    
39.                 for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);
    
40.                  
    
41.                 // X[k]的虚部
    
42.                 Xim[k] = 0;
    
43.                 for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);
    
44.                  
    
45.                 // 第k个频率点的功率
    
46.                 P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
    
47.         }
    
48.          
    
49.         // 输出结果到MATLAB / Octave的M文件来绘制图表
    
50.         FILE *f = fopen("dftplots.m", "w");
    
51.         fprintf(f, "n = [0:%d];\n", N-1);
    
52.         fprintf(f, "x = [ ");
    
53.         for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
    
54.         fprintf(f, "];\n");
    
55.         fprintf(f, "Xre = [ ");
    
56.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
    
57.         fprintf(f, "];\n");
    
58.         fprintf(f, "Xim = [ ");
    
59.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
    
60.         fprintf(f, "];\n");
    
61.         fprintf(f, "P = [ ");
    
62.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
    
63.         fprintf(f, "];\n");
    
64.         fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
    
65.         fprintf(f, "xlim([0 %d])\n", N-1);
    
66.         fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
    
67.         fprintf(f, "xlim([0 %d])\n", N-1);
    
68.         fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
    
69.         fprintf(f, "xlim([0 %d])\n", N-1);
    
70.         fclose(f);
    
71.          
    
72.         // 正常退出
    
73.         return 0;
    
74. }

复制代码

用下图的方式来编译并运行上面的代码：

当你运行编译好的程序dft.exe后，它会生成一个名为dftplots.m的M文件，然后可以在MATLAB或Octave中运行来生成结果图。这是dft.exe生成的M文件：

1. n = [0:15];
   
2. x = [ 0.672957 -0.453061 -0.835088 0.980334 0.972232 0.640295 0.791619 -0.042803 0.282745 0.153629 0.939992 0.588169 0.189058 0.461301 -0.667901 -0.314791 ];
   
3. Xre = [ 4.358686 -2.627209 -2.558252 2.144204 1.888348 3.210599 2.147089 -1.166725 0.332542 -1.166852 2.146972 3.210663 1.888432 2.144283 -2.558007 -2.627225 ];
   
4. Xim = [ 0.000000 -0.959613 -0.352430 1.534066 0.408726 -0.478590 -0.390162 0.160532 0.000022 -0.160326 0.389984 0.478336 -0.408804 -1.534269 0.352723 0.959997 ];
   
5. P = [ 18.998148 7.823084 6.668859 6.950971 3.732913 10.536995 4.762218 1.387017 0.110584 1.387247 4.761577 10.537162 3.733295 6.951930 6.667812 7.823905 ];
   
6. subplot(3,1,1)
   
7. plot(n,x)
   
8. xlim([0 15])
   
9. subplot(3,1,2)
   
10. plot(n,Xre,n,Xim)
    
11. xlim([0 15])
    
12. subplot(3,1,3)
    
13. stem(n,P)
    
14. xlim([0 15])

复制代码

我在Octave中打开了dftplots.m，如下所示：

它产生了如下图表：

最上面的是原始随机离散时间序列x［n］。第二个图是X［k］的实部和虚部。注意X［k］的对称性，正如我们对实数x［n］所预想的那样。最后那个图是P［k］也就是x［n］的功率图。P［k］的每个值是对应的复数值X［k］的大小的平方。

作为上述示例的后续，我稍微修改了程序来进行更快的DFT变换，并转换更长（64个样本）的数据帧。我还稍微修改了随机时域信号生成过程以添加主频分量（在k = 5.7，即在区间5和6之间），用于测试算法的正确性。

这是修改过的C代码：

1. //
   
2. // dft2.c - 基本，但效率更高的DFT
   
3. // 作者：Ted Burke
   
4. // 最后更新：10-12-2013
   
5. //
   
6. // 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
   
7. //    gcc dft2.c -o dft2.exe
   
8. //
   
9. // 运行：
   
10. //    dft2.exe
    
11. // （Linux上则用“./dft2”来运行）
    
12. //
    
13. 
    
14. #include <stdio.h>
    
15. #include <stdlib.h>
    
16. #include <math.h>
    
17. 
    
18. // 假设N大于4并且是2的幂
    
19. #define N 64
    
20. #define PI 3.1415926535897932384626433832795
    
21. #define PI2 (PI * 2)
    
22. 
    
23. // 旋转因子（64次单位根）
    
24. // https://zh.wikipedia.org/wiki/%E6%97%8B%E8%BD%89%E5%9B%A0%E5%AD%90
    
25. // https://zh.wikipedia.org/wiki/%E5%8D%95%E4%BD%8D%E6%A0%B9
    
26. const float W[] = {
    
27.          1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
    
28.          0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
    
29.         -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
    
30.         -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
    
31.         -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
    
32.         -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
    
33.          0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
    
34.          0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
    
35. };
    
36. 
    
37. int main()
    
38. {
    
39.         // 时域和频域数据的数组
    
40.         int n, k;             // 时域和频域的索引
    
41.         float x[N];           // 离散时间信号x
    
42.         float Xre[N], Xim[N]; // DFT变换后的x（实部和虚部）
    
43.         float P[N];           // x的功率图
    
44.          
    
45.         // 生成在(-1,+1)区间内的随机离散时间信号x
    
46.         srand(time(0));
    
47.         for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);
    
48.          
    
49.         // Calculate DFT and power spectrum up to Nyquist frequency
    
50.         // 计算DFT和功率图，上至奈奎斯特频率（在不引入误差的情况下对信号进行采样的最小速率，这是信号中存在的最高频率的两倍）
    
51.         int to_sin = 3*N/4; // sin的索引偏移
    
52.         int a, b;
    
53.         for (k=0 ; k<=N/2 ; ++k)
    
54.         {
    
55.                 Xre[k] = 0; Xim[k] = 0;
    
56.                 a = 0; b = to_sin;
    
57.                 for (n=0 ; n<N ; ++n)
    
58.                 {
    
59.                         Xre[k] += x[n] * W[a%N];
    
60.                         Xim[k] -= x[n] * W[b%N];
    
61.                         a += k; b += k;
    
62.                 }
    
63.                 P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
    
64.         }
    
65.          
    
66.         // 输出结果到MATLAB / Octave的M文件来绘制图表
    
67.         FILE *f = fopen("dftplots.m", "w");
    
68.         fprintf(f, "n = [0:%d];\n", N-1);
    
69.         fprintf(f, "x = [ ");
    
70.         for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
    
71.         fprintf(f, "];\n");
    
72.         fprintf(f, "Xre = [ ");
    
73.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
    
74.         fprintf(f, "];\n");
    
75.         fprintf(f, "Xim = [ ");
    
76.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
    
77.         fprintf(f, "];\n");
    
78.         fprintf(f, "P = [ ");
    
79.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
    
80.         fprintf(f, "];\n");
    
81.         fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
    
82.         fprintf(f, "xlim([0 %d])\n", N-1);
    
83.         fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
    
84.         fprintf(f, "xlim([0 %d])\n", N/2);  
    
85.         fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
    
86.         fprintf(f, "xlim([0 %d])\n", N/2);
    
87.         fclose(f);
    
88.          
    
89.         // 正常退出
    
90.         return 0;
    
91. }

复制代码

（译者注：有关“旋转因子”请看下图）

当我在Octave中打开生成的M文件时，生成了以下图表：

引入随机信号的主要频率分量，让上图中产生了单个突出峰值，证明了算法的正确性。

顺便说一句，我在上面的例子中使用以下C程序（在PC而不是dsPIC上运行）生成全局旋转因子数组的值：

1. //
   
2. // twiddle.c - 打印旋转因子数组
   
3. // 作者：Ted Burke
   
4. // 最后更新：10-12-2013
   
5. //
   
6. // 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
   
7. //    gcc twiddle.c -o twiddle.exe
   
8. //
   
9. // 运行：
   
10. //    twiddle.exe
    
11. // （Linux上则用“./twiddle”来运行）
    
12. //
    
13. 
    
14. #include <stdio.h>
    
15. #include <math.h>
    
16. 
    
17. #define N 64
    
18. #define PI 3.1415926535897932384626433832795
    
19. #define PI2 (PI * 2)
    
20. 
    
21. int main()
    
22. {
    
23.     int n;
    
24.      
    
25.     for (n=0 ; n<N ; ++n)
    
26.     {
    
27.         printf("%8.5lf", cos(n*PI2/N));
    
28.         if (n<N-1) printf(",");
    
29.         if ((n+1)%8==0) printf("\n");
    
30.     }
    
31.     return 0;
    
32. }

复制代码

青山彦樹

说到快速傅里叶变换算法，我知道傅里叶这个热力学大佬其实是被热死的。那个年代大概没有中暑的概念吧。

0xAA55

根据DFT算法的实现，对于大量数据的DFT变换似乎是可以进行GPU加速的。毕竟它就是针对每个单独的数据做一个求和的运算。

顺带一提谷歌浏览器的拼音插件真是个好东西，可以让我在只有日文输入法的机器上用Chrome输入中文。