再来一个fft的源码,免移植
#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) & 0xcccccccc;
right = (index >> 2) & 0x33333333;
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;
dat = dat;
dat = temp;
}
for (i = 0; i < k; i++) {
fft_k(dat, w, i, k);
}
}
调用例程
#define K 4
#define N (1 << K)
static complex w;
static complex dat;
int main() {
uint32_t i;
for (i = 0; i < N; i++) {
dat.r = 0.0f;
dat.i = 0.0f;
}
dat.r = 1.0f;
dat.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.r, dat.i);
}
return 0;
} 赞一个,非常好! 这是zhouzichao大神吗? 大师,村长和火锅 欠你1024:D 请问这是多少点数的fft
983462736 发表于 2021-1-4 09:26
请问这是多少点数的fft
#define K 4
#define N (1 << K)
上面是2^4=16点的fft, 改K就能改变展开点数了
页:
[1]