| 1425 | |
| 1426 | |
| 1427 | void cftmdl2(int n, float *a) |
| 1428 | { |
| 1429 | int i, i0, j, j0, j1, j2, j3, m, mh; |
| 1430 | float ew, w1r, w1i, wn4r, wk1r, wk1i, wk3r, wk3i, |
| 1431 | wl1r, wl1i, wl3r, wl3i, wd1r, wd1i, wd3r, wd3i, |
| 1432 | we1r, we1i, we3r, we3i, ss1, ss3; |
| 1433 | float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i; |
| 1434 | |
| 1435 | mh = n >> 3; |
| 1436 | m = 2 * mh; |
| 1437 | wn4r = WR5000; |
| 1438 | j1 = m; |
| 1439 | j2 = j1 + m; |
| 1440 | j3 = j2 + m; |
| 1441 | x0r = a[0] - a[j2 + 1]; |
| 1442 | x0i = a[1] + a[j2]; |
| 1443 | x1r = a[0] + a[j2 + 1]; |
| 1444 | x1i = a[1] - a[j2]; |
| 1445 | x2r = a[j1] - a[j3 + 1]; |
| 1446 | x2i = a[j1 + 1] + a[j3]; |
| 1447 | x3r = a[j1] + a[j3 + 1]; |
| 1448 | x3i = a[j1 + 1] - a[j3]; |
| 1449 | y0r = wn4r * (x2r - x2i); |
| 1450 | y0i = wn4r * (x2i + x2r); |
| 1451 | a[0] = x0r + y0r; |
| 1452 | a[1] = x0i + y0i; |
| 1453 | a[j1] = x0r - y0r; |
| 1454 | a[j1 + 1] = x0i - y0i; |
| 1455 | y0r = wn4r * (x3r - x3i); |
| 1456 | y0i = wn4r * (x3i + x3r); |
| 1457 | a[j2] = x1r - y0i; |
| 1458 | a[j2 + 1] = x1i + y0r; |
| 1459 | a[j3] = x1r + y0i; |
| 1460 | a[j3 + 1] = x1i - y0r; |
| 1461 | wl1r = 1; |
| 1462 | wl1i = 0; |
| 1463 | wl3r = 1; |
| 1464 | wl3i = 0; |
| 1465 | we1r = wn4r; |
| 1466 | we1i = wn4r; |
| 1467 | we3r = -wn4r; |
| 1468 | we3i = -wn4r; |
| 1469 | ew = M_PI_2 / (2 * m); |
| 1470 | w1r = (float)cos(2 * ew); |
| 1471 | w1i = (float)sin(2 * ew); |
| 1472 | wk1r = w1r; |
| 1473 | wk1i = w1i; |
| 1474 | wd1r = wn4r * (w1r - w1i); |
| 1475 | wd1i = wn4r * (w1i + w1r); |
| 1476 | ss1 = 2 * w1i; |
| 1477 | wk3i = 2 * ss1 * wk1r; |
| 1478 | wk3r = wk1r - wk3i * wk1i; |
| 1479 | wk3i = wk1i - wk3i * wk1r; |
| 1480 | ss3 = 2 * wk3i; |
| 1481 | wd3r = -wn4r * (wk3r - wk3i); |
| 1482 | wd3i = -wn4r * (wk3i + wk3r); |
| 1483 | i = 0; |
| 1484 | for (;;) { |