參考原始碼: https://cnx.org/contents/qAa9OhlP@2.44:zmcmahhR@7/Decimation-in-time-DIT-Radix-2-FFT
範例修改: mfft.c
#include "math.h"
#include "stdio.h"
void mfft(double *real, double *imag, int nT) {
int m = 0;
for (int i = 0; i < 30; i++) if (nT == (1 << i)) {
m = i;
break;
}
if ( m < 1) return;
const int halfT = nT >> 1;
for (int i = 1, j = 0; i < nT - 1; i++) {// bit-reverse
int k = halfT;
while (j >= k) {
j -= k;
k >>= 1;
}
j += k;
if (i < j) {
double temp = real[i];
real[i] = real[j];
real[j] = temp;
temp = imag[i];
imag[i] = imag[j];
imag[j] = temp;
}
}
double PIx2 = - M_PI * 2;
for (int i = 0, nx2 = 1; i < m; i++) { // radix-2 DIT FFT
int k = nx2;
nx2 += nx2;
double delta = PIx2 / nx2;
double theta = 0.0;
for (int j = 0; j < k; j++) {
double c = cos(theta);
double s = sin(theta);
theta += delta;
for (int n = j; n < nT; n += nx2) {
int h = n + k;
double w1 = c * real[h] - s * imag[h];
double w2 = s * real[h] + c * imag[h];
real[h] = real[n] - w1;
imag[h] = imag[n] - w2;
real[n] += w1;
imag[n] += w2;
}
}
}
}
int main() {
double x[16];
int n = sizeof(x) / sizeof(x[0]);
double y[n];
for (int i = 0; i < n; i++) {
x[i] = i;
y[i] = 0;
}
mfft(x, y , n);
for (int i = 0; i < n; i++) printf("%.3f,\t %.3f\n", x[i], y[i]);
return 0;
}
編譯並執行: g++ mfft.c && ./a.out
後記: 改成 Javascript 版本 mfft.js :
mfft = (real, imag, nT) => {
let m = 0;
for (let i = 0; i < 30; i++) if (nT == (1 << i)) {
m = i;
break;
}
if ( m < 1) return;
const halfT = nT >> 1;
for (let i = 1, j = 0; i < nT - 1; i++) {// bit-reverse
let k = halfT;
while (j >= k) {
j -= k;
k >>= 1;
}
j += k;
if (i < j) {
let temp = real[i];
real[i] = real[j];
real[j] = temp;
temp = imag[i];
imag[i] = imag[j];
imag[j] = temp;
}
}
let PIx2 = - Math.PI * 2;
for (let i = 0, nx2 = 1; i < m; i++) { // radix-2 DIT FFT
let k = nx2;
nx2 += nx2;
let delta = PIx2 / nx2;
let theta = 0.0;
for (let j = 0; j < k; j++) {
let c = Math.cos(theta);
let s = Math.sin(theta);
theta += delta;
for (let n = j; n < nT; n += nx2) {
let h = n + k;
let w1 = c * real[h] - s * imag[h];
let w2 = s * real[h] + c * imag[h];
real[h] = real[n] - w1;
imag[h] = imag[n] - w2;
real[n] += w1;
imag[n] += w2;
}
}
}
};
let n = 16;
let x = new Float32Array(n);
let y = new Float32Array(n);
for (let i = 0; i < n; i++) {
x[i] = i;
y[i] = 0;
}
mfft(x, y , n);
for (let i = 0; i < n; i++) console.log(`${Math.floor(x[i]*1e3 + 0.5)/1e3},\t ${Math.floor(y[i]*1e3 + 0.5)/1e3}`);
