OLD | NEW |
---|---|
1 /* | 1 /* |
2 * http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html | 2 * http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html |
3 * Copyright Takuya OOURA, 1996-2001 | 3 * Copyright Takuya OOURA, 1996-2001 |
4 * | 4 * |
5 * You may use, copy, modify and distribute this code for any purpose (include | 5 * You may use, copy, modify and distribute this code for any purpose (include |
6 * commercial use) and without fee. Please refer to this package when you modify | 6 * commercial use) and without fee. Please refer to this package when you modify |
7 * this code. | 7 * this code. |
8 * | 8 * |
9 * Changes: | 9 * Changes: |
10 * Trivial type modifications by the WebRTC authors. | 10 * Trivial type modifications by the WebRTC authors. |
(...skipping 268 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
279 a[j] *= 2.0 / n; | 279 a[j] *= 2.0 / n; |
280 } | 280 } |
281 . | 281 . |
282 | 282 |
283 | 283 |
284 Appendix : | 284 Appendix : |
285 The cos/sin table is recalculated when the larger table required. | 285 The cos/sin table is recalculated when the larger table required. |
286 w[] and ip[] are compatible with all routines. | 286 w[] and ip[] are compatible with all routines. |
287 */ | 287 */ |
288 | 288 |
289 static void makewt(int nw, int *ip, float *w); | 289 #include <stddef.h> |
290 static void makect(int nc, int *ip, float *c); | 290 |
291 static void bitrv2(int n, int *ip, float *a); | 291 static void makewt(size_t nw, size_t *ip, float *w); |
Andrew MacDonald
2015/07/24 19:01:55
It's not totally clear to me that this "work area
Peter Kasting
2015/07/27 23:23:39
The type of *ip has to match the type of |n| in th
Andrew MacDonald
2015/07/28 07:20:05
OK thanks for the explanation.
| |
292 static void makect(size_t nc, size_t *ip, float *c); | |
293 static void bitrv2(size_t n, size_t *ip, float *a); | |
292 #if 0 // Not used. | 294 #if 0 // Not used. |
293 static void bitrv2conj(int n, int *ip, float *a); | 295 static void bitrv2conj(int n, int *ip, float *a); |
294 #endif | 296 #endif |
295 static void cftfsub(int n, float *a, float *w); | 297 static void cftfsub(size_t n, float *a, float *w); |
296 static void cftbsub(int n, float *a, float *w); | 298 static void cftbsub(size_t n, float *a, float *w); |
297 static void cft1st(int n, float *a, float *w); | 299 static void cft1st(size_t n, float *a, float *w); |
298 static void cftmdl(int n, int l, float *a, float *w); | 300 static void cftmdl(size_t n, size_t l, float *a, float *w); |
299 static void rftfsub(int n, float *a, int nc, float *c); | 301 static void rftfsub(size_t n, float *a, size_t nc, float *c); |
300 static void rftbsub(int n, float *a, int nc, float *c); | 302 static void rftbsub(size_t n, float *a, size_t nc, float *c); |
301 #if 0 // Not used. | 303 #if 0 // Not used. |
302 static void dctsub(int n, float *a, int nc, float *c) | 304 static void dctsub(int n, float *a, int nc, float *c) |
303 static void dstsub(int n, float *a, int nc, float *c) | 305 static void dstsub(int n, float *a, int nc, float *c) |
304 #endif | 306 #endif |
305 | 307 |
306 | 308 |
307 #if 0 // Not used. | 309 #if 0 // Not used. |
308 void WebRtc_cdft(int n, int isgn, float *a, int *ip, float *w) | 310 void WebRtc_cdft(int n, int isgn, float *a, int *ip, float *w) |
309 { | 311 { |
310 if (n > (ip[0] << 2)) { | 312 if (n > (ip[0] << 2)) { |
311 makewt(n >> 2, ip, w); | 313 makewt(n >> 2, ip, w); |
312 } | 314 } |
313 if (n > 4) { | 315 if (n > 4) { |
314 if (isgn >= 0) { | 316 if (isgn >= 0) { |
315 bitrv2(n, ip + 2, a); | 317 bitrv2(n, ip + 2, a); |
316 cftfsub(n, a, w); | 318 cftfsub(n, a, w); |
317 } else { | 319 } else { |
318 bitrv2conj(n, ip + 2, a); | 320 bitrv2conj(n, ip + 2, a); |
319 cftbsub(n, a, w); | 321 cftbsub(n, a, w); |
320 } | 322 } |
321 } else if (n == 4) { | 323 } else if (n == 4) { |
322 cftfsub(n, a, w); | 324 cftfsub(n, a, w); |
323 } | 325 } |
324 } | 326 } |
325 #endif | 327 #endif |
326 | 328 |
327 | 329 |
328 void WebRtc_rdft(int n, int isgn, float *a, int *ip, float *w) | 330 void WebRtc_rdft(size_t n, int isgn, float *a, size_t *ip, float *w) |
329 { | 331 { |
330 int nw, nc; | 332 size_t nw, nc; |
331 float xi; | 333 float xi; |
332 | 334 |
333 nw = ip[0]; | 335 nw = ip[0]; |
334 if (n > (nw << 2)) { | 336 if (n > (nw << 2)) { |
335 nw = n >> 2; | 337 nw = n >> 2; |
336 makewt(nw, ip, w); | 338 makewt(nw, ip, w); |
337 } | 339 } |
338 nc = ip[1]; | 340 nc = ip[1]; |
339 if (n > (nc << 2)) { | 341 if (n > (nc << 2)) { |
340 nc = n >> 2; | 342 nc = n >> 2; |
(...skipping 295 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
636 a[0] = 0; | 638 a[0] = 0; |
637 } | 639 } |
638 #endif // Not used. | 640 #endif // Not used. |
639 | 641 |
640 | 642 |
641 /* -------- initializing routines -------- */ | 643 /* -------- initializing routines -------- */ |
642 | 644 |
643 | 645 |
644 #include <math.h> | 646 #include <math.h> |
645 | 647 |
646 static void makewt(int nw, int *ip, float *w) | 648 static void makewt(size_t nw, size_t *ip, float *w) |
647 { | 649 { |
648 int j, nwh; | 650 size_t j, nwh; |
649 float delta, x, y; | 651 float delta, x, y; |
650 | 652 |
651 ip[0] = nw; | 653 ip[0] = nw; |
652 ip[1] = 1; | 654 ip[1] = 1; |
653 if (nw > 2) { | 655 if (nw > 2) { |
654 nwh = nw >> 1; | 656 nwh = nw >> 1; |
655 delta = atanf(1.0f) / nwh; | 657 delta = atanf(1.0f) / nwh; |
656 w[0] = 1; | 658 w[0] = 1; |
657 w[1] = 0; | 659 w[1] = 0; |
658 w[nwh] = (float)cos(delta * nwh); | 660 w[nwh] = (float)cos(delta * nwh); |
659 w[nwh + 1] = w[nwh]; | 661 w[nwh + 1] = w[nwh]; |
660 if (nwh > 2) { | 662 if (nwh > 2) { |
661 for (j = 2; j < nwh; j += 2) { | 663 for (j = 2; j < nwh; j += 2) { |
662 x = (float)cos(delta * j); | 664 x = (float)cos(delta * j); |
663 y = (float)sin(delta * j); | 665 y = (float)sin(delta * j); |
664 w[j] = x; | 666 w[j] = x; |
665 w[j + 1] = y; | 667 w[j + 1] = y; |
666 w[nw - j] = y; | 668 w[nw - j] = y; |
667 w[nw - j + 1] = x; | 669 w[nw - j + 1] = x; |
668 } | 670 } |
669 bitrv2(nw, ip + 2, w); | 671 bitrv2(nw, ip + 2, w); |
670 } | 672 } |
671 } | 673 } |
672 } | 674 } |
673 | 675 |
674 | 676 |
675 static void makect(int nc, int *ip, float *c) | 677 static void makect(size_t nc, size_t *ip, float *c) |
676 { | 678 { |
677 int j, nch; | 679 size_t j, nch; |
678 float delta; | 680 float delta; |
679 | 681 |
680 ip[1] = nc; | 682 ip[1] = nc; |
681 if (nc > 1) { | 683 if (nc > 1) { |
682 nch = nc >> 1; | 684 nch = nc >> 1; |
683 delta = atanf(1.0f) / nch; | 685 delta = atanf(1.0f) / nch; |
684 c[0] = (float)cos(delta * nch); | 686 c[0] = (float)cos(delta * nch); |
685 c[nch] = 0.5f * c[0]; | 687 c[nch] = 0.5f * c[0]; |
686 for (j = 1; j < nch; j++) { | 688 for (j = 1; j < nch; j++) { |
687 c[j] = 0.5f * (float)cos(delta * j); | 689 c[j] = 0.5f * (float)cos(delta * j); |
688 c[nc - j] = 0.5f * (float)sin(delta * j); | 690 c[nc - j] = 0.5f * (float)sin(delta * j); |
689 } | 691 } |
690 } | 692 } |
691 } | 693 } |
692 | 694 |
693 | 695 |
694 /* -------- child routines -------- */ | 696 /* -------- child routines -------- */ |
695 | 697 |
696 | 698 |
697 static void bitrv2(int n, int *ip, float *a) | 699 static void bitrv2(size_t n, size_t *ip, float *a) |
698 { | 700 { |
699 int j, j1, k, k1, l, m, m2; | 701 size_t j, j1, k, k1, l, m, m2; |
700 float xr, xi, yr, yi; | 702 float xr, xi, yr, yi; |
701 | 703 |
702 ip[0] = 0; | 704 ip[0] = 0; |
703 l = n; | 705 l = n; |
704 m = 1; | 706 m = 1; |
705 while ((m << 3) < l) { | 707 while ((m << 3) < l) { |
706 l >>= 1; | 708 l >>= 1; |
707 for (j = 0; j < m; j++) { | 709 for (j = 0; j < m; j++) { |
708 ip[m + j] = ip[j] + l; | 710 ip[m + j] = ip[j] + l; |
709 } | 711 } |
(...skipping 186 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
896 a[k1 + 1] = xi; | 898 a[k1 + 1] = xi; |
897 } | 899 } |
898 k1 = 2 * k + ip[k]; | 900 k1 = 2 * k + ip[k]; |
899 a[k1 + 1] = -a[k1 + 1]; | 901 a[k1 + 1] = -a[k1 + 1]; |
900 a[k1 + m2 + 1] = -a[k1 + m2 + 1]; | 902 a[k1 + m2 + 1] = -a[k1 + m2 + 1]; |
901 } | 903 } |
902 } | 904 } |
903 } | 905 } |
904 #endif | 906 #endif |
905 | 907 |
906 static void cftfsub(int n, float *a, float *w) | 908 static void cftfsub(size_t n, float *a, float *w) |
907 { | 909 { |
908 int j, j1, j2, j3, l; | 910 size_t j, j1, j2, j3, l; |
909 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; | 911 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
910 | 912 |
911 l = 2; | 913 l = 2; |
912 if (n > 8) { | 914 if (n > 8) { |
913 cft1st(n, a, w); | 915 cft1st(n, a, w); |
914 l = 8; | 916 l = 8; |
915 while ((l << 2) < n) { | 917 while ((l << 2) < n) { |
916 cftmdl(n, l, a, w); | 918 cftmdl(n, l, a, w); |
917 l <<= 2; | 919 l <<= 2; |
918 } | 920 } |
(...skipping 27 matching lines...) Expand all Loading... | |
946 x0i = a[j + 1] - a[j1 + 1]; | 948 x0i = a[j + 1] - a[j1 + 1]; |
947 a[j] += a[j1]; | 949 a[j] += a[j1]; |
948 a[j + 1] += a[j1 + 1]; | 950 a[j + 1] += a[j1 + 1]; |
949 a[j1] = x0r; | 951 a[j1] = x0r; |
950 a[j1 + 1] = x0i; | 952 a[j1 + 1] = x0i; |
951 } | 953 } |
952 } | 954 } |
953 } | 955 } |
954 | 956 |
955 | 957 |
956 static void cftbsub(int n, float *a, float *w) | 958 static void cftbsub(size_t n, float *a, float *w) |
957 { | 959 { |
958 int j, j1, j2, j3, l; | 960 size_t j, j1, j2, j3, l; |
959 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; | 961 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
960 | 962 |
961 l = 2; | 963 l = 2; |
962 if (n > 8) { | 964 if (n > 8) { |
963 cft1st(n, a, w); | 965 cft1st(n, a, w); |
964 l = 8; | 966 l = 8; |
965 while ((l << 2) < n) { | 967 while ((l << 2) < n) { |
966 cftmdl(n, l, a, w); | 968 cftmdl(n, l, a, w); |
967 l <<= 2; | 969 l <<= 2; |
968 } | 970 } |
(...skipping 27 matching lines...) Expand all Loading... | |
996 x0i = -a[j + 1] + a[j1 + 1]; | 998 x0i = -a[j + 1] + a[j1 + 1]; |
997 a[j] += a[j1]; | 999 a[j] += a[j1]; |
998 a[j + 1] = -a[j + 1] - a[j1 + 1]; | 1000 a[j + 1] = -a[j + 1] - a[j1 + 1]; |
999 a[j1] = x0r; | 1001 a[j1] = x0r; |
1000 a[j1 + 1] = x0i; | 1002 a[j1 + 1] = x0i; |
1001 } | 1003 } |
1002 } | 1004 } |
1003 } | 1005 } |
1004 | 1006 |
1005 | 1007 |
1006 static void cft1st(int n, float *a, float *w) | 1008 static void cft1st(size_t n, float *a, float *w) |
1007 { | 1009 { |
1008 int j, k1, k2; | 1010 size_t j, k1, k2; |
1009 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; | 1011 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
1010 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; | 1012 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
1011 | 1013 |
1012 x0r = a[0] + a[2]; | 1014 x0r = a[0] + a[2]; |
1013 x0i = a[1] + a[3]; | 1015 x0i = a[1] + a[3]; |
1014 x1r = a[0] - a[2]; | 1016 x1r = a[0] - a[2]; |
1015 x1i = a[1] - a[3]; | 1017 x1i = a[1] - a[3]; |
1016 x2r = a[4] + a[6]; | 1018 x2r = a[4] + a[6]; |
1017 x2i = a[5] + a[7]; | 1019 x2i = a[5] + a[7]; |
1018 x3r = a[4] - a[6]; | 1020 x3r = a[4] - a[6]; |
(...skipping 82 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
1101 a[j + 10] = wk1r * x0r - wk1i * x0i; | 1103 a[j + 10] = wk1r * x0r - wk1i * x0i; |
1102 a[j + 11] = wk1r * x0i + wk1i * x0r; | 1104 a[j + 11] = wk1r * x0i + wk1i * x0r; |
1103 x0r = x1r + x3i; | 1105 x0r = x1r + x3i; |
1104 x0i = x1i - x3r; | 1106 x0i = x1i - x3r; |
1105 a[j + 14] = wk3r * x0r - wk3i * x0i; | 1107 a[j + 14] = wk3r * x0r - wk3i * x0i; |
1106 a[j + 15] = wk3r * x0i + wk3i * x0r; | 1108 a[j + 15] = wk3r * x0i + wk3i * x0r; |
1107 } | 1109 } |
1108 } | 1110 } |
1109 | 1111 |
1110 | 1112 |
1111 static void cftmdl(int n, int l, float *a, float *w) | 1113 static void cftmdl(size_t n, size_t l, float *a, float *w) |
1112 { | 1114 { |
1113 int j, j1, j2, j3, k, k1, k2, m, m2; | 1115 size_t j, j1, j2, j3, k, k1, k2, m, m2; |
1114 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; | 1116 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
1115 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; | 1117 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
1116 | 1118 |
1117 m = l << 2; | 1119 m = l << 2; |
1118 for (j = 0; j < l; j += 2) { | 1120 for (j = 0; j < l; j += 2) { |
1119 j1 = j + l; | 1121 j1 = j + l; |
1120 j2 = j1 + l; | 1122 j2 = j1 + l; |
1121 j3 = j2 + l; | 1123 j3 = j2 + l; |
1122 x0r = a[j] + a[j1]; | 1124 x0r = a[j] + a[j1]; |
1123 x0i = a[j + 1] + a[j1 + 1]; | 1125 x0i = a[j + 1] + a[j1 + 1]; |
(...skipping 104 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
1228 a[j1 + 1] = wk1r * x0i + wk1i * x0r; | 1230 a[j1 + 1] = wk1r * x0i + wk1i * x0r; |
1229 x0r = x1r + x3i; | 1231 x0r = x1r + x3i; |
1230 x0i = x1i - x3r; | 1232 x0i = x1i - x3r; |
1231 a[j3] = wk3r * x0r - wk3i * x0i; | 1233 a[j3] = wk3r * x0r - wk3i * x0i; |
1232 a[j3 + 1] = wk3r * x0i + wk3i * x0r; | 1234 a[j3 + 1] = wk3r * x0i + wk3i * x0r; |
1233 } | 1235 } |
1234 } | 1236 } |
1235 } | 1237 } |
1236 | 1238 |
1237 | 1239 |
1238 static void rftfsub(int n, float *a, int nc, float *c) | 1240 static void rftfsub(size_t n, float *a, size_t nc, float *c) |
1239 { | 1241 { |
1240 int j, k, kk, ks, m; | 1242 size_t j, k, kk, ks, m; |
1241 float wkr, wki, xr, xi, yr, yi; | 1243 float wkr, wki, xr, xi, yr, yi; |
1242 | 1244 |
1243 m = n >> 1; | 1245 m = n >> 1; |
1244 ks = 2 * nc / m; | 1246 ks = 2 * nc / m; |
1245 kk = 0; | 1247 kk = 0; |
1246 for (j = 2; j < m; j += 2) { | 1248 for (j = 2; j < m; j += 2) { |
1247 k = n - j; | 1249 k = n - j; |
1248 kk += ks; | 1250 kk += ks; |
1249 wkr = 0.5f - c[nc - kk]; | 1251 wkr = 0.5f - c[nc - kk]; |
1250 wki = c[kk]; | 1252 wki = c[kk]; |
1251 xr = a[j] - a[k]; | 1253 xr = a[j] - a[k]; |
1252 xi = a[j + 1] + a[k + 1]; | 1254 xi = a[j + 1] + a[k + 1]; |
1253 yr = wkr * xr - wki * xi; | 1255 yr = wkr * xr - wki * xi; |
1254 yi = wkr * xi + wki * xr; | 1256 yi = wkr * xi + wki * xr; |
1255 a[j] -= yr; | 1257 a[j] -= yr; |
1256 a[j + 1] -= yi; | 1258 a[j + 1] -= yi; |
1257 a[k] += yr; | 1259 a[k] += yr; |
1258 a[k + 1] -= yi; | 1260 a[k + 1] -= yi; |
1259 } | 1261 } |
1260 } | 1262 } |
1261 | 1263 |
1262 | 1264 |
1263 static void rftbsub(int n, float *a, int nc, float *c) | 1265 static void rftbsub(size_t n, float *a, size_t nc, float *c) |
1264 { | 1266 { |
1265 int j, k, kk, ks, m; | 1267 size_t j, k, kk, ks, m; |
1266 float wkr, wki, xr, xi, yr, yi; | 1268 float wkr, wki, xr, xi, yr, yi; |
1267 | 1269 |
1268 a[1] = -a[1]; | 1270 a[1] = -a[1]; |
1269 m = n >> 1; | 1271 m = n >> 1; |
1270 ks = 2 * nc / m; | 1272 ks = 2 * nc / m; |
1271 kk = 0; | 1273 kk = 0; |
1272 for (j = 2; j < m; j += 2) { | 1274 for (j = 2; j < m; j += 2) { |
1273 k = n - j; | 1275 k = n - j; |
1274 kk += ks; | 1276 kk += ks; |
1275 wkr = 0.5f - c[nc - kk]; | 1277 wkr = 0.5f - c[nc - kk]; |
(...skipping 45 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... | |
1321 kk += ks; | 1323 kk += ks; |
1322 wkr = c[kk] - c[nc - kk]; | 1324 wkr = c[kk] - c[nc - kk]; |
1323 wki = c[kk] + c[nc - kk]; | 1325 wki = c[kk] + c[nc - kk]; |
1324 xr = wki * a[k] - wkr * a[j]; | 1326 xr = wki * a[k] - wkr * a[j]; |
1325 a[k] = wkr * a[k] + wki * a[j]; | 1327 a[k] = wkr * a[k] + wki * a[j]; |
1326 a[j] = xr; | 1328 a[j] = xr; |
1327 } | 1329 } |
1328 a[m] *= c[0]; | 1330 a[m] *= c[0]; |
1329 } | 1331 } |
1330 #endif // Not used. | 1332 #endif // Not used. |
OLD | NEW |