/**
* @brief Fast, tiny, portable and generic floating point 2d FFT function
*
* Author:   Stefan Moebius (mail@stefanmoebius.de)
*
* Date: 2012-10-03
*
* Licence: Released to the PUBLIC DOMAIN
*
* THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
* KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
* PARTICULAR PURPOSE.
*/

#include <math.h>

// a fast and generic 2d FFT algorithm for power of 2 sizes (2,4,8,16,...). No temporary array required (transformation is done "in-place")
void fft2d(float** real, float** img, int reverse, int ld_width, int ld_height) {
    float d, c_r, c_i;
    int i, j, width, height, sz;

    // determine number of complex samples
    width = 1 << ld_width; // n = 2 ^ ld_width
    height = 1 << ld_height; // n = 2 ^ ld_height

    if (reverse)
        d = 0.70710678118654757f; // sqrtf(0.5f)
    else {
        // scale at forward transformation (could be also done after radix-2 algorithm)
        d = 1.0f / (float)(width * height);
        int y = height;
        while (y--) {
            int x = width;
            while(x--) {
                real[y][x] *= d; // multiplication instead of division (for speed reason)
                img[y][x]  *= d;
            }
        }
        d = -0.70710678118654757f;  // -sqrtf(0.5f)
    }

    //1.) Transform rows
    // sort array elements by reversed bit order (0000, 1000, 0100, 1100, 0010, ...)
    j = width >> 1; // set leftmost bit
    for (i=1; i < width-1; i++) {
        if (i < j) {
            int y = height;
            while (y--) {
                float tmp = real[y][i];
                real[y][i] = real[y][j];
                real[y][j] = tmp;
                tmp = img[y][i];
                img[y][i] = img[y][j];
                img[y][j] = tmp;
            }
        }
        j ^= width - (unsigned int)width / (unsigned int) (1 + (i ^ (i + 1) ) ); // j contains the reversed bit order
    }
    // radix-2 algorithm
    c_r = -1.0f;
    c_i =  0.0f;
    sz = 1;
    while (ld_width--) {
        float u_r = 1.0f, u_i = 0.0f;
        int half = sz;
        sz <<= 1;
        for (j=0; j < half; j++) {
            float tmp_r, tmp_i;
            for (i=j; i < width; i+=sz) {
                int ih = i + half;
                int y = height;
                while (y--) {
                    tmp_r    = u_r * real[y][ih] - u_i * img[y][ih];
                    tmp_i    = u_r * img[y][ih]  + u_i * real[y][ih];
                    real[y][ih] = real[y][i] - tmp_r;
                    img[y][ih]  = img[y][i]  - tmp_i;
                    real[y][i] += tmp_r;
                    img[y][i]  += tmp_i;
                }
            }
            tmp_r = u_r;
            u_r   = u_r   * c_r - u_i * c_i;
            u_i   = tmp_r * c_i + u_i * c_r;
        }
        c_i = sqrtf( 1.0f - c_r) * d;
        c_r = sqrtf( 0.5f * (1.0f + c_r));
    }

    //2.) Transform columns
    // sort array elements by reversed bit order (0000, 1000, 0100, 1100, 0010, ...)
    j = height >> 1; // set leftmost bit
    for (i=1; i < height-1; i++) {
        if (i < j) {
            int x = width;
            while (x--) {
                float tmp = real[i][x];
                real[i][x] = real[j][x];
                real[j][x] = tmp;
                tmp = img[i][x];
                img[i][x] = img[j][x];
                img[j][x] = tmp;
            }
        }
        j ^= height - (unsigned int)height / (unsigned int) (1 + (i ^ (i + 1) ) ); // j contains the reversed bit order
    }
    // radix-2 algorithm
    c_r = -1.0f;
    c_i =  0.0f;
    sz = 1;
    while (ld_height--) {
        float u_r = 1.0f, u_i = 0.0f;
        int half = sz;
        sz <<= 1;
        for (j=0; j < half; j++) {
            float tmp_r, tmp_i;
            for (i=j; i < height; i+=sz) {
                int ih = i + half;
                int x = width;
                while (x--) {
                    tmp_r    = u_r * real[ih][x] - u_i * img[ih][x];
                    tmp_i    = u_r * img[ih][x]  + u_i * real[ih][x];
                    real[ih][x] = real[i][x] - tmp_r;
                    img[ih][x]  = img[i][x]  - tmp_i;
                    real[i][x] += tmp_r;
                    img[i][x]  += tmp_i;
                }
            }
            tmp_r = u_r;
            u_r   = u_r   * c_r - u_i * c_i;
            u_i   = tmp_r * c_i + u_i * c_r;
        }
        c_i = sqrtf( 1.0f - c_r) * d;
        c_r = sqrtf( 0.5f * (1.0f + c_r));
    }
}