| 970 | } |
| 971 | |
| 972 | void cftb1st(int n, float *a) |
| 973 | { |
| 974 | int i, i0, j, j0, j1, j2, j3, m, mh; |
| 975 | float ew, w1r, w1i, wk1r, wk1i, wk3r, wk3i, |
| 976 | wd1r, wd1i, wd3r, wd3i, ss1, ss3; |
| 977 | float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
| 978 | |
| 979 | mh = n >> 3; |
| 980 | m = 2 * mh; |
| 981 | j1 = m; |
| 982 | j2 = j1 + m; |
| 983 | j3 = j2 + m; |
| 984 | x0r = a[0] + a[j2]; |
| 985 | x0i = -a[1] - a[j2 + 1]; |
| 986 | x1r = a[0] - a[j2]; |
| 987 | x1i = -a[1] + a[j2 + 1]; |
| 988 | x2r = a[j1] + a[j3]; |
| 989 | x2i = a[j1 + 1] + a[j3 + 1]; |
| 990 | x3r = a[j1] - a[j3]; |
| 991 | x3i = a[j1 + 1] - a[j3 + 1]; |
| 992 | a[0] = x0r + x2r; |
| 993 | a[1] = x0i - x2i; |
| 994 | a[j1] = x0r - x2r; |
| 995 | a[j1 + 1] = x0i + x2i; |
| 996 | a[j2] = x1r + x3i; |
| 997 | a[j2 + 1] = x1i + x3r; |
| 998 | a[j3] = x1r - x3i; |
| 999 | a[j3 + 1] = x1i - x3r; |
| 1000 | wd1r = 1; |
| 1001 | wd1i = 0; |
| 1002 | wd3r = 1; |
| 1003 | wd3i = 0; |
| 1004 | ew = M_PI_2 / m; |
| 1005 | w1r = (float)cos(2 * ew); |
| 1006 | w1i = (float)sin(2 * ew); |
| 1007 | wk1r = w1r; |
| 1008 | wk1i = w1i; |
| 1009 | ss1 = 2 * w1i; |
| 1010 | wk3i = 2 * ss1 * wk1r; |
| 1011 | wk3r = wk1r - wk3i * wk1i; |
| 1012 | wk3i = wk1i - wk3i * wk1r; |
| 1013 | ss3 = 2 * wk3i; |
| 1014 | i = 0; |
| 1015 | for (;;) { |
| 1016 | i0 = i + 4 * CDFT_LOOP_DIV; |
| 1017 | if (i0 > mh - 4) { |
| 1018 | i0 = mh - 4; |
| 1019 | } |
| 1020 | for (j = i + 2; j < i0; j += 4) { |
| 1021 | wd1r -= ss1 * wk1i; |
| 1022 | wd1i += ss1 * wk1r; |
| 1023 | wd3r -= ss3 * wk3i; |
| 1024 | wd3i += ss3 * wk3r; |
| 1025 | j1 = j + m; |
| 1026 | j2 = j1 + m; |
| 1027 | j3 = j2 + m; |
| 1028 | x0r = a[j] + a[j2]; |
| 1029 | x0i = -a[j + 1] - a[j2 + 1]; |