// C言語で高速フーリエ変換（FFT）
#include <stdio.h>
#include <malloc.h>        // malloc/callocといった動的メモリ確保のために必要
#include <math.h>        // 数学関数を使用するためのヘッダ

#define SIZE 256        // 信号サイズ
#define PI 3.14159265    // 円周率
#define varsize float

// ビットを左右反転した配列を返す
int* BitScrollArray(int arraySize)
{
    int i, j;
    int* reBitArray = (int*)calloc(arraySize, sizeof(int));
    int arraySizeHarf = arraySize >> 1;

    reBitArray[0] = 0;
    for (i = 1; i < arraySize; i <<= 1)
    {
        for (j = 0; j < i; j++)
            reBitArray[j + i] = reBitArray[j] + arraySizeHarf;
        arraySizeHarf >>= 1;
    }
    return reBitArray;
}

// FFT
void FFT(float* inputRe,
         float* inputIm,
         float* outputRe,
         float* outputIm,
         int bitSize)
{
    int i, j, stage, type;
    int dataSize = 1 << bitSize;
    int butterflyDistance;
    int numType;
    int butterflySize;
    int jp;
    int* reverseBitArray = BitScrollArray(dataSize);
    float wRe, wIm, uRe, uIm, tempRe, tempIm, tempWRe, tempWIm;

    // バタフライ演算のための置き換え
    for (i = 0; i < dataSize; i++)
    {
        outputRe[i] = inputRe[reverseBitArray[i]];
        outputIm[i] = inputIm[reverseBitArray[i]];
    }

    // バタフライ演算
    for (stage = 1; stage <= bitSize; stage++)
    {
        butterflyDistance = 1 << stage;
        numType = butterflyDistance >> 1;
        butterflySize = butterflyDistance >> 1;

        wRe = 1.0;
        wIm = 0.0;
        uRe = cos(PI / butterflySize);
        uIm = -sin(PI / butterflySize);

        for (type = 0; type < numType; type++)
        {
            for (j = type; j < dataSize; j += butterflyDistance)
            {
                jp = j + butterflySize;
                tempRe = outputRe[jp] * wRe - outputIm[jp] * wIm;
                tempIm = outputRe[jp] * wIm + outputIm[jp] * wRe;
                outputRe[jp] = outputRe[j] - tempRe;
                outputIm[jp] = outputIm[j] - tempIm;
                outputRe[j] += tempRe;
                outputIm[j] += tempIm;
            }
            tempWRe = wRe * uRe - wIm * uIm;
            tempWIm = wRe * uIm + wIm * uRe;
            wRe = tempWRe;
            wIm = tempWIm;
        }
    }
}

int main(int argc, int** argv)
{
    int i;

    float* data_re = (float*)calloc(SIZE, sizeof(float));
    float* data_im = (float*)calloc(SIZE, sizeof(float));
    float* re = (float*)calloc(SIZE, sizeof(float));
    float* im = (float*)calloc(SIZE, sizeof(float));

    // 初期化
    for(i = 0; i < SIZE; i++)
    {
        data_re[i] = i;
        data_im[i] = 0;
        re[i] = 0;
        im[i] = 0;
    }

    // フーリエ変換
    FFT(data_re, data_im, re, im, 8);

    // 出力
    for(i = 0; i < SIZE; i++)
    {
        printf("%lf %lf %lf\n", data_re[i], re[i], im[i]);
    }
    
    free(data_re);
    free(data_im);
    free(re);
    free(im);
}