• API
• FAQ
• Tools
• Archive
daily pastebin goal
1%
SHARE
TWEET

# Untitled

a guest Jan 21st, 2019 119 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.     }
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.     }
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.

Top