野火电子论坛

 找回密码
 注册

QQ登录

只需一步,快速开始

查看: 11502|回复: 6

再来一个fft的源码,免移植

[复制链接]
发表于 2020-5-21 12:39:02 | 显示全部楼层 |阅读模式
  1. #include <stdio.h>
  2. #include <stdint.h>
  3. #include <math.h>

  4. typedef struct {
  5.     float r;
  6.     float i;
  7. }complex;

  8. static void butter_fly(complex* a, complex* b, const complex* c) {
  9.     complex bc;
  10.     bc.r = b->r * c->r - b->i * c->i;
  11.     bc.i = b->r * c->i + b->i * c->r;
  12.     b->r = a->r - bc.r;
  13.     b->i = a->i - bc.i;
  14.     a->r += bc.r;
  15.     a->i += bc.i;
  16. }

  17. static uint32_t bits_reverse(uint32_t index, uint32_t bits) {
  18.     uint32_t left, right;
  19.     left = index << 16;
  20.     right = index >> 16;
  21.     index = left | right;
  22.     left = (index << 8) & 0xff00ff00;
  23.     right = (index >> 8) & 0x00ff00ff;
  24.     index = left | right;
  25.     left = (index << 4) & 0xf0f0f0f0;
  26.     right = (index >> 4) & 0x0f0f0f0f;
  27.     index = left | right;
  28.     left = (index << 2) & 0xcccccccc;
  29.     right = (index >> 2) & 0x33333333;
  30.     index = left | right;
  31.     left = (index << 1) & 0xaaaaaaaa;
  32.     right = (index >> 1) & 0x55555555;
  33.     index = left | right;
  34.     return index >> (32 - bits);
  35. }

  36. static void fft_k(complex* dat, const complex* w, uint32_t k, uint32_t k_all) {
  37.     uint32_t i, j;
  38.     complex* dat1;
  39.     k_all = 1 << (k_all - k - 1);
  40.     k = 1 << k;
  41.     dat1 = dat + k;
  42.     for (i = 0; i < k_all; i++) {
  43.         for (j = 0; j < k; j++) {
  44.             butter_fly(dat++, dat1++, w + j * k_all);
  45.         }
  46.         dat += k;
  47.         dat1 += k;
  48.     }
  49. }

  50. void fft_init(complex* w, uint32_t k) {
  51.     float beta = 0.0f, dbeta = 3.1415926535f / k ;
  52.     for ( ; k; k--) {
  53.         w->r = cosf(beta);
  54.         w->i = sinf(beta);
  55.         beta += dbeta;
  56.         w++;
  57.     }
  58. }

  59. void fft(complex* dat, const complex* w, uint32_t k) {
  60.     uint32_t i, j, n;
  61.     complex temp;
  62.     n = 1 << k;
  63.     for (i = 0; i < n; i++) {
  64.         j = bits_reverse(i, k);
  65.         if (i <= j) {
  66.             continue;
  67.         }
  68.         temp = dat[i];
  69.         dat[i] = dat[j];
  70.         dat[j] = temp;
  71.     }
  72.     for (i = 0; i < k; i++) {
  73.         fft_k(dat, w, i, k);
  74.     }
  75. }
复制代码



评分

参与人数 1火花 +1024 收起 理由
flyleaf + 1024 1024!

查看全部评分

回复

使用道具 举报

 楼主| 发表于 2020-5-21 12:39:59 | 显示全部楼层
调用例程
  1. #define K 4
  2. #define N (1 << K)
  3. static complex w[N / 2];
  4. static complex dat[N];

  5. int main() {
  6.     uint32_t i;
  7.     for (i = 0; i < N; i++) {
  8.         dat[i].r = 0.0f;
  9.         dat[i].i = 0.0f;
  10.     }
  11.     dat[0].r = 1.0f;
  12.     dat[1].r = 1.0f;
  13.     fft_init((complex*)w, N / 2);
  14.     fft((complex*)dat, (const complex*)w, K);
  15.     for (i = 0; i < N; i++) {
  16.         printf((const char*)"dat[%d] = %f + %f * i\n", i, dat[i].r, dat[i].i);
  17.     }
  18.     return 0;
  19. }
复制代码
回复 支持 反对

使用道具 举报

发表于 2020-5-22 01:17:03 | 显示全部楼层
赞一个,非常好!      
回复 支持 反对

使用道具 举报

发表于 2020-5-22 07:18:49 | 显示全部楼层
这是zhouzichao大神吗?
回复 支持 反对

使用道具 举报

发表于 2020-5-22 09:24:41 | 显示全部楼层
大师,村长和火锅 欠你1024
回复 支持 反对

使用道具 举报

发表于 2021-1-4 09:26:22 | 显示全部楼层
请问这是多少点数的fft
回复 支持 反对

使用道具 举报

发表于 2021-12-19 22:31:42 | 显示全部楼层
983462736 发表于 2021-1-4 09:26
请问这是多少点数的fft

#define K 4
#define N (1 << K)
上面是2^4=16点的fft, 改K就能改变展开点数了
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

联系站长|手机版|野火电子官网|野火淘宝店铺|野火电子论坛 ( 粤ICP备14069197号 ) 大学生ARM嵌入式2群

GMT+8, 2024-3-29 03:29 , Processed in 0.056959 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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