Skip to content

Commit 68b6ab8

Browse files
redchenjsredchenjs
committed
fft: use real-valued sequence FFT
1 parent 132e752 commit 68b6ab8

File tree

2 files changed

+56
-45
lines changed

2 files changed

+56
-45
lines changed

‎main/inc/user/fft.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,9 @@
1717
#define BAND_FADE (2)
1818
#define BAND_DELAY (2)
1919

20-
#define FFT_N (512)
20+
#define FFT_N (256)
2121
#define FFT_OUT_N (64)
22-
#define FFT_BLOCK_SIZE (FFT_N * 4)
22+
#define FFT_BLOCK_SIZE (FFT_N * 8)
2323

2424
#define MIN(a, b) ((a) < (b) ? (a) : (b))
2525
#define MAX(a, b) ((a) > (b) ? (a) : (b))

‎main/src/user/fft.c‎

Lines changed: 54 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313
#include "user/fft.h"
1414

1515
static int bitrev[FFT_N] = {0.0};
16-
static float window[FFT_N] = {0.0};
16+
static float window[FFT_N * 2] = {0.0};
1717
static float complex data[FFT_N] = {0.0};
1818
static float complex root[FFT_N / 2] = {0.0};
1919

@@ -35,24 +35,26 @@ static int bit_reverse(int x)
3535

3636
static void compute_fft_tables(void)
3737
{
38+
// bit reversal
3839
for (int i = 0; i < FFT_N; i++) {
3940
bitrev[i] = bit_reverse(i);
4041
}
4142

42-
// Hanning window
43-
for (int i = 0; i < FFT_N; i++) {
44-
window[i] = 0.5 * (1 - cosf(i * (TWO_PI / FFT_N)));
43+
// hamming window
44+
for (int i = 0; i < FFT_N * 2; i++) {
45+
window[i] = 0.53836 - 0.46164 * cosf(i * TWO_PI / (FFT_N * 2 - 1));
4546
}
4647

48+
// twiddle factor
4749
for (int i = 0; i < FFT_N / 2; i++) {
48-
root[i] = cexpf(-i * (TWO_PI / FFT_N) * I);
50+
root[i] = cexpf(-i * TWO_PI / (FFT_N - 1) * I);
4951
}
5052
}
5153

5254
static void compute_log_xscale(float *xscale, int bands)
5355
{
5456
for (int i = 0; i <= bands; i++) {
55-
xscale[i] = powf(FFT_N / 2, (float)i / bands) - 0.5f;
57+
xscale[i] = powf(FFT_N, (float)i / bands) - 0.5;
5658
}
5759
}
5860

@@ -71,30 +73,30 @@ static float compute_freq_band(const float *freq, const float *xscale, int band,
7173
for (; a < b; a++) {
7274
n += freq[a];
7375
}
74-
if (b < FFT_N / 2) {
76+
if (b < FFT_N) {
7577
n += freq[b] * (xscale[band + 1] - b);
7678
}
7779
}
7880

79-
n *= (float)bands / 12;
81+
n *= bands / 12.0;
8082

8183
return 20 * log10f(n);
8284
}
8385

8486
void fft_compute_lin(uint16_t *data_out, uint16_t scale_factor, uint16_t max_val, uint16_t min_val)
8587
{
86-
data_out[0] += cabsf(data[0]) / FFT_N * (max_val * scale_factor / 16384.0);
87-
data_out[0] /= 2;
88+
float freq[FFT_N] = {0.0};
8889

89-
if (data_out[0] > max_val) {
90-
data_out[0] = max_val;
91-
} else if (data_out[0] < min_val) {
92-
data_out[0] = min_val;
90+
for (int i = 0; i < FFT_N / 2; i++) {
91+
freq[i * 2] = cabsf(data[i] + conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0);
92+
freq[i * 2 + 1] = cabsf(data[i] - conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0);
9393
}
9494

95-
for (int i = 1; i < FFT_OUT_N; i++) {
96-
data_out[i] += cabsf(data[(FFT_N / 2 / FFT_OUT_N) * i]) / FFT_N * 2 * (max_val * scale_factor / 16384.0);
97-
data_out[i] /= 2;
95+
freq[0] /= 2.0;
96+
97+
for (int i = 0; i < FFT_OUT_N; i++) {
98+
data_out[i] += freq[FFT_N / FFT_OUT_N * i] * (max_val / 32.0);
99+
data_out[i] /= 2.0;
98100

99101
if (data_out[i] > max_val) {
100102
data_out[i] = max_val;
@@ -106,18 +108,18 @@ void fft_compute_lin(uint16_t *data_out, uint16_t scale_factor, uint16_t max_val
106108

107109
void fft_compute_log(uint16_t *data_out, uint16_t scale_factor, uint16_t max_val, uint16_t min_val)
108110
{
109-
data_out[0] += 20 * log10f(1 + cabsf(data[0]) / FFT_N) * (max_val * scale_factor / 16384.0);
110-
data_out[0] /= 2;
111+
float freq[FFT_N] = {0.0};
111112

112-
if (data_out[0] > max_val) {
113-
data_out[0] = max_val;
114-
} else if (data_out[0] < min_val) {
115-
data_out[0] = min_val;
113+
for (int i = 0; i < FFT_N / 2; i++) {
114+
freq[i * 2] = cabsf(data[i] + conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0);
115+
freq[i * 2 + 1] = cabsf(data[i] - conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0);
116116
}
117117

118-
for (int i = 1; i < FFT_OUT_N; i++) {
119-
data_out[i] += 20 * log10f(1 + cabsf(data[(FFT_N / 2 / FFT_OUT_N) * i]) / FFT_N * 2) * (max_val * scale_factor / 16384.0);
120-
data_out[i] /= 2;
118+
freq[0] /= 2.0;
119+
120+
for (int i = 0; i < FFT_OUT_N; i++) {
121+
data_out[i] += 20 * log10f(1 + freq[FFT_N / FFT_OUT_N * i]) * (max_val / 32.0);
122+
data_out[i] /= 2.0;
121123

122124
if (data_out[i] > max_val) {
123125
data_out[i] = max_val;
@@ -129,14 +131,15 @@ void fft_compute_log(uint16_t *data_out, uint16_t scale_factor, uint16_t max_val
129131

130132
void fft_compute_bands(uint16_t *data_out, uint16_t scale_factor, uint16_t max_val, uint16_t min_val)
131133
{
132-
float freq[FFT_N / 2] = {0.0};
134+
float freq[FFT_N] = {0.0};
133135

134-
freq[0] = cabsf(data[0]) / FFT_N * (scale_factor / 2.0 / 65536.0);
135-
136-
for (int i = 1; i < FFT_N / 2; i++) {
137-
freq[i] = 2 * cabsf(data[i]) / FFT_N * (scale_factor / 2.0 / 65536.0);
136+
for (int i = 0; i < FFT_N / 2; i++) {
137+
freq[i * 2] = cabsf(data[i] + conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0 / FFT_N);
138+
freq[i * 2 + 1] = cabsf(data[i] - conjf(data[FFT_N - 1 - i])) / 4.0 / FFT_N * (scale_factor / 512.0 / FFT_N);
138139
}
139140

141+
freq[0] /= 2.0;
142+
140143
for (int i = 0; i < BAND_N; i++) {
141144
float x = (40 + compute_freq_band(freq, xscale, i, BAND_N)) * (max_val / 64.0);
142145

@@ -164,14 +167,14 @@ void fft_execute(void)
164167
{
165168
int half = 1;
166169

167-
// Cooley–Tukey Fast Fourier Transform, radix-2 case
170+
// fast fourier transform, radix-2 case
168171
for (int i = FFT_N >> 1; i > 0; i >>= 1) {
169172
for (int g = 0; g < FFT_N; g += half << 1) {
170173
for (int b = 0, r = 0; b < half; b++, r += i) {
171174
float complex even = data[g + b];
172175
float complex odd = data[g + b + half] * root[r];
173176

174-
data[g + b] = even + odd;
177+
data[g + b] = even + odd;
175178
data[g + b + half] = even - odd;
176179
}
177180
}
@@ -182,29 +185,37 @@ void fft_execute(void)
182185

183186
void fft_load_data(const uint8_t *data_in, fft_channel_t channel)
184187
{
185-
int16_t data_l = 0, data_r = 0;
186-
187188
switch (channel) {
188189
case FFT_CHANNEL_L:
189190
for (int i = 0; i < FFT_N; i++) {
190-
data_l = data_in[i * 4 + 1] << 8 | data_in[i * 4];
191+
int16_t data_re_l = data_in[i * 8 + 1] << 8 | data_in[i * 8];
192+
int16_t data_im_l = data_in[i * 8 + 5] << 8 | data_in[i * 8 + 4];
193+
float data_re = data_re_l * window[i * 2];
194+
float data_im = data_im_l * window[i * 2 + 1];
191195

192-
data[bitrev[i]] = (float)data_l * window[i];
196+
data[bitrev[i]] = data_re + data_im * I;
193197
}
194198
break;
195199
case FFT_CHANNEL_R:
196200
for (int i = 0; i < FFT_N; i++) {
197-
data_r = data_in[i * 4 + 3] << 8 | data_in[i * 4 + 2];
201+
int16_t data_re_r = data_in[i * 8 + 3] << 8 | data_in[i * 8 + 2];
202+
int16_t data_im_r = data_in[i * 8 + 7] << 8 | data_in[i * 8 + 6];
203+
float data_re = data_re_r * window[i * 2];
204+
float data_im = data_im_r * window[i * 2 + 1];
198205

199-
data[bitrev[i]] = (float)data_r * window[i];
206+
data[bitrev[i]] = data_re + data_im * I;
200207
}
201208
break;
202209
case FFT_CHANNEL_LR:
203210
for (int i = 0; i < FFT_N; i++) {
204-
data_l = data_in[i * 4 + 1] << 8 | data_in[i * 4];
205-
data_r = data_in[i * 4 + 3] << 8 | data_in[i * 4 + 2];
206-
207-
data[bitrev[i]] = (float)((data_l + data_r) / 2) * window[i];
211+
int16_t data_re_l = data_in[i * 8 + 1] << 8 | data_in[i * 8];
212+
int16_t data_re_r = data_in[i * 8 + 3] << 8 | data_in[i * 8 + 2];
213+
int16_t data_im_l = data_in[i * 8 + 5] << 8 | data_in[i * 8 + 4];
214+
int16_t data_im_r = data_in[i * 8 + 7] << 8 | data_in[i * 8 + 6];
215+
float data_re = (data_re_l + data_re_r) / 2.0 * window[i * 2];
216+
float data_im = (data_im_l + data_im_r) / 2.0 * window[i * 2 + 1];
217+
218+
data[bitrev[i]] = data_re + data_im * I;
208219
}
209220
break;
210221
default:

0 commit comments

Comments
 (0)