DIT FFT 範例修改

參考原始碼: 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}`);