递归版的快速傅里叶转换 (FFT),C语言实现

Dear 丶 2022-12-01 03:57 104阅读 0赞

看了好多个版本的快速傅里叶转换 (FFT),看懂了理论部分,但是却无法把理论联系到代码。也许是我的水平太差了吧。网上看到的版本中大量使用位操作符<<和>>。开始的时候还真的反应不过来。这里,我不讲FFT的理论,只是实现一下递归版本的FFT实现。有网友说FFT是玄科学,我深表赞同。相同的代码,在A电脑里运行结果是正确的,用B电脑就错误。找了好久才发现只需要把两个数组申请很大的空间即可,但这是为什么呢?反正没人告诉我原因。下面是具体代码:

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. double PI = 3.1415926;
  5. //定义复数
  6. struct Complex{
  7. double x; //实部
  8. double y; //虚部
  9. };
  10. //网上大部分都是采用C++实现FFT,使用了操作符重载。
  11. //但是其它语言没有操作符重载,所以我们还是用函数的方法实现复数的加法、减法和乘法
  12. struct Complex add(struct Complex a, struct Complex b)//复数的加法
  13. {
  14. struct Complex c;
  15. c.x = a.x + b.x;
  16. c.y = a.y + b.y;
  17. return c;
  18. }
  19. //复数的减法,不能用minus,好像是与关键字冲突
  20. struct Complex myMinus(struct Complex a, struct Complex b)
  21. {
  22. struct Complex c;
  23. c.x = a.x - b.x;
  24. c.y = a.y - b.y;
  25. return c;
  26. }
  27. struct Complex multiply(struct Complex a, struct Complex b) //复数的乘法
  28. {
  29. struct Complex c;
  30. c.x = a.x * b.x - a.y * b.y;
  31. c.y = a.x * b.y + b.x * a.y;
  32. return c;
  33. }
  34. void fft(int n, Complex * arr, int type)
  35. {
  36. if(n == 1) //递归的出口
  37. {
  38. return;
  39. }
  40. int m = n / 2; //由于保证了n是足够大的偶数,所以m一定是n的一半
  41. struct Complex *arr1; //储存arr的偶数部分
  42. arr1 = (struct Complex *)(malloc(m * sizeof(struct Complex))); //动态申请内存
  43. struct Complex *arr2; //储存arr的奇数部分
  44. arr2 = (struct Complex *)(malloc(m * sizeof(struct Complex))); //动态申请内存
  45. for (int i = 0; i < m; i++)
  46. {
  47. arr1[i] = arr[i * 2]; //在arr1中储存arr的偶数部分
  48. arr2[i] = arr[i * 2 + 1]; //在arr2中储存arr的奇数部分
  49. }
  50. fft(m, arr1, type); //对偶数部分进行快速傅里叶变换
  51. fft(m, arr2, type); //对奇数部分进行快速傅里叶变换
  52. struct Complex wn; //把360度等分成n份,这时的复数。
  53. wn.x = cos(2.0 * PI / n); //复数的x
  54. wn.y = type * sin(2.0 * PI / n); //复数的y,注意逆运算的时候,type为-1
  55. struct Complex w; //x = 1,y = 0的复数
  56. w.x = 1;
  57. w.y = 0;
  58. //计算arr中的每一个值
  59. for (int i = 0; i < m; i++)
  60. {
  61. arr[i] = add(arr1[i], multiply(w, arr2[i]));
  62. arr[i + m] = myMinus(arr1[i], multiply(w, arr2[i]));
  63. w = multiply(w, wn);
  64. }
  65. }
  66. int main()
  67. {
  68. int n = 9; //数组a的长度为9
  69. int m = 5; //数组b的长度为5
  70. int complexSize = sizeof(struct Complex);
  71. //要申请很大的空间,不知道为什么
  72. //如果只乘以100或者更小的空间,就会出错。
  73. struct Complex * a = (struct Complex *)(malloc(1000 * complexSize));
  74. struct Complex * b = (struct Complex *)(malloc(1000 * complexSize));
  75. struct Complex * c = (struct Complex *)(malloc(1000 * complexSize));
  76. for(int i = 0; i < n; i++) //数组a的初始化
  77. {
  78. a[i].x = 1;
  79. a[i].y = 0;
  80. }
  81. for(int i = 0; i < m; i++) //数组b的初始化
  82. {
  83. b[i].x = 1;
  84. b[i].y = 0;
  85. }
  86. int lim = 1;
  87. //lim为1, 2, 4, 8, 16, 32....除了1之外,其余都是偶数
  88. //lim为1时,fft直接就递归返回了,所以实际上lim都是偶数
  89. while (lim < n + m - 1)
  90. {
  91. lim = lim * 2;
  92. }
  93. //即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换,
  94. //能使傅里叶分析中许多问题的处理得到简化。
  95. //F(g(x) * f(x)) = F(g(x)) * F(f(x))
  96. //其中F表示的是傅里叶变换
  97. fft(lim, a, 1); //快速傅里叶转换
  98. fft(lim, b, 1); //快速傅里叶转换
  99. for (int i = 0; i < lim; i++)
  100. {
  101. c[i] = multiply(a[i], b[i]); //把a和b中相乘的结果都放到c中
  102. }
  103. fft(lim, c, -1); //对c进行快速傅里叶转换的逆转换
  104. //取实数四舍五入,此时虚数部分应当为0或由于浮点误差接近0
  105. //理论上计算结果应该是虚部为0,但实际上虚部可能有一点误差
  106. //为了抵消这个误差,用下面的方法进行处理
  107. //最终显示c的内容
  108. for (int i = 0; i < n + m - 1; i++)
  109. {
  110. printf("%d, ", (int)(c[i].x / lim + 0.5));
  111. }
  112. return 0;
  113. }

运行结果如下:

在这里插入图片描述
如果main函数中,struct Complex * a、b、c申请的空间比较小 (100或更小),则会出错,不知道为什么。

在这里插入图片描述
下次再写迭代版本的FFT吧。

发表评论

表情:
评论列表 (有 0 条评论,104人围观)

还没有评论,来说两句吧...

相关阅读

    相关 转换(FT),C语言实现

    最近在看傅里叶转换(FT)和快速傅里叶转换(FFT),感觉很难,包括号称“小学生都能看懂的快速傅里叶转换”,也是一样晦涩难懂,当然,这不排除现在的小学生智力都很强,而我的智商已