CodeSnippet
C, fft, 代码

代码超级精简的 fft 源码

Catalogue 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <stdio.h>
#include <stdint.h>
#include <math.h>

typedef struct {
    float r;
    float i;
}complex;

static void butter_fly(complex* a, complex* b, const complex* c) {
    complex bc;
    bc.r = b->r * c->r - b->i * c->i;
    bc.i = b->r * c->i + b->i * c->r;
    b->r = a->r - bc.r;
    b->i = a->i - bc.i;
    a->r += bc.r;
    a->i += bc.i;
}

static uint32_t bits_reverse(uint32_t index, uint32_t bits) {
    uint32_t left, right;
    left = index << 16;
    right = index >> 16;
    index = left | right;
    left = (index << 8) & 0xff00ff00;
    right = (index >> 8) & 0x00ff00ff;
    index = left | right;
    left = (index << 4) & 0xf0f0f0f0;
    right = (index >> 4) & 0x0f0f0f0f;
    index = left | right;
    left = (index << 2) & 0xc3c3c3c3;
    right = (index >> 2) & 0x3c3c3c3c;
    index = left | right;
    left = (index << 1) & 0xaaaaaaaa;
    right = (index >> 1) & 0x55555555;
    index = left | right;
    return index >> (32 - bits);
}

static void fft_k(complex* dat, const complex* w, uint32_t k, uint32_t k_all) {
    uint32_t i, j;
    complex* dat1;
    k_all = 1 << (k_all - k - 1);
    k = 1 << k;
    dat1 = dat + k;
    for (i = 0; i < k_all; i++) {
        for (j = 0; j < k; j++) {
            butter_fly(dat++, dat1++, w + j * k_all);
        }
        dat += k;
        dat1 += k;
    }
}

void fft_init(complex* w, uint32_t k) {
    float beta = 0.0f, dbeta = 3.1415926535f / k ;
    for ( ; k; k--) {
        w->r = cosf(beta);
        w->i = sinf(beta);
        beta += dbeta;
        w++;
    }
}

void fft(complex* dat, const complex* w, uint32_t k) {
    uint32_t i, j, n;
    complex temp;
    n = 1 << k;
    for (i = 0; i < n; i++) {
        j = bits_reverse(i, k);
        if (i <= j) {
            continue;
        }
        temp = dat[i];
        dat[i] = dat[j];
        dat[j] = temp;
    }
    for (i = 0; i < k; i++) {
        fft_k(dat, w, i, k);
    }
}

使用方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#define K 4
#define N (1 << K)
static complex w[N / 2];
static complex dat[N];

int main() {
    uint32_t i;
    for (i = 0; i < N; i++) {
        dat[i].r = 0.0f;
        dat[i].i = 0.0f;
    }
    dat[0].r = 1.0f;
    dat[1].r = 1.0f;
    fft_init((complex*)w, N / 2);
    fft((complex*)dat, (const complex*)w, K);
    for (i = 0; i < N; i++) {
        printf((const char*)"dat[%d] = %f + %f * i\n", i, dat[i].r, dat[i].i);
    }
    return 0;
}

输出结果
fft_output
对比 matlab 输出结果
matlab_output
运算速度

stm32f030 16 点 fft 400us
stm32f030 64 点 fft 3900us
stm32f030 128 点 fft 9200us
开优化后 128 点 fft 8800us

https://debugdump.com/t_1535.html