| 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. |
| 11 */ | 11 */ |
| 12 | 12 |
| 13 /* | 13 /* |
| 14 Fast Fourier/Cosine/Sine Transform | 14 Fast Fourier/Cosine/Sine Transform |
| 15 dimension :one | 15 dimension :one |
| 16 data length :power of 2 | 16 data length :power of 2 |
| 17 decimation :frequency | 17 decimation :frequency |
| 18 radix :4, 2 | 18 radix :4, 2 |
| 19 data :inplace | 19 data :inplace |
| 20 table :use | 20 table :use |
| 21 functions | 21 functions |
| 22 cdft: Complex Discrete Fourier Transform | 22 cdft: Complex Discrete Fourier Transform |
| 23 rdft: Real Discrete Fourier Transform | 23 rdft: Real Discrete Fourier Transform |
| 24 ddct: Discrete Cosine Transform | 24 ddct: Discrete Cosine Transform |
| 25 ddst: Discrete Sine Transform | 25 ddst: Discrete Sine Transform |
| 26 dfct: Cosine Transform of RDFT (Real Symmetric DFT) | 26 dfct: Cosine Transform of RDFT (Real Symmetric DFT) |
| 27 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) | 27 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) |
| 28 function prototypes | 28 function prototypes |
| 29 void cdft(int, int, float *, int *, float *); | 29 void cdft(int, int, float *, int *, float *); |
| 30 void rdft(int, int, float *, int *, float *); | 30 void rdft(size_t, int, float *, size_t *, float *); |
| 31 void ddct(int, int, float *, int *, float *); | 31 void ddct(int, int, float *, int *, float *); |
| 32 void ddst(int, int, float *, int *, float *); | 32 void ddst(int, int, float *, int *, float *); |
| 33 void dfct(int, float *, float *, int *, float *); | 33 void dfct(int, float *, float *, int *, float *); |
| 34 void dfst(int, float *, float *, int *, float *); | 34 void dfst(int, float *, float *, int *, float *); |
| 35 | 35 |
| 36 | 36 |
| 37 -------- Complex DFT (Discrete Fourier Transform) -------- | 37 -------- Complex DFT (Discrete Fourier Transform) -------- |
| 38 [definition] | 38 [definition] |
| 39 <case1> | 39 <case1> |
| 40 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n | 40 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n |
| (...skipping 46 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 87 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + | 87 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + |
| 88 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n | 88 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n |
| 89 [usage] | 89 [usage] |
| 90 <case1> | 90 <case1> |
| 91 ip[0] = 0; // first time only | 91 ip[0] = 0; // first time only |
| 92 rdft(n, 1, a, ip, w); | 92 rdft(n, 1, a, ip, w); |
| 93 <case2> | 93 <case2> |
| 94 ip[0] = 0; // first time only | 94 ip[0] = 0; // first time only |
| 95 rdft(n, -1, a, ip, w); | 95 rdft(n, -1, a, ip, w); |
| 96 [parameters] | 96 [parameters] |
| 97 n :data length (int) | 97 n :data length (size_t) |
| 98 n >= 2, n = power of 2 | 98 n >= 2, n = power of 2 |
| 99 a[0...n-1] :input/output data (float *) | 99 a[0...n-1] :input/output data (float *) |
| 100 <case1> | 100 <case1> |
| 101 output data | 101 output data |
| 102 a[2*k] = R[k], 0<=k<n/2 | 102 a[2*k] = R[k], 0<=k<n/2 |
| 103 a[2*k+1] = I[k], 0<k<n/2 | 103 a[2*k+1] = I[k], 0<k<n/2 |
| 104 a[1] = R[n/2] | 104 a[1] = R[n/2] |
| 105 <case2> | 105 <case2> |
| 106 input data | 106 input data |
| 107 a[2*j] = R[j], 0<=j<n/2 | 107 a[2*j] = R[j], 0<=j<n/2 |
| 108 a[2*j+1] = I[j], 0<j<n/2 | 108 a[2*j+1] = I[j], 0<j<n/2 |
| 109 a[1] = R[n/2] | 109 a[1] = R[n/2] |
| 110 ip[0...*] :work area for bit reversal (int *) | 110 ip[0...*] :work area for bit reversal (size_t *) |
| 111 length of ip >= 2+sqrt(n/2) | 111 length of ip >= 2+sqrt(n/2) |
| 112 strictly, | 112 strictly, |
| 113 length of ip >= | 113 length of ip >= |
| 114 2+(1<<(int)(log(n/2+0.5)/log(2))/2). | 114 2+(1<<(int)(log(n/2+0.5)/log(2))/2). |
| 115 ip[0],ip[1] are pointers of the cos/sin table. | 115 ip[0],ip[1] are pointers of the cos/sin table. |
| 116 w[0...n/2-1] :cos/sin table (float *) | 116 w[0...n/2-1] :cos/sin table (float *) |
| 117 w[],ip[] are initialized if ip[0] == 0. | 117 w[],ip[] are initialized if ip[0] == 0. |
| 118 [remark] | 118 [remark] |
| 119 Inverse of | 119 Inverse of |
| 120 rdft(n, 1, a, ip, w); | 120 rdft(n, 1, a, ip, w); |
| (...skipping 158 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); |
| 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 |