Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /**
- * @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));
- }
- }
Add Comment
Please, Sign In to add comment