找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 6914|回复: 2

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

[复制链接]

1109

主题

1649

回帖

7万

积分

用户组: 管理员

一只技术宅

UID
1
精华
244
威望
743 点
宅币
24180 个
贡献
46222 次
宅之契约
0 份
在线时间
2294 小时
注册时间
2014-1-26
发表于 2018-10-24 10:10:30 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有账号?立即注册→加入我们

×
翻译原文: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
1.png

De Moivre的定理表明:
2.png

因此,我们可以将DFT求和重写为X[k]的实部和虚部的表达式,如下所示(假设x[n]为实数)。
3.png

其中:
4.png


5.png

作为我学习过程的一部分,我一直在尝试编写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. #include <stdio.h>
  14. #include <stdlib.h>
  15. #include <math.h>

  16. #define N 16
  17. #define PI 3.1415926535897932384626433832795
  18. #define PI2 (PI * 2)

  19. int main()
  20. {
  21.         // 时域和频域数据的数组
  22.         int n, k;             // 时域和频域的索引
  23.         float x[N];           // 离散时间信号x
  24.         float Xre[N], Xim[N]; // DFT变换后的x(实部和虚部)
  25.         float P[N];           // x的功率图
  26.          
  27.         // 生成在(-1,+1)区间内的随机离散时间信号x
  28.         srand(time(0));
  29.         for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
  30.          
  31.         // 用暴力算法计算x的DFT变换
  32.         for (k=0 ; k<N ; ++k)
  33.         {
  34.                 // X[k]的实部
  35.                 Xre[k] = 0;
  36.                 for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);
  37.                  
  38.                 // X[k]的虚部
  39.                 Xim[k] = 0;
  40.                 for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);
  41.                  
  42.                 // 第k个频率点的功率
  43.                 P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
  44.         }
  45.          
  46.         // 输出结果到MATLAB / Octave的M文件来绘制图表
  47.         FILE *f = fopen("dftplots.m", "w");
  48.         fprintf(f, "n = [0:%d];\n", N-1);
  49.         fprintf(f, "x = [ ");
  50.         for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
  51.         fprintf(f, "];\n");
  52.         fprintf(f, "Xre = [ ");
  53.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
  54.         fprintf(f, "];\n");
  55.         fprintf(f, "Xim = [ ");
  56.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
  57.         fprintf(f, "];\n");
  58.         fprintf(f, "P = [ ");
  59.         for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
  60.         fprintf(f, "];\n");
  61.         fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
  62.         fprintf(f, "xlim([0 %d])\n", N-1);
  63.         fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
  64.         fprintf(f, "xlim([0 %d])\n", N-1);
  65.         fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
  66.         fprintf(f, "xlim([0 %d])\n", N-1);
  67.         fclose(f);
  68.          
  69.         // 正常退出
  70.         return 0;
  71. }
复制代码
用下图的方式来编译并运行上面的代码:

当你运行编译好的程序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. #include <stdio.h>
  14. #include <stdlib.h>
  15. #include <math.h>

  16. // 假设N大于4并且是2的幂
  17. #define N 64
  18. #define PI 3.1415926535897932384626433832795
  19. #define PI2 (PI * 2)

  20. // 旋转因子(64次单位根)
  21. // https://zh.wikipedia.org/wiki/%E6%97%8B%E8%BD%89%E5%9B%A0%E5%AD%90
  22. // https://zh.wikipedia.org/wiki/%E5%8D%95%E4%BD%8D%E6%A0%B9
  23. const float W[] = {
  24.          1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
  25.          0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
  26.         -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
  27.         -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
  28.         -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
  29.         -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
  30.          0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
  31.          0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
  32. };

  33. int main()
  34. {
  35.         // 时域和频域数据的数组
  36.         int n, k;             // 时域和频域的索引
  37.         float x[N];           // 离散时间信号x
  38.         float Xre[N], Xim[N]; // DFT变换后的x(实部和虚部)
  39.         float P[N];           // x的功率图
  40.          
  41.         // 生成在(-1,+1)区间内的随机离散时间信号x
  42.         srand(time(0));
  43.         for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);
  44.          
  45.         // Calculate DFT and power spectrum up to Nyquist frequency
  46.         // 计算DFT和功率图,上至奈奎斯特频率(在不引入误差的情况下对信号进行采样的最小速率,这是信号中存在的最高频率的两倍)
  47.         int to_sin = 3*N/4; // sin的索引偏移
  48.         int a, b;
  49.         for (k=0 ; k<=N/2 ; ++k)
  50.         {
  51.                 Xre[k] = 0; Xim[k] = 0;
  52.                 a = 0; b = to_sin;
  53.                 for (n=0 ; n<N ; ++n)
  54.                 {
  55.                         Xre[k] += x[n] * W[a%N];
  56.                         Xim[k] -= x[n] * W[b%N];
  57.                         a += k; b += k;
  58.                 }
  59.                 P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
  60.         }
  61.          
  62.         // 输出结果到MATLAB / Octave的M文件来绘制图表
  63.         FILE *f = fopen("dftplots.m", "w");
  64.         fprintf(f, "n = [0:%d];\n", N-1);
  65.         fprintf(f, "x = [ ");
  66.         for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
  67.         fprintf(f, "];\n");
  68.         fprintf(f, "Xre = [ ");
  69.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
  70.         fprintf(f, "];\n");
  71.         fprintf(f, "Xim = [ ");
  72.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
  73.         fprintf(f, "];\n");
  74.         fprintf(f, "P = [ ");
  75.         for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
  76.         fprintf(f, "];\n");
  77.         fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
  78.         fprintf(f, "xlim([0 %d])\n", N-1);
  79.         fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
  80.         fprintf(f, "xlim([0 %d])\n", N/2);  
  81.         fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
  82.         fprintf(f, "xlim([0 %d])\n", N/2);
  83.         fclose(f);
  84.          
  85.         // 正常退出
  86.         return 0;
  87. }
复制代码
(译者注:有关“旋转因子”请看下图)
TwiddleFactor%20Diagram.jpg
当我在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. #include <stdio.h>
  14. #include <math.h>

  15. #define N 64
  16. #define PI 3.1415926535897932384626433832795
  17. #define PI2 (PI * 2)

  18. int main()
  19. {
  20.     int n;
  21.      
  22.     for (n=0 ; n<N ; ++n)
  23.     {
  24.         printf("%8.5lf", cos(n*PI2/N));
  25.         if (n<N-1) printf(",");
  26.         if ((n+1)%8==0) printf("\n");
  27.     }
  28.     return 0;
  29. }
复制代码

本帖被以下淘专辑推荐:

回复

使用道具 举报

0

主题

3

回帖

27

积分

用户组: 初·技术宅

UID
4330
精华
0
威望
0 点
宅币
24 个
贡献
0 次
宅之契约
0 份
在线时间
2 小时
注册时间
2018-9-29
发表于 2018-10-24 18:34:16 | 显示全部楼层
说到快速傅里叶变换算法,我知道傅里叶这个热力学大佬其实是被热死的。那个年代大概没有中暑的概念吧。
回复 赞! 靠!

使用道具 举报

1109

主题

1649

回帖

7万

积分

用户组: 管理员

一只技术宅

UID
1
精华
244
威望
743 点
宅币
24180 个
贡献
46222 次
宅之契约
0 份
在线时间
2294 小时
注册时间
2014-1-26
 楼主| 发表于 2018-10-25 09:14:42 | 显示全部楼层
根据DFT算法的实现,对于大量数据的DFT变换似乎是可以进行GPU加速的。毕竟它就是针对每个单独的数据做一个求和的运算。

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

使用道具 举报

QQ|Archiver|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )|网站地图

GMT+8, 2024-3-19 12:49 , Processed in 0.042138 second(s), 37 queries , Gzip On.

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表