#include <cmath>
#include "FFTHelper.h"
#include "const.h"

void FFTHelper::bitrp(double* xreal, double* ximag, int n)
{
    // λ·´×ªÖû» Bit-reversal Permutation
    int i, j, a, b, p;
    for (i = 1, p = 0; i < n; i *= 2)
    {
        p++;
    }
    for (i = 0; i < n; i++)
    {
        a = i;
        b = 0;
        for (j = 0; j < p; j++)
        {
            b = b * 2 + a % 2;
            a = a / 2;
        }
        if (b > i)
        {
            double t = xreal[i];
            xreal[i] = xreal[b];
            xreal[b] = t;
            t = ximag[i];
            ximag[i] = ximag[b];
            ximag[b] = t;
        }
    }
}


int FFTHelper::FFT(double* xreal, double* ximag, int len)
{
    //nֵΪ2µÄN´Î·½
    int n = 2;
    while (n <= len)
    {
        n *= 2;
    }
    n /= 2;
    // ¿ìËÙ¸µÁ¢Ò¶±ä»»£¬½«¸´Êý x ±ä»»ºóÈÔ±£´æÔÚ x ÖУ¬xreal, ximag ·Ö±ðÊÇ x µÄʵ²¿ºÍÐ鲿
    double* wreal = new double[n / 2];
    double* wimag = new double[n / 2];
    double treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;

    bitrp(xreal, ximag, n);
    // ¼ÆËã 1 µÄÇ° n / 2 ¸ö n ´Î·½¸ùµÄ¹²éÊý W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = -2 * pi / n;
    treal = cos(arg);
    timag = sin(arg);
    wreal[0] = 1.0f;
    wimag[0] = 0.0f;
    for (j = 1; j < n / 2; j++)
    {
        wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
        wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    }
    for (m = 2; m <= n; m *= 2)
    {
        for (k = 0; k < n; k += m)
        {
            for (j = 0; j < m / 2; j++)
            {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // ÐýתÒò×Ó w µÄʵ²¿ÔÚ wreal [] ÖеÄϱêΪ t
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
    }
    delete[] wreal;
    delete[] wimag;
    return n;
}

int FFTHelper::IFFT(double* xreal, double* ximag,int len)
{
    //nֵΪ2µÄN´Î·½
    int n = 2;
    while (n <= len)
    {
        n *= 2;
    }
    n /= 2;
    // ¿ìËÙ¸µÁ¢Ò¶Äæ±ä»»
    double* wreal = new double[n / 2];
    double* wimag = new double[n / 2];
    double treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
    bitrp(xreal, ximag, n);
    // ¼ÆËã 1 µÄÇ° n / 2 ¸ö n ´Î·½¸ù Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = 2 * pi / n;
    treal = cos(arg);
    timag = sin(arg);
    wreal[0] = 1.0;
    wimag[0] = 0.0;
    for (j = 1; j < n / 2; j++)
    {
        wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
        wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    }
    for (m = 2; m <= n; m *= 2)
    {
        for (k = 0; k < n; k += m)
        {
            for (j = 0; j < m / 2; j++)
            {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // ÐýתÒò×Ó w µÄʵ²¿ÔÚ wreal [] ÖеÄϱêΪ t
                treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                ureal = xreal[index1];
                uimag = ximag[index1];
                xreal[index1] = ureal + treal;
                ximag[index1] = uimag + timag;
                xreal[index2] = ureal - treal;
                ximag[index2] = uimag - timag;
            }
        }
    }
    for (j = 0; j < n; j++)
    {
        xreal[j] /= n;
        ximag[j] /= n;
    }

    delete[] wreal;
    delete[] wimag;
    return n;
}