FFT例子

最新推荐文章于 2021-05-14 08:45:16 发布

原创最新推荐文章于 2021-05-14 08:45:16 发布 · 745 阅读

0 ·

CC 4.0 BY-SA版权

CV 专栏收录该内容

53 篇文章

订阅专栏

#include <iostream>
#include <cmath>
#include <complex>
#include <iterator>

using namespace std;

const double PI = 3.1415926536;

// 位反转置换
unsigned int bitReverse(unsigned int x, int log2n)
{
    int n = 0;
    int mask = 0x1;
    for (int i = 0; i < log2n; i++)  // log2n = 3
 { n <<= 1; // n = 0, 0, 0  | n = 0, 2, 4 | n = 0, 0, 2 | n = 0, 2, 6
        n |= (x & 1); // n= 0, 0, 0 x = 0, 0, 0 | n =1, 2, 4 x = 1, 0, 0 | n = 0, 1, 2 x = 2, 1, 0 | n =1, 3, 6  x = 3, 1, 0
        x >>= 1; // x = 0, 0, 0 | x = 0, 0, 0 | x = 1, 0, 0 | x = 1, 0, 0
    }
    return n; // n = 0, 4, 2, 6
}

template<class Iter_T>
void fft(Iter_T a, Iter_T b, int log2n)
{
    typedef typename iterator_traits<Iter_T>::value_type complex;
    const complex J(0, 1);
    int n = 1 << log2n; // log2n = 3, n = 8
    for (unsigned int i = 0; i < n; ++i)
        b[bitReverse(i, log2n)] = a[i];
    for (int s = 1; s <= log2n; ++s)
    {
        int m = 1 << s; // 2
        int m2 = m >> 1; // 1
        complex w(1, 0);
        complex wm = exp(-J * (PI / m2)); // 2*PI/m
        for (int j = 0; j < m2; ++j)
        {
            for (int k = j; k < n; k+=m)
            {
                complex t = w * b[k + m2];
                complex u = b[k];
                b[k] = u + t;
                b[k + m2] = u - t;
            }
            w *= wm;
        }
    }
}

int main(int argc, char** argv)
{
    typedef complex<double> cx;
    cx a[] = { cx(0, 0), cx(1, 1), cx(3, 3), cx(4, 4),
    cx(4, 4), cx(3, 3), cx(1, 1), cx(0, 0)};
    cx b[8];
    fft(a, b, 3);
    for (int i = 0; i < 8; ++i)
        cout << b[i] << "\n";
    system("pause");
    return 0;
}