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

fucks

1. #include <stdio.h>
   
2. #include <stdint.h>
   
3. #include <math.h>
   
4. 
   
5. typedef struct {
   
6.     float r;
   
7.     float i;
   
8. }complex;
   
9. 
   
10. static void butter_fly(complex* a, complex* b, const complex* c) {
    
11.     complex bc;
    
12.     bc.r = b->r * c->r - b->i * c->i;
    
13.     bc.i = b->r * c->i + b->i * c->r;
    
14.     b->r = a->r - bc.r;
    
15.     b->i = a->i - bc.i;
    
16.     a->r += bc.r;
    
17.     a->i += bc.i;
    
18. }
    
19. 
    
20. static uint32_t bits_reverse(uint32_t index, uint32_t bits) {
    
21.     uint32_t left, right;
    
22.     left = index << 16;
    
23.     right = index >> 16;
    
24.     index = left | right;
    
25.     left = (index << 8) & 0xff00ff00;
    
26.     right = (index >> 8) & 0x00ff00ff;
    
27.     index = left | right;
    
28.     left = (index << 4) & 0xf0f0f0f0;
    
29.     right = (index >> 4) & 0x0f0f0f0f;
    
30.     index = left | right;
    
31.     left = (index << 2) & 0xcccccccc;
    
32.     right = (index >> 2) & 0x33333333;
    
33.     index = left | right;
    
34.     left = (index << 1) & 0xaaaaaaaa;
    
35.     right = (index >> 1) & 0x55555555;
    
36.     index = left | right;
    
37.     return index >> (32 - bits);
    
38. }
    
39. 
    
40. static void fft_k(complex* dat, const complex* w, uint32_t k, uint32_t k_all) {
    
41.     uint32_t i, j;
    
42.     complex* dat1;
    
43.     k_all = 1 << (k_all - k - 1);
    
44.     k = 1 << k;
    
45.     dat1 = dat + k;
    
46.     for (i = 0; i < k_all; i++) {
    
47.         for (j = 0; j < k; j++) {
    
48.             butter_fly(dat++, dat1++, w + j * k_all);
    
49.         }
    
50.         dat += k;
    
51.         dat1 += k;
    
52.     }
    
53. }
    
54. 
    
55. void fft_init(complex* w, uint32_t k) {
    
56.     float beta = 0.0f, dbeta = 3.1415926535f / k ;
    
57.     for ( ; k; k--) {
    
58.         w->r = cosf(beta);
    
59.         w->i = sinf(beta);
    
60.         beta += dbeta;
    
61.         w++;
    
62.     }
    
63. }
    
64. 
    
65. void fft(complex* dat, const complex* w, uint32_t k) {
    
66.     uint32_t i, j, n;
    
67.     complex temp;
    
68.     n = 1 << k;
    
69.     for (i = 0; i < n; i++) {
    
70.         j = bits_reverse(i, k);
    
71.         if (i <= j) {
    
72.             continue;
    
73.         }
    
74.         temp = dat[i];
    
75.         dat[i] = dat[j];
    
76.         dat[j] = temp;
    
77.     }
    
78.     for (i = 0; i < k; i++) {
    
79.         fft_k(dat, w, i, k);
    
80.     }
    
81. }

复制代码

fucks

调用例程

1. #define K 4
   
2. #define N (1 << K)
   
3. static complex w[N / 2];
   
4. static complex dat[N];
   
5. 
   
6. int main() {
   
7.     uint32_t i;
   
8.     for (i = 0; i < N; i++) {
   
9.         dat[i].r = 0.0f;
   
10.         dat[i].i = 0.0f;
    
11.     }
    
12.     dat[0].r = 1.0f;
    
13.     dat[1].r = 1.0f;
    
14.     fft_init((complex*)w, N / 2);
    
15.     fft((complex*)dat, (const complex*)w, K);
    
16.     for (i = 0; i < N; i++) {
    
17.         printf((const char*)"dat[%d] = %f + %f * i\n", i, dat[i].r, dat[i].i);
    
18.     }
    
19.     return 0;
    
20. }

复制代码

cgbabc

赞一个，非常好！      

netman

这是zhouzichao大神吗？

a123123

大师，村长和火锅 欠你1024

983462736

请问这是多少点数的fft

CrazySuiJi

983462736 发表于 2021-1-4 09:26
请问这是多少点数的fft

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