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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
// Fast Fourier Transform (adapted from spectral_editor/dct.js)
// Radix-2 Cooley-Tukey algorithm
// Bit-reversal permutation (in-place)
function bitReversePermute(real, imag, N) {
let temp_bits = N;
let num_bits = 0;
while (temp_bits > 1) {
temp_bits >>= 1;
num_bits++;
}
for (let i = 0; i < N; ++i) {
let j = 0;
let temp = i;
for (let b = 0; b < num_bits; ++b) {
j = (j << 1) | (temp & 1);
temp >>= 1;
}
if (j > i) {
const tmp_real = real[i];
const tmp_imag = imag[i];
real[i] = real[j];
imag[i] = imag[j];
real[j] = tmp_real;
imag[j] = tmp_imag;
}
}
}
// In-place radix-2 FFT
// direction: +1 for forward, -1 for inverse
function fftRadix2(real, imag, N, direction) {
const PI = Math.PI;
for (let stage_size = 2; stage_size <= N; stage_size *= 2) {
const half_stage = stage_size / 2;
const angle = direction * 2.0 * PI / stage_size;
let wr = 1.0;
let wi = 0.0;
const wr_delta = Math.cos(angle);
const wi_delta = Math.sin(angle);
for (let k = 0; k < half_stage; ++k) {
for (let group_start = k; group_start < N; group_start += stage_size) {
const i = group_start;
const j = group_start + half_stage;
const temp_real = real[j] * wr - imag[j] * wi;
const temp_imag = real[j] * wi + imag[j] * wr;
real[j] = real[i] - temp_real;
imag[j] = imag[i] - temp_imag;
real[i] = real[i] + temp_real;
imag[i] = imag[i] + temp_imag;
}
const wr_old = wr;
wr = wr_old * wr_delta - wi * wi_delta;
wi = wr_old * wi_delta + wi * wr_delta;
}
}
}
// Forward FFT: Time domain → Frequency domain
function fftForward(real, imag, N) {
bitReversePermute(real, imag, N);
fftRadix2(real, imag, N, +1);
}
// Real FFT wrapper for MQ extraction
// Input: Float32Array (time-domain signal)
// Output: Float32Array (interleaved [re0, im0, re1, im1, ...])
function realFFT(signal) {
const N = signal.length;
// Must be power of 2
if ((N & (N - 1)) !== 0) {
throw new Error('FFT size must be power of 2');
}
const real = new Float32Array(N);
const imag = new Float32Array(N);
// Copy input to real part
for (let i = 0; i < N; ++i) {
real[i] = signal[i];
}
// Compute FFT
fftForward(real, imag, N);
// Interleave output
const spectrum = new Float32Array(N * 2);
for (let i = 0; i < N; ++i) {
spectrum[i * 2] = real[i];
spectrum[i * 2 + 1] = imag[i];
}
return spectrum;
}
|
