FFT实现

前端时间搞前端处理，遇到RealFFT算法，挺感兴趣的，网上资料太少，就自己实现了一遍。
考虑大部分情况下，FFT输入都是实信号，所以一般的实现方法，都是在虚部补0之后，转成虚信号，进行变换，RealFFT就是快速处理实信号的一类算法，这个以后再谈。先是最基本的实现。

void ComplexFFT(float *R, float *I, int invert)
{
    int n, nn, i, j, m, cnt, inc, k;
    float tmpR, tmpI, WR, WI, Ri, Ii, Rj, Ij, dR, dI;

    n = R[0], nn = n >> 1;
    for(j = 0, i = 0; i < n - 1; i++)   
    {  
        if(i < j)
        {  
            tmpR = R[j + 1], tmpI = I[j + 1];
            R[j + 1] = R[i + 1], I[j + 1] = I[i + 1];
            R[i + 1] = tmpR, I[i + 1] = tmpI;
        }
        m = n >> 1;
        while(j >= m)
        {
            j = j - m;
            m = m >> 1;
        } 
        j = j + m;
    }

    m = 1;
    // 1, 2, 4 级
    while(m < n)
    {
        /*
            m = 1: [1, 2], [3, 4], [5, 6], [7, 8] 4
            m = 2: [1, 3], [2, 4], [5, 7], [6, 8] 2
            m = 4: [1, 5], [2, 6], [3, 7], [4, 8] 1
         */
        //printf("M = %d\n", m);
        cnt = 0, inc = n / (m << 1);
        // inc: 4 2 1
        // m  : 1 2 4
        // W递增inc
        while(cnt < inc)
        {
            // m = 1: 1 3 5 7
            // m = 2: 1 5
            // m = 4: 1
            i = cnt * m * 2 + 1;
            // W[0, n]: inc
            // 计算m次 迭代inc次
            for(int t = 0; t < m; t++, i++)
            {
                j = i + m;
                k = t * inc;
                // printf("[%3d, %3d] W[%3d, %3d]\n", i, j, k, nn);
                k == 0 ? WR = 1.0, WI = 0.0: WR = cos(PI * k / nn), WI = sin(PI * k / nn);
                if(invert) WI = - WI;
                //(R[i], I[i]) = (Ri, Ii) + W * (Rj, Ij)
                //(R[j], I[j]) = (Ri, Ii) - W * (Rj, Ij)
                Rj = R[j], Ij = I[j], Ri = R[i], Ii = I[i];
                R[i] = Ri + WR * Rj - WI * Ij, I[i] = Ii + WR * Ij + WI * Rj;
                R[j] = Ri - WR * Rj + WI * Ij, I[j] = Ii - WR * Ij - WI * Rj;
            }
            cnt++;
        }
        m = m << 1;
    }

    if (invert)
        for (i = 1; i <= n; i++)
            R[i] = R[i] / n, I[i] = I[i] / n;
}

用的是最基本的基于时间抽取的算法，写的时候参照下图：
下面对一个窗函数进行16点采样，测试结果如下，一次是输入数据，正变换和逆变换的结果。

最后顺便用C++也写了一下，用了一下虚数库。

#include <iostream>
#include <vector>
#include <complex>
#include <algorithm>

using namespace std;

const float PI = 3.14159265;

void FFT(vector< complex<float> > &sig, bool invert) {

    int N = sig.size();

    for (int i = 0, j = 0; i < N - 1; i++) {
        if(i < j) {
            swap(sig[i], sig[j]);
        }
        int k = N >> 1;
        while(j >= k) {
            j -= k;
            k = k >> 1;
        }
        j += k;
    }

    vector< complex<float> > W(N / 2);
    for (int i = 0; i < W.size(); i++) {
        W[i] = polar(1.0f, -2 * PI * i / N);
        if (invert) W[i] = conj(W[i]);
    }

    for (int m = 1; m < N; m = m << 1) {
        for (int n = 0; n < N; n += (m << 1)) {
            int k = 0;
            for(int j = n; j < n + m; j++, k += N / (m << 1)) {
                complex<float> A = sig[j], B = sig[j + m];
                sig[j] = A + W[k] * B, sig[j + m] = A - W[k] * B;
            }
        }
    }

    if (invert) {
        for_each(sig.begin(), sig.end(), [=] (complex<float> &s) {s /= N;});
    }
}

int main() {
    vector< complex<float> > signal = {1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

    FFT(signal, 0);
    FFT(signal, 1);

    for_each(signal.begin(), signal.end(), [](complex<float> s) {cout << abs(s) << endl;});
    return 0;
}