daily pastebin goal
58%
SHARE
TWEET

Untitled

a guest Jan 21st, 2019 99 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /**
  2. * @brief Fast, tiny, portable and generic floating point 2d FFT function
  3. *
  4. * Author:   Stefan Moebius (mail@stefanmoebius.de)
  5. *
  6. * Date: 2012-10-03
  7. *
  8. * Licence: Released to the PUBLIC DOMAIN
  9. *
  10. * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
  11. * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
  12. * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
  13. * PARTICULAR PURPOSE.
  14. */
  15.  
  16. #include <math.h>
  17.  
  18. // 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")
  19. void fft2d(float** real, float** img, int reverse, int ld_width, int ld_height) {
  20.     float d, c_r, c_i;
  21.     int i, j, width, height, sz;
  22.  
  23.     // determine number of complex samples
  24.     width = 1 << ld_width; // n = 2 ^ ld_width
  25.     height = 1 << ld_height; // n = 2 ^ ld_height
  26.  
  27.     if (reverse)
  28.         d = 0.70710678118654757f; // sqrtf(0.5f)
  29.     else {  
  30.         // scale at forward transformation (could be also done after radix-2 algorithm)
  31.         d = 1.0f / (float)(width * height);
  32.         int y = height;
  33.         while (y--) {
  34.             int x = width;
  35.             while(x--) {
  36.                 real[y][x] *= d; // multiplication instead of division (for speed reason)
  37.                 img[y][x]  *= d;
  38.             }
  39.         }
  40.         d = -0.70710678118654757f;  // -sqrtf(0.5f)
  41.     }
  42.  
  43.     //1.) Transform rows
  44.     // sort array elements by reversed bit order (0000, 1000, 0100, 1100, 0010, ...)
  45.     j = width >> 1; // set leftmost bit
  46.     for (i=1; i < width-1; i++) {
  47.         if (i < j) {
  48.             int y = height;
  49.             while (y--) {
  50.                 float tmp = real[y][i];
  51.                 real[y][i] = real[y][j];
  52.                 real[y][j] = tmp;
  53.                 tmp = img[y][i];  
  54.                 img[y][i] = img[y][j];     
  55.                 img[y][j] = tmp;
  56.             }
  57.         }    
  58.         j ^= width - (unsigned int)width / (unsigned int) (1 + (i ^ (i + 1) ) ); // j contains the reversed bit order
  59.     }
  60.     // radix-2 algorithm
  61.     c_r = -1.0f;
  62.     c_i =  0.0f;
  63.     sz = 1;
  64.     while (ld_width--) {
  65.         float u_r = 1.0f, u_i = 0.0f;
  66.         int half = sz;
  67.         sz <<= 1;          
  68.         for (j=0; j < half; j++) {
  69.             float tmp_r, tmp_i;
  70.             for (i=j; i < width; i+=sz) {
  71.                 int ih = i + half;
  72.                 int y = height;
  73.                 while (y--) {
  74.                     tmp_r    = u_r * real[y][ih] - u_i * img[y][ih];
  75.                     tmp_i    = u_r * img[y][ih]  + u_i * real[y][ih];
  76.                     real[y][ih] = real[y][i] - tmp_r;
  77.                     img[y][ih]  = img[y][i]  - tmp_i;
  78.                     real[y][i] += tmp_r;
  79.                     img[y][i]  += tmp_i;
  80.                 }
  81.             }              
  82.             tmp_r = u_r;
  83.             u_r   = u_r   * c_r - u_i * c_i;
  84.             u_i   = tmp_r * c_i + u_i * c_r;
  85.         }
  86.         c_i = sqrtf( 1.0f - c_r) * d;
  87.         c_r = sqrtf( 0.5f * (1.0f + c_r));
  88.     }
  89.  
  90.     //2.) Transform columns
  91.     // sort array elements by reversed bit order (0000, 1000, 0100, 1100, 0010, ...)
  92.     j = height >> 1; // set leftmost bit
  93.     for (i=1; i < height-1; i++) {
  94.         if (i < j) {
  95.             int x = width;
  96.             while (x--) {
  97.                 float tmp = real[i][x];
  98.                 real[i][x] = real[j][x];
  99.                 real[j][x] = tmp;
  100.                 tmp = img[i][x];  
  101.                 img[i][x] = img[j][x];     
  102.                 img[j][x] = tmp;
  103.             }
  104.         }    
  105.         j ^= height - (unsigned int)height / (unsigned int) (1 + (i ^ (i + 1) ) ); // j contains the reversed bit order
  106.     }
  107.     // radix-2 algorithm
  108.     c_r = -1.0f;
  109.     c_i =  0.0f;
  110.     sz = 1;
  111.     while (ld_height--) {
  112.         float u_r = 1.0f, u_i = 0.0f;
  113.         int half = sz;
  114.         sz <<= 1;          
  115.         for (j=0; j < half; j++) {
  116.             float tmp_r, tmp_i;
  117.             for (i=j; i < height; i+=sz) {
  118.                 int ih = i + half;
  119.                 int x = width;
  120.                 while (x--) {
  121.                     tmp_r    = u_r * real[ih][x] - u_i * img[ih][x];
  122.                     tmp_i    = u_r * img[ih][x]  + u_i * real[ih][x];
  123.                     real[ih][x] = real[i][x] - tmp_r;
  124.                     img[ih][x]  = img[i][x]  - tmp_i;
  125.                     real[i][x] += tmp_r;
  126.                     img[i][x]  += tmp_i;
  127.                 }
  128.             }              
  129.             tmp_r = u_r;
  130.             u_r   = u_r   * c_r - u_i * c_i;
  131.             u_i   = tmp_r * c_i + u_i * c_r;
  132.         }
  133.         c_i = sqrtf( 1.0f - c_r) * d;
  134.         c_r = sqrtf( 0.5f * (1.0f + c_r));
  135.     }
  136. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top