免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
12下一页
最近访问板块 发新帖
查看: 3152 | 回复: 12

这个代码用g++编译OK,gcc编译出错,谁帮我改改? [复制链接]

论坛徽章:
0
发表于 2010-11-12 10:18 |显示全部楼层
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>

  4. const int N = 1024;
  5. const float PI = 3.1416;

  6. inline void swap(float &a, float &b)
  7. {
  8.   float t;
  9.   t = a;
  10.   a = b;
  11.   b = t;
  12. }

  13. void bitrp (float xreal[], float ximag[], int n)
  14. {
  15.   // 位反转置换 Bit-reversal Permutation
  16.   int i, j, a, b, p;

  17.   for (i = 1, p = 0; i < n; i *= 2) {
  18.     p++;
  19.   }
  20.   for (i = 0; i < n; i++) {
  21.     a = i;
  22.     b = 0;
  23.     for (j = 0; j < p; j++) {
  24.       b = (b << 1) + (a & 1);     // b = b * 2 + a % 2;
  25.       a >>= 1;         // a = a / 2;
  26.     }
  27.     if ( b > i) {
  28.       swap (xreal[i], xreal[b]);
  29.       swap (ximag[i], ximag[b]);
  30.     }
  31.   }
  32. }

  33. void FFT(float xreal [], float ximag [], int n)
  34. {
  35.   // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分别是 x 的实部和虚部
  36.   float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
  37.   int m, k, j, t, index1, index2;

  38.   bitrp (xreal, ximag, n);

  39.   // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal[j] + i * wimag[j] , j = 0, 1,   , n / 2 - 1
  40.   arg = - 2 * PI / n;
  41.   treal = cos (arg);
  42.   timag = sin (arg);
  43.   wreal [0] = 1.0;
  44.   wimag [0] = 0.0;
  45.   for (j = 1; j < n / 2; j++) {
  46.     wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
  47.     wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
  48.   }

  49.   for (m = 2; m <= n; m *= 2) {
  50.     for (k = 0; k < n; k += m) {
  51.       for (j = 0; j < m / 2; j++) {
  52.         index1 = k + j;
  53.         index2 = index1 + m / 2;
  54.         t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
  55.         treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
  56.         timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
  57.         ureal = xreal[index1];
  58.         uimag = ximag[index1];
  59.         xreal[index1] = ureal + treal;
  60.         ximag[index1] = uimag + timag;
  61.         xreal[index2] = ureal - treal;
  62.         ximag[index2] = uimag - timag;
  63.       }
  64.     }
  65.   }
  66. }

  67. void   IFFT (float xreal [], float ximag [], int n)
  68. {
  69.   // 快速傅立叶逆变换
  70.   float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
  71.   int m, k, j, t, index1, index2;

  72.   bitrp (xreal, ximag, n);

  73.   // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1,   , n / 2 - 1
  74.   arg = 2 * PI / n;
  75.   treal = cos (arg);
  76.   timag = sin (arg);
  77.   wreal [0] = 1.0;
  78.   wimag [0] = 0.0;
  79.   for (j = 1; j < n / 2; j++) {
  80.     wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
  81.     wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
  82.   }

  83.   for (m = 2; m <= n; m *= 2) {
  84.     for (k = 0; k < n; k += m) {
  85.       for (j = 0; j < m / 2; j++) {
  86.         index1 = k + j;
  87.         index2 = index1 + m / 2;
  88.         t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
  89.         treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
  90.         timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
  91.         ureal = xreal[index1];
  92.         uimag = ximag[index1];
  93.         xreal[index1] = ureal + treal;
  94.         ximag[index1] = uimag + timag;
  95.         xreal[index2] = ureal - treal;
  96.         ximag[index2] = uimag - timag;
  97.       }
  98.     }
  99.   }

  100.   for (j=0; j < n; j++) {
  101.     xreal[j] /= n;
  102.     ximag[j] /= n;
  103.   }
  104. }

  105. void FFT_test ()
  106. {
  107.   char inputfile[] = "input.txt";
  108.   char outputfile[] = "output.txt";
  109.   float xreal[N] = {}, ximag[N] = {};
  110.   int n, i;
  111.   FILE *input, *output;

  112.   if (!(input = fopen (inputfile, "r"))) {
  113.     printf ("Cannot open file.\n");
  114.     exit (1);
  115.   }
  116.   if (!(output = fopen (outputfile, "w"))) {
  117.     printf ("Cannot open file.\n");
  118.     exit (1);
  119.   }

  120.   i = 0;
  121.   while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF) {
  122.     i++;
  123.   }
  124.   n = i;     // 要求 n 为 2 的整数幂
  125.   while (i > 1) {
  126.     if (i % 2) {
  127.       fprintf (output, "%d is not a power of 2!\n", n);
  128.       exit (1);
  129.     }
  130.     i /= 2;
  131.   }

  132.   FFT (xreal, ximag, n);
  133.   fprintf (output, "FFT:     i      real imag\n");
  134.   for (i = 0; i < n; i++) {
  135.     fprintf (output, "%4d     %8.4f     %8.4f\n", i, xreal[i], ximag[i]);
  136.   }
  137.   fprintf (output, "=================================\n");

  138.   IFFT (xreal, ximag, n);
  139.   fprintf (output, "IFFT:     i      real imag\n");
  140.   for (i = 0; i < n; i++) {
  141.     fprintf (output, "%4d     %8.4f     %8.4f\n", i, xreal[i], ximag[i]);
  142.   }

  143.   if ( fclose (input)) {
  144.     printf ("File close error.\n");
  145.     exit (1);
  146.   }
  147.   if ( fclose (output)) {
  148.     printf ("File close error.\n");
  149.     exit (1);
  150.   }
  151. }

  152. int main (void)
  153. {
  154.   FFT_test ();

  155.   return 0;
  156. }
复制代码
该死的FFT

论坛徽章:
278
射手座
日期:2013-08-23 12:04:38射手座
日期:2013-08-23 16:18:12未羊
日期:2013-08-30 14:33:15水瓶座
日期:2013-09-02 16:44:31摩羯座
日期:2013-09-25 09:33:52双子座
日期:2013-09-26 12:21:10金牛座
日期:2013-10-14 09:08:49申猴
日期:2013-10-16 13:09:43子鼠
日期:2013-10-17 23:23:19射手座
日期:2013-10-18 13:00:27金牛座
日期:2013-10-18 15:47:57午马
日期:2013-10-18 21:43:38
发表于 2010-11-12 10:26 |显示全部楼层
const int 定义改define

float &引用改成指针

论坛徽章:
0
发表于 2010-11-12 10:29 |显示全部楼层
float &引用改成指针可以么?
传进去的参数就是数组的下标啊。我改了不行的。

论坛徽章:
278
射手座
日期:2013-08-23 12:04:38射手座
日期:2013-08-23 16:18:12未羊
日期:2013-08-30 14:33:15水瓶座
日期:2013-09-02 16:44:31摩羯座
日期:2013-09-25 09:33:52双子座
日期:2013-09-26 12:21:10金牛座
日期:2013-10-14 09:08:49申猴
日期:2013-10-16 13:09:43子鼠
日期:2013-10-17 23:23:19射手座
日期:2013-10-18 13:00:27金牛座
日期:2013-10-18 15:47:57午马
日期:2013-10-18 21:43:38
发表于 2010-11-12 10:48 |显示全部楼层
inline void swap(float *a, float *b)
{
  float t;

  t = *a;
  *a = *b;
  *b = t;
}

      swap (&xreal[i], &xreal[b]);

      swap (&ximag[i], &ximag[b]);

论坛徽章:
0
发表于 2010-11-12 10:52 |显示全部楼层
特别坏,你编译通过了么?

论坛徽章:
278
射手座
日期:2013-08-23 12:04:38射手座
日期:2013-08-23 16:18:12未羊
日期:2013-08-30 14:33:15水瓶座
日期:2013-09-02 16:44:31摩羯座
日期:2013-09-25 09:33:52双子座
日期:2013-09-26 12:21:10金牛座
日期:2013-10-14 09:08:49申猴
日期:2013-10-16 13:09:43子鼠
日期:2013-10-17 23:23:19射手座
日期:2013-10-18 13:00:27金牛座
日期:2013-10-18 15:47:57午马
日期:2013-10-18 21:43:38
发表于 2010-11-12 11:00 |显示全部楼层
特别坏,你编译通过了么?
prolj 发表于 2010-11-12 10:52



    编译没问题,链接错,加上 -lm就可以了


gcc fft.c -lm

论坛徽章:
0
发表于 2010-11-12 11:08 |显示全部楼层
回复 6# hellioncu


    数学运算库不是默认加载的。需要显示指定。

  这个大牛就是大牛哦

论坛徽章:
0
发表于 2010-11-12 11:14 |显示全部楼层
-lm这个我肯定知道啊。
可是我这里还是不对啊。

论坛徽章:
278
射手座
日期:2013-08-23 12:04:38射手座
日期:2013-08-23 16:18:12未羊
日期:2013-08-30 14:33:15水瓶座
日期:2013-09-02 16:44:31摩羯座
日期:2013-09-25 09:33:52双子座
日期:2013-09-26 12:21:10金牛座
日期:2013-10-14 09:08:49申猴
日期:2013-10-16 13:09:43子鼠
日期:2013-10-17 23:23:19射手座
日期:2013-10-18 13:00:27金牛座
日期:2013-10-18 15:47:57午马
日期:2013-10-18 21:43:38
发表于 2010-11-12 11:20 |显示全部楼层
回复 8# prolj
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>

  4. #define N        1024
  5. #define PI        3.1416

  6. inline void swap(float *a, float *b)
  7. {
  8.   float t;
  9.   t = *a;
  10.   *a = *b;
  11.   *b = t;
  12. }

  13. void bitrp (float xreal[], float ximag[], int n)
  14. {
  15.   // 位反转置换 Bit-reversal Permutation
  16.   int i, j, a, b, p;

  17.   for (i = 1, p = 0; i < n; i *= 2) {
  18.     p++;
  19.   }
  20.   for (i = 0; i < n; i++) {
  21.     a = i;
  22.     b = 0;
  23.     for (j = 0; j < p; j++) {
  24.       b = (b << 1) + (a & 1);     // b = b * 2 + a % 2;
  25.       a >>= 1;         // a = a / 2;
  26.     }
  27.     if ( b > i) {
  28.       swap (&xreal[i], &xreal[b]);
  29.       swap (&ximag[i], &ximag[b]);
  30.     }
  31.   }
  32. }

  33. void FFT(float xreal [], float ximag [], int n)
  34. {
  35.   // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分别是 x 的实部和虚部
  36.   float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
  37.   int m, k, j, t, index1, index2;

  38.   bitrp (xreal, ximag, n);

  39.   // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal[j] + i * wimag[j] , j = 0, 1,   , n / 2 - 1
  40.   arg = - 2 * PI / n;
  41.   treal = cos (arg);
  42.   timag = sin (arg);
  43.   wreal [0] = 1.0;
  44.   wimag [0] = 0.0;
  45.   for (j = 1; j < n / 2; j++) {
  46.     wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
  47.     wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
  48.   }

  49.   for (m = 2; m <= n; m *= 2) {
  50.     for (k = 0; k < n; k += m) {
  51.       for (j = 0; j < m / 2; j++) {
  52.         index1 = k + j;
  53.         index2 = index1 + m / 2;
  54.         t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
  55.         treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
  56.         timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
  57.         ureal = xreal[index1];
  58.         uimag = ximag[index1];
  59.         xreal[index1] = ureal + treal;
  60.         ximag[index1] = uimag + timag;
  61.         xreal[index2] = ureal - treal;
  62.         ximag[index2] = uimag - timag;
  63.       }
  64.     }
  65.   }
  66. }

  67. void   IFFT (float xreal [], float ximag [], int n)
  68. {
  69.   // 快速傅立叶逆变换
  70.   float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
  71.   int m, k, j, t, index1, index2;

  72.   bitrp (xreal, ximag, n);

  73.   // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1,   , n / 2 - 1
  74.   arg = 2 * PI / n;
  75.   treal = cos (arg);
  76.   timag = sin (arg);
  77.   wreal [0] = 1.0;
  78.   wimag [0] = 0.0;
  79.   for (j = 1; j < n / 2; j++) {
  80.     wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
  81.     wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
  82.   }

  83.   for (m = 2; m <= n; m *= 2) {
  84.     for (k = 0; k < n; k += m) {
  85.       for (j = 0; j < m / 2; j++) {
  86.         index1 = k + j;
  87.         index2 = index1 + m / 2;
  88.         t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
  89.         treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
  90.         timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
  91.         ureal = xreal[index1];
  92.         uimag = ximag[index1];
  93.         xreal[index1] = ureal + treal;
  94.         ximag[index1] = uimag + timag;
  95.         xreal[index2] = ureal - treal;
  96.         ximag[index2] = uimag - timag;
  97.       }
  98.     }
  99.   }

  100.   for (j=0; j < n; j++) {
  101.     xreal[j] /= n;
  102.     ximag[j] /= n;
  103.   }
  104. }

  105. void FFT_test ()
  106. {
  107.   char inputfile[] = "input.txt";
  108.   char outputfile[] = "output.txt";
  109.   float xreal[N] = {}, ximag[N] = {};
  110.   int n, i;
  111.   FILE *input, *output;

  112.   if (!(input = fopen (inputfile, "r"))) {
  113.     printf ("Cannot open file.\n");
  114.     exit (1);
  115.   }
  116.   if (!(output = fopen (outputfile, "w"))) {
  117.     printf ("Cannot open file.\n");
  118.     exit (1);
  119.   }

  120.   i = 0;
  121.   while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF) {
  122.     i++;
  123.   }
  124.   n = i;     // 要求 n 为 2 的整数幂
  125.   while (i > 1) {
  126.     if (i % 2) {
  127.       fprintf (output, "%d is not a power of 2!\n", n);
  128.       exit (1);
  129.     }
  130.     i /= 2;
  131.   }

  132.   FFT (xreal, ximag, n);
  133.   fprintf (output, "FFT:     i      real imag\n");
  134.   for (i = 0; i < n; i++) {
  135.     fprintf (output, "%4d     %8.4f     %8.4f\n", i, xreal[i], ximag[i]);
  136.   }
  137.   fprintf (output, "=================================\n");

  138.   IFFT (xreal, ximag, n);
  139.   fprintf (output, "IFFT:     i      real imag\n");
  140.   for (i = 0; i < n; i++) {
  141.     fprintf (output, "%4d     %8.4f     %8.4f\n", i, xreal[i], ximag[i]);
  142.   }

  143.   if ( fclose (input)) {
  144.     printf ("File close error.\n");
  145.     exit (1);
  146.   }
  147.   if ( fclose (output)) {
  148.     printf ("File close error.\n");
  149.     exit (1);
  150.   }
  151. }

  152. int main (void)
  153. {
  154.   FFT_test ();

  155.   return 0;
  156. }

复制代码

论坛徽章:
0
发表于 2010-11-12 11:23 |显示全部楼层
哦,不知道我哪里改错了。
您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

SACC2019中国系统架构师大会

【数字转型 架构演进】SACC2019中国系统架构师大会,8.5折限时优惠重磅来袭!
2019年10月31日~11月2日第11届中国系统架构师大会(SACC2019)将在北京隆重召开。四大主线并行的演讲模式,1个主会场、20个技术专场、超千人参与的会议规模,100+来自互联网、金融、制造业、电商等领域的嘉宾阵容,将为广大参会者提供一场最具价值的技术交流盛会。

限时8.5折扣期:2019年9月30日前


----------------------------------------

大会官网>>
  

北京盛拓优讯信息技术有限公司. 版权所有 16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122
中国互联网协会会员  联系我们:huangweiwei@it168.com
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP