285 extern "C" void cdft(
int n,
int isgn,
double *a,
int *ip,
double *w)
287 void makewt(
int nw,
int *ip,
double *w);
288 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
289 void cftbsub(
int n,
double *a,
int *ip,
int nw,
double *w);
298 cftfsub(n, a, ip, nw, w);
300 cftbsub(n, a, ip, nw, w);
305 extern "C" void rdft(
int n,
int isgn,
double *a,
int *ip,
double *w)
307 void makewt(
int nw,
int *ip,
double *w);
308 void makect(
int nc,
int *ip,
double *c);
309 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
310 void cftbsub(
int n,
double *a,
int *ip,
int nw,
double *w);
311 void rftfsub(
int n,
double *a,
int nc,
double *c);
312 void rftbsub(
int n,
double *a,
int nc,
double *c);
324 makect(nc, ip, w + nw);
328 cftfsub(n, a, ip, nw, w);
329 rftfsub(n, a, nc, w + nw);
331 cftfsub(n, a, ip, nw, w);
337 a[1] = 0.5 * (a[0] - a[1]);
340 rftbsub(n, a, nc, w + nw);
341 cftbsub(n, a, ip, nw, w);
343 cftbsub(n, a, ip, nw, w);
349 extern "C" void ddct(
int n,
int isgn,
double *a,
int *ip,
double *w)
351 void makewt(
int nw,
int *ip,
double *w);
352 void makect(
int nc,
int *ip,
double *c);
353 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
354 void cftbsub(
int n,
double *a,
int *ip,
int nw,
double *w);
355 void rftfsub(
int n,
double *a,
int nc,
double *c);
356 void rftbsub(
int n,
double *a,
int nc,
double *c);
357 void dctsub(
int n,
double *a,
int nc,
double *c);
369 makect(nc, ip, w + nw);
373 for (j = n - 2; j >= 2; j -= 2) {
374 a[j + 1] = a[j] - a[j - 1];
380 rftbsub(n, a, nc, w + nw);
381 cftbsub(n, a, ip, nw, w);
383 cftbsub(n, a, ip, nw, w);
386 dctsub(n, a, nc, w + nw);
389 cftfsub(n, a, ip, nw, w);
390 rftfsub(n, a, nc, w + nw);
392 cftfsub(n, a, ip, nw, w);
396 for (j = 2; j < n; j += 2) {
397 a[j - 1] = a[j] - a[j + 1];
405 extern "C" void ddst(
int n,
int isgn,
double *a,
int *ip,
double *w)
407 void makewt(
int nw,
int *ip,
double *w);
408 void makect(
int nc,
int *ip,
double *c);
409 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
410 void cftbsub(
int n,
double *a,
int *ip,
int nw,
double *w);
411 void rftfsub(
int n,
double *a,
int nc,
double *c);
412 void rftbsub(
int n,
double *a,
int nc,
double *c);
413 void dstsub(
int n,
double *a,
int nc,
double *c);
425 makect(nc, ip, w + nw);
429 for (j = n - 2; j >= 2; j -= 2) {
430 a[j + 1] = -a[j] - a[j - 1];
436 rftbsub(n, a, nc, w + nw);
437 cftbsub(n, a, ip, nw, w);
439 cftbsub(n, a, ip, nw, w);
442 dstsub(n, a, nc, w + nw);
445 cftfsub(n, a, ip, nw, w);
446 rftfsub(n, a, nc, w + nw);
448 cftfsub(n, a, ip, nw, w);
452 for (j = 2; j < n; j += 2) {
453 a[j - 1] = -a[j] - a[j + 1];
461 extern "C" void dfct(
int n,
double *a,
double *t,
int *ip,
double *w)
463 void makewt(
int nw,
int *ip,
double *w);
464 void makect(
int nc,
int *ip,
double *c);
465 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
466 void rftfsub(
int n,
double *a,
int nc,
double *c);
467 void dctsub(
int n,
double *a,
int nc,
double *c);
468 int j, k, l, m, mh, nw, nc;
469 double xr, xi, yr, yi;
479 makect(nc, ip, w + nw);
489 for (j = 1; j < mh; j++) {
491 xr = a[j] - a[n - j];
492 xi = a[j] + a[n - j];
493 yr = a[k] - a[n - k];
494 yi = a[k] + a[n - k];
500 t[mh] = a[mh] + a[n - mh];
502 dctsub(m, a, nc, w + nw);
504 cftfsub(m, a, ip, nw, w);
505 rftfsub(m, a, nc, w + nw);
507 cftfsub(m, a, ip, nw, w);
509 a[n - 1] = a[0] - a[1];
511 for (j = m - 2; j >= 2; j -= 2) {
512 a[2 * j + 1] = a[j] + a[j + 1];
513 a[2 * j - 1] = a[j] - a[j + 1];
518 dctsub(m, t, nc, w + nw);
520 cftfsub(m, t, ip, nw, w);
521 rftfsub(m, t, nc, w + nw);
523 cftfsub(m, t, ip, nw, w);
525 a[n - l] = t[0] - t[1];
528 for (j = 2; j < m; j += 2) {
530 a[k - l] = t[j] - t[j + 1];
531 a[k + l] = t[j] + t[j + 1];
535 for (j = 0; j < mh; j++) {
537 t[j] = t[m + k] - t[m + j];
538 t[k] = t[m + k] + t[m + j];
554 extern "C" void dfst(
int n,
double *a,
double *t,
int *ip,
double *w)
556 void makewt(
int nw,
int *ip,
double *w);
557 void makect(
int nc,
int *ip,
double *c);
558 void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w);
559 void rftfsub(
int n,
double *a,
int nc,
double *c);
560 void dstsub(
int n,
double *a,
int nc,
double *c);
561 int j, k, l, m, mh, nw, nc;
562 double xr, xi, yr, yi;
572 makect(nc, ip, w + nw);
577 for (j = 1; j < mh; j++) {
579 xr = a[j] + a[n - j];
580 xi = a[j] - a[n - j];
581 yr = a[k] + a[n - k];
582 yi = a[k] - a[n - k];
588 t[0] = a[mh] - a[n - mh];
591 dstsub(m, a, nc, w + nw);
593 cftfsub(m, a, ip, nw, w);
594 rftfsub(m, a, nc, w + nw);
596 cftfsub(m, a, ip, nw, w);
598 a[n - 1] = a[1] - a[0];
600 for (j = m - 2; j >= 2; j -= 2) {
601 a[2 * j + 1] = a[j] - a[j + 1];
602 a[2 * j - 1] = -a[j] - a[j + 1];
607 dstsub(m, t, nc, w + nw);
609 cftfsub(m, t, ip, nw, w);
610 rftfsub(m, t, nc, w + nw);
612 cftfsub(m, t, ip, nw, w);
614 a[n - l] = t[1] - t[0];
617 for (j = 2; j < m; j += 2) {
619 a[k - l] = -t[j] - t[j + 1];
620 a[k + l] = t[j] - t[j + 1];
624 for (j = 1; j < mh; j++) {
626 t[j] = t[m + k] + t[m + j];
627 t[k] = t[m + k] - t[m + j];
643 extern "C" void makewt(
int nw,
int *ip,
double *w)
645 void makeipt(
int nw,
int *ip);
646 int j, nwh, nw0, nw1;
647 double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
653 delta = atan(1.0) / nwh;
654 wn4r = cos(delta * nwh);
658 w[2] = cos(delta * 2);
659 w[3] = sin(delta * 2);
660 }
else if (nwh > 4) {
662 w[2] = 0.5 / cos(delta * 2);
663 w[3] = 0.5 / cos(delta * 6);
664 for (j = 4; j < nwh; j += 4) {
665 w[j] = cos(delta * j);
666 w[j + 1] = sin(delta * j);
667 w[j + 2] = cos(3 * delta * j);
668 w[j + 3] = -sin(3 * delta * j);
682 }
else if (nwh > 4) {
685 w[nw1 + 2] = 0.5 / wk1r;
686 w[nw1 + 3] = 0.5 / wk3r;
687 for (j = 4; j < nwh; j += 4) {
688 wk1r = w[nw0 + 2 * j];
689 wk1i = w[nw0 + 2 * j + 1];
690 wk3r = w[nw0 + 2 * j + 2];
691 wk3i = w[nw0 + 2 * j + 3];
693 w[nw1 + j + 1] = wk1i;
694 w[nw1 + j + 2] = wk3r;
695 w[nw1 + j + 3] = wk3i;
704 extern "C" void makeipt(
int nw,
int *ip)
706 int j, l, m, m2, p, q;
711 for (l = nw; l > 32; l >>= 2) {
714 for (j = m; j < m2; j++) {
724 extern "C" void makect(
int nc,
int *ip,
double *c)
732 delta = atan(1.0) / nch;
733 c[0] = cos(delta * nch);
735 for (j = 1; j < nch; j++) {
736 c[j] = 0.5 * cos(delta * j);
737 c[nc - j] = 0.5 * sin(delta * j);
746 #ifdef USE_CDFT_PTHREADS
747 #define USE_CDFT_THREADS
748 #ifndef CDFT_THREADS_BEGIN_N
749 #define CDFT_THREADS_BEGIN_N 8192
751 #ifndef CDFT_4THREADS_BEGIN_N
752 #define CDFT_4THREADS_BEGIN_N 65536
757 #define cdft_thread_t pthread_t
758 #define cdft_thread_create(thp,func,argp) { \
759 if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
760 fprintf(stderr, "cdft thread error\n"); \
764 #define cdft_thread_wait(th) { \
765 if (pthread_join(th, NULL) != 0) { \
766 fprintf(stderr, "cdft thread error\n"); \
773 #ifdef USE_CDFT_WINTHREADS
774 #define USE_CDFT_THREADS
775 #ifndef CDFT_THREADS_BEGIN_N
776 #define CDFT_THREADS_BEGIN_N 32768
778 #ifndef CDFT_4THREADS_BEGIN_N
779 #define CDFT_4THREADS_BEGIN_N 524288
784 #define cdft_thread_t HANDLE
785 #define cdft_thread_create(thp,func,argp) { \
787 *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
789 fprintf(stderr, "cdft thread error\n"); \
793 #define cdft_thread_wait(th) { \
794 WaitForSingleObject(th, INFINITE); \
800 extern "C" void cftfsub(
int n,
double *a,
int *ip,
int nw,
double *w)
802 void bitrv2(
int n,
int *ip,
double *a);
803 void bitrv216(
double *a);
804 void bitrv208(
double *a);
805 void cftf1st(
int n,
double *a,
double *w);
806 void cftrec4(
int n,
double *a,
int nw,
double *w);
807 void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w);
808 void cftfx41(
int n,
double *a,
int nw,
double *w);
809 void cftf161(
double *a,
double *w);
810 void cftf081(
double *a,
double *w);
811 void cftf040(
double *a);
812 void cftx020(
double *a);
813 #ifdef USE_CDFT_THREADS
814 void cftrec4_th(
int n,
double *a,
int nw,
double *w);
819 cftf1st(n, a, &w[nw - (n >> 2)]);
820 #ifdef USE_CDFT_THREADS
821 if (n > CDFT_THREADS_BEGIN_N) {
822 cftrec4_th(n, a, nw, w);
826 cftrec4(n, a, nw, w);
827 }
else if (n > 128) {
828 cftleaf(n, 1, a, nw, w);
830 cftfx41(n, a, nw, w);
833 }
else if (n == 32) {
834 cftf161(a, &w[nw - 8]);
848 extern "C" void cftbsub(
int n,
double *a,
int *ip,
int nw,
double *w)
850 void bitrv2conj(
int n,
int *ip,
double *a);
851 void bitrv216neg(
double *a);
852 void bitrv208neg(
double *a);
853 void cftb1st(
int n,
double *a,
double *w);
854 void cftrec4(
int n,
double *a,
int nw,
double *w);
855 void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w);
856 void cftfx41(
int n,
double *a,
int nw,
double *w);
857 void cftf161(
double *a,
double *w);
858 void cftf081(
double *a,
double *w);
859 void cftb040(
double *a);
860 void cftx020(
double *a);
861 #ifdef USE_CDFT_THREADS
862 void cftrec4_th(
int n,
double *a,
int nw,
double *w);
867 cftb1st(n, a, &w[nw - (n >> 2)]);
868 #ifdef USE_CDFT_THREADS
869 if (n > CDFT_THREADS_BEGIN_N) {
870 cftrec4_th(n, a, nw, w);
874 cftrec4(n, a, nw, w);
875 }
else if (n > 128) {
876 cftleaf(n, 1, a, nw, w);
878 cftfx41(n, a, nw, w);
880 bitrv2conj(n, ip, a);
881 }
else if (n == 32) {
882 cftf161(a, &w[nw - 8]);
896 extern "C" void bitrv2(
int n,
int *ip,
double *a)
898 int j, j1, k, k1, l, m, nh, nm;
899 double xr, xi, yr, yi;
902 for (l = n >> 2; l > 8; l >>= 2) {
908 for (k = 0; k < m; k++) {
909 for (j = 0; j < k; j++) {
910 j1 = 4 * j + 2 * ip[m + k];
911 k1 = 4 * k + 2 * ip[m + j];
1071 k1 = 4 * k + 2 * ip[m + k];
1134 for (k = 0; k < m; k++) {
1135 for (j = 0; j < k; j++) {
1136 j1 = 4 * j + ip[m + k];
1137 k1 = 4 * k + ip[m + j];
1217 k1 = 4 * k + ip[m + k];
1243 extern "C" void bitrv2conj(
int n,
int *ip,
double *a)
1245 int j, j1, k, k1, l, m, nh, nm;
1246 double xr, xi, yr, yi;
1249 for (l = n >> 2; l > 8; l >>= 2) {
1255 for (k = 0; k < m; k++) {
1256 for (j = 0; j < k; j++) {
1257 j1 = 4 * j + 2 * ip[m + k];
1258 k1 = 4 * k + 2 * ip[m + j];
1418 k1 = 4 * k + 2 * ip[m + k];
1421 a[j1 - 1] = -a[j1 - 1];
1430 a[k1 + 3] = -a[k1 + 3];
1473 a[j1 - 1] = -a[j1 - 1];
1482 a[k1 + 3] = -a[k1 + 3];
1485 for (k = 0; k < m; k++) {
1486 for (j = 0; j < k; j++) {
1487 j1 = 4 * j + ip[m + k];
1488 k1 = 4 * k + ip[m + j];
1568 k1 = 4 * k + ip[m + k];
1571 a[j1 - 1] = -a[j1 - 1];
1580 a[k1 + 3] = -a[k1 + 3];
1583 a[j1 - 1] = -a[j1 - 1];
1592 a[k1 + 3] = -a[k1 + 3];
1598 extern "C" void bitrv216(
double *a)
1600 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1601 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
1602 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1655 extern "C" void bitrv216neg(
double *a)
1657 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1658 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
1659 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
1660 x13r, x13i, x14r, x14i, x15r, x15i;
1725 extern "C" void bitrv208(
double *a)
1727 double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1748 extern "C" void bitrv208neg(
double *a)
1750 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1751 x5r, x5i, x6r, x6i, x7r, x7i;
1784 extern "C" void cftf1st(
int n,
double *a,
double *w)
1786 int j, j0, j1, j2, j3, k, m, mh;
1787 double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
1788 wd1r, wd1i, wd3r, wd3i;
1789 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1790 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1798 x0i = a[1] + a[j2 + 1];
1800 x1i = a[1] - a[j2 + 1];
1801 x2r = a[j1] + a[j3];
1802 x2i = a[j1 + 1] + a[j3 + 1];
1803 x3r = a[j1] - a[j3];
1804 x3i = a[j1 + 1] - a[j3 + 1];
1808 a[j1 + 1] = x0i - x2i;
1810 a[j2 + 1] = x1i + x3r;
1812 a[j3 + 1] = x1i - x3r;
1821 for (j = 2; j < mh - 2; j += 4) {
1823 wk1r = csc1 * (wd1r + w[k]);
1824 wk1i = csc1 * (wd1i + w[k + 1]);
1825 wk3r = csc3 * (wd3r + w[k + 2]);
1826 wk3i = csc3 * (wd3i + w[k + 3]);
1835 x0i = a[j + 1] + a[j2 + 1];
1837 x1i = a[j + 1] - a[j2 + 1];
1838 y0r = a[j + 2] + a[j2 + 2];
1839 y0i = a[j + 3] + a[j2 + 3];
1840 y1r = a[j + 2] - a[j2 + 2];
1841 y1i = a[j + 3] - a[j2 + 3];
1842 x2r = a[j1] + a[j3];
1843 x2i = a[j1 + 1] + a[j3 + 1];
1844 x3r = a[j1] - a[j3];
1845 x3i = a[j1 + 1] - a[j3 + 1];
1846 y2r = a[j1 + 2] + a[j3 + 2];
1847 y2i = a[j1 + 3] + a[j3 + 3];
1848 y3r = a[j1 + 2] - a[j3 + 2];
1849 y3i = a[j1 + 3] - a[j3 + 3];
1851 a[j + 1] = x0i + x2i;
1852 a[j + 2] = y0r + y2r;
1853 a[j + 3] = y0i + y2i;
1855 a[j1 + 1] = x0i - x2i;
1856 a[j1 + 2] = y0r - y2r;
1857 a[j1 + 3] = y0i - y2i;
1860 a[j2] = wk1r * x0r - wk1i * x0i;
1861 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
1864 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
1865 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
1868 a[j3] = wk3r * x0r + wk3i * x0i;
1869 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
1872 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
1873 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
1878 x0r = a[j0] + a[j2];
1879 x0i = a[j0 + 1] + a[j2 + 1];
1880 x1r = a[j0] - a[j2];
1881 x1i = a[j0 + 1] - a[j2 + 1];
1882 y0r = a[j0 - 2] + a[j2 - 2];
1883 y0i = a[j0 - 1] + a[j2 - 1];
1884 y1r = a[j0 - 2] - a[j2 - 2];
1885 y1i = a[j0 - 1] - a[j2 - 1];
1886 x2r = a[j1] + a[j3];
1887 x2i = a[j1 + 1] + a[j3 + 1];
1888 x3r = a[j1] - a[j3];
1889 x3i = a[j1 + 1] - a[j3 + 1];
1890 y2r = a[j1 - 2] + a[j3 - 2];
1891 y2i = a[j1 - 1] + a[j3 - 1];
1892 y3r = a[j1 - 2] - a[j3 - 2];
1893 y3i = a[j1 - 1] - a[j3 - 1];
1895 a[j0 + 1] = x0i + x2i;
1896 a[j0 - 2] = y0r + y2r;
1897 a[j0 - 1] = y0i + y2i;
1899 a[j1 + 1] = x0i - x2i;
1900 a[j1 - 2] = y0r - y2r;
1901 a[j1 - 1] = y0i - y2i;
1904 a[j2] = wk1i * x0r - wk1r * x0i;
1905 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
1908 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
1909 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
1912 a[j3] = wk3i * x0r + wk3r * x0i;
1913 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
1916 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
1917 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
1919 wk1r = csc1 * (wd1r + wn4r);
1920 wk1i = csc1 * (wd1i + wn4r);
1921 wk3r = csc3 * (wd3r - wn4r);
1922 wk3i = csc3 * (wd3i - wn4r);
1927 x0r = a[j0 - 2] + a[j2 - 2];
1928 x0i = a[j0 - 1] + a[j2 - 1];
1929 x1r = a[j0 - 2] - a[j2 - 2];
1930 x1i = a[j0 - 1] - a[j2 - 1];
1931 x2r = a[j1 - 2] + a[j3 - 2];
1932 x2i = a[j1 - 1] + a[j3 - 1];
1933 x3r = a[j1 - 2] - a[j3 - 2];
1934 x3i = a[j1 - 1] - a[j3 - 1];
1935 a[j0 - 2] = x0r + x2r;
1936 a[j0 - 1] = x0i + x2i;
1937 a[j1 - 2] = x0r - x2r;
1938 a[j1 - 1] = x0i - x2i;
1941 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
1942 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
1945 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
1946 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
1947 x0r = a[j0] + a[j2];
1948 x0i = a[j0 + 1] + a[j2 + 1];
1949 x1r = a[j0] - a[j2];
1950 x1i = a[j0 + 1] - a[j2 + 1];
1951 x2r = a[j1] + a[j3];
1952 x2i = a[j1 + 1] + a[j3 + 1];
1953 x3r = a[j1] - a[j3];
1954 x3i = a[j1 + 1] - a[j3 + 1];
1956 a[j0 + 1] = x0i + x2i;
1958 a[j1 + 1] = x0i - x2i;
1961 a[j2] = wn4r * (x0r - x0i);
1962 a[j2 + 1] = wn4r * (x0i + x0r);
1965 a[j3] = -wn4r * (x0r + x0i);
1966 a[j3 + 1] = -wn4r * (x0i - x0r);
1967 x0r = a[j0 + 2] + a[j2 + 2];
1968 x0i = a[j0 + 3] + a[j2 + 3];
1969 x1r = a[j0 + 2] - a[j2 + 2];
1970 x1i = a[j0 + 3] - a[j2 + 3];
1971 x2r = a[j1 + 2] + a[j3 + 2];
1972 x2i = a[j1 + 3] + a[j3 + 3];
1973 x3r = a[j1 + 2] - a[j3 + 2];
1974 x3i = a[j1 + 3] - a[j3 + 3];
1975 a[j0 + 2] = x0r + x2r;
1976 a[j0 + 3] = x0i + x2i;
1977 a[j1 + 2] = x0r - x2r;
1978 a[j1 + 3] = x0i - x2i;
1981 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
1982 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
1985 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
1986 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
1990 extern "C" void cftb1st(
int n,
double *a,
double *w)
1992 int j, j0, j1, j2, j3, k, m, mh;
1993 double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
1994 wd1r, wd1i, wd3r, wd3i;
1995 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1996 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
2004 x0i = -a[1] - a[j2 + 1];
2006 x1i = -a[1] + a[j2 + 1];
2007 x2r = a[j1] + a[j3];
2008 x2i = a[j1 + 1] + a[j3 + 1];
2009 x3r = a[j1] - a[j3];
2010 x3i = a[j1 + 1] - a[j3 + 1];
2014 a[j1 + 1] = x0i + x2i;
2016 a[j2 + 1] = x1i + x3r;
2018 a[j3 + 1] = x1i - x3r;
2027 for (j = 2; j < mh - 2; j += 4) {
2029 wk1r = csc1 * (wd1r + w[k]);
2030 wk1i = csc1 * (wd1i + w[k + 1]);
2031 wk3r = csc3 * (wd3r + w[k + 2]);
2032 wk3i = csc3 * (wd3i + w[k + 3]);
2041 x0i = -a[j + 1] - a[j2 + 1];
2043 x1i = -a[j + 1] + a[j2 + 1];
2044 y0r = a[j + 2] + a[j2 + 2];
2045 y0i = -a[j + 3] - a[j2 + 3];
2046 y1r = a[j + 2] - a[j2 + 2];
2047 y1i = -a[j + 3] + a[j2 + 3];
2048 x2r = a[j1] + a[j3];
2049 x2i = a[j1 + 1] + a[j3 + 1];
2050 x3r = a[j1] - a[j3];
2051 x3i = a[j1 + 1] - a[j3 + 1];
2052 y2r = a[j1 + 2] + a[j3 + 2];
2053 y2i = a[j1 + 3] + a[j3 + 3];
2054 y3r = a[j1 + 2] - a[j3 + 2];
2055 y3i = a[j1 + 3] - a[j3 + 3];
2057 a[j + 1] = x0i - x2i;
2058 a[j + 2] = y0r + y2r;
2059 a[j + 3] = y0i - y2i;
2061 a[j1 + 1] = x0i + x2i;
2062 a[j1 + 2] = y0r - y2r;
2063 a[j1 + 3] = y0i + y2i;
2066 a[j2] = wk1r * x0r - wk1i * x0i;
2067 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2070 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2071 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2074 a[j3] = wk3r * x0r + wk3i * x0i;
2075 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2078 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2079 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2084 x0r = a[j0] + a[j2];
2085 x0i = -a[j0 + 1] - a[j2 + 1];
2086 x1r = a[j0] - a[j2];
2087 x1i = -a[j0 + 1] + a[j2 + 1];
2088 y0r = a[j0 - 2] + a[j2 - 2];
2089 y0i = -a[j0 - 1] - a[j2 - 1];
2090 y1r = a[j0 - 2] - a[j2 - 2];
2091 y1i = -a[j0 - 1] + a[j2 - 1];
2092 x2r = a[j1] + a[j3];
2093 x2i = a[j1 + 1] + a[j3 + 1];
2094 x3r = a[j1] - a[j3];
2095 x3i = a[j1 + 1] - a[j3 + 1];
2096 y2r = a[j1 - 2] + a[j3 - 2];
2097 y2i = a[j1 - 1] + a[j3 - 1];
2098 y3r = a[j1 - 2] - a[j3 - 2];
2099 y3i = a[j1 - 1] - a[j3 - 1];
2101 a[j0 + 1] = x0i - x2i;
2102 a[j0 - 2] = y0r + y2r;
2103 a[j0 - 1] = y0i - y2i;
2105 a[j1 + 1] = x0i + x2i;
2106 a[j1 - 2] = y0r - y2r;
2107 a[j1 - 1] = y0i + y2i;
2110 a[j2] = wk1i * x0r - wk1r * x0i;
2111 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2114 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2115 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2118 a[j3] = wk3i * x0r + wk3r * x0i;
2119 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2122 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2123 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2125 wk1r = csc1 * (wd1r + wn4r);
2126 wk1i = csc1 * (wd1i + wn4r);
2127 wk3r = csc3 * (wd3r - wn4r);
2128 wk3i = csc3 * (wd3i - wn4r);
2133 x0r = a[j0 - 2] + a[j2 - 2];
2134 x0i = -a[j0 - 1] - a[j2 - 1];
2135 x1r = a[j0 - 2] - a[j2 - 2];
2136 x1i = -a[j0 - 1] + a[j2 - 1];
2137 x2r = a[j1 - 2] + a[j3 - 2];
2138 x2i = a[j1 - 1] + a[j3 - 1];
2139 x3r = a[j1 - 2] - a[j3 - 2];
2140 x3i = a[j1 - 1] - a[j3 - 1];
2141 a[j0 - 2] = x0r + x2r;
2142 a[j0 - 1] = x0i - x2i;
2143 a[j1 - 2] = x0r - x2r;
2144 a[j1 - 1] = x0i + x2i;
2147 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2148 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2151 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2152 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2153 x0r = a[j0] + a[j2];
2154 x0i = -a[j0 + 1] - a[j2 + 1];
2155 x1r = a[j0] - a[j2];
2156 x1i = -a[j0 + 1] + a[j2 + 1];
2157 x2r = a[j1] + a[j3];
2158 x2i = a[j1 + 1] + a[j3 + 1];
2159 x3r = a[j1] - a[j3];
2160 x3i = a[j1 + 1] - a[j3 + 1];
2162 a[j0 + 1] = x0i - x2i;
2164 a[j1 + 1] = x0i + x2i;
2167 a[j2] = wn4r * (x0r - x0i);
2168 a[j2 + 1] = wn4r * (x0i + x0r);
2171 a[j3] = -wn4r * (x0r + x0i);
2172 a[j3 + 1] = -wn4r * (x0i - x0r);
2173 x0r = a[j0 + 2] + a[j2 + 2];
2174 x0i = -a[j0 + 3] - a[j2 + 3];
2175 x1r = a[j0 + 2] - a[j2 + 2];
2176 x1i = -a[j0 + 3] + a[j2 + 3];
2177 x2r = a[j1 + 2] + a[j3 + 2];
2178 x2i = a[j1 + 3] + a[j3 + 3];
2179 x3r = a[j1 + 2] - a[j3 + 2];
2180 x3i = a[j1 + 3] - a[j3 + 3];
2181 a[j0 + 2] = x0r + x2r;
2182 a[j0 + 3] = x0i - x2i;
2183 a[j1 + 2] = x0r - x2r;
2184 a[j1 + 3] = x0i + x2i;
2187 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2188 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2191 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2192 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2196 #ifdef USE_CDFT_THREADS
2197 struct cdft_arg_st {
2204 typedef struct cdft_arg_st cdft_arg_t;
2207 extern "C" void cftrec4_th(
int n,
double *a,
int nw,
double *w)
2209 void *cftrec1_th(
void *p);
2210 void *cftrec2_th(
void *p);
2211 int i, idiv4, m, nthread;
2212 cdft_thread_t th[4];
2218 if (n > CDFT_4THREADS_BEGIN_N) {
2223 for (i = 0; i < nthread; i++) {
2226 ag[i].a = &a[i * m];
2230 cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
2232 cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
2235 for (i = 0; i < nthread; i++) {
2236 cdft_thread_wait(th[i]);
2241 extern "C" void *cftrec1_th(
void *p)
2243 int cfttree(
int n,
int j,
int k,
double *a,
int nw,
double *w);
2244 void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w);
2245 void cftmdl1(
int n,
double *a,
double *w);
2246 int isplt, j, k, m, n, n0, nw;
2249 n0 = ((cdft_arg_t *) p)->n0;
2250 n = ((cdft_arg_t *) p)->n;
2251 a = ((cdft_arg_t *) p)->a;
2252 nw = ((cdft_arg_t *) p)->nw;
2253 w = ((cdft_arg_t *) p)->w;
2257 cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2259 cftleaf(m, 1, &a[n - m], nw, w);
2261 for (j = n - m; j > 0; j -= m) {
2263 isplt = cfttree(m, j, k, a, nw, w);
2264 cftleaf(m, isplt, &a[j - m], nw, w);
2270 extern "C" void *cftrec2_th(
void *p)
2272 int cfttree(
int n,
int j,
int k,
double *a,
int nw,
double *w);
2273 void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w);
2274 void cftmdl2(
int n,
double *a,
double *w);
2275 int isplt, j, k, m, n, n0, nw;
2278 n0 = ((cdft_arg_t *) p)->n0;
2279 n = ((cdft_arg_t *) p)->n;
2280 a = ((cdft_arg_t *) p)->a;
2281 nw = ((cdft_arg_t *) p)->nw;
2282 w = ((cdft_arg_t *) p)->w;
2288 cftmdl2(m, &a[n - m], &w[nw - m]);
2290 cftleaf(m, 0, &a[n - m], nw, w);
2292 for (j = n - m; j > 0; j -= m) {
2294 isplt = cfttree(m, j, k, a, nw, w);
2295 cftleaf(m, isplt, &a[j - m], nw, w);
2302 extern "C" void cftrec4(
int n,
double *a,
int nw,
double *w)
2304 int cfttree(
int n,
int j,
int k,
double *a,
int nw,
double *w);
2305 void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w);
2306 void cftmdl1(
int n,
double *a,
double *w);
2312 cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2314 cftleaf(m, 1, &a[n - m], nw, w);
2316 for (j = n - m; j > 0; j -= m) {
2318 isplt = cfttree(m, j, k, a, nw, w);
2319 cftleaf(m, isplt, &a[j - m], nw, w);
2324 extern "C" int cfttree(
int n,
int j,
int k,
double *a,
int nw,
double *w)
2326 void cftmdl1(
int n,
double *a,
double *w);
2327 void cftmdl2(
int n,
double *a,
double *w);
2333 cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
2335 cftmdl2(n, &a[j - n], &w[nw - n]);
2339 for (i = k; (i & 3) == 0; i >>= 2) {
2345 cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
2350 cftmdl2(m, &a[j - m], &w[nw - m]);
2359 extern "C" void cftleaf(
int n,
int isplt,
double *a,
int nw,
double *w)
2361 void cftmdl1(
int n,
double *a,
double *w);
2362 void cftmdl2(
int n,
double *a,
double *w);
2363 void cftf161(
double *a,
double *w);
2364 void cftf162(
double *a,
double *w);
2365 void cftf081(
double *a,
double *w);
2366 void cftf082(
double *a,
double *w);
2369 cftmdl1(128, a, &w[nw - 64]);
2370 cftf161(a, &w[nw - 8]);
2371 cftf162(&a[32], &w[nw - 32]);
2372 cftf161(&a[64], &w[nw - 8]);
2373 cftf161(&a[96], &w[nw - 8]);
2374 cftmdl2(128, &a[128], &w[nw - 128]);
2375 cftf161(&a[128], &w[nw - 8]);
2376 cftf162(&a[160], &w[nw - 32]);
2377 cftf161(&a[192], &w[nw - 8]);
2378 cftf162(&a[224], &w[nw - 32]);
2379 cftmdl1(128, &a[256], &w[nw - 64]);
2380 cftf161(&a[256], &w[nw - 8]);
2381 cftf162(&a[288], &w[nw - 32]);
2382 cftf161(&a[320], &w[nw - 8]);
2383 cftf161(&a[352], &w[nw - 8]);
2385 cftmdl1(128, &a[384], &w[nw - 64]);
2386 cftf161(&a[480], &w[nw - 8]);
2388 cftmdl2(128, &a[384], &w[nw - 128]);
2389 cftf162(&a[480], &w[nw - 32]);
2391 cftf161(&a[384], &w[nw - 8]);
2392 cftf162(&a[416], &w[nw - 32]);
2393 cftf161(&a[448], &w[nw - 8]);
2395 cftmdl1(64, a, &w[nw - 32]);
2396 cftf081(a, &w[nw - 8]);
2397 cftf082(&a[16], &w[nw - 8]);
2398 cftf081(&a[32], &w[nw - 8]);
2399 cftf081(&a[48], &w[nw - 8]);
2400 cftmdl2(64, &a[64], &w[nw - 64]);
2401 cftf081(&a[64], &w[nw - 8]);
2402 cftf082(&a[80], &w[nw - 8]);
2403 cftf081(&a[96], &w[nw - 8]);
2404 cftf082(&a[112], &w[nw - 8]);
2405 cftmdl1(64, &a[128], &w[nw - 32]);
2406 cftf081(&a[128], &w[nw - 8]);
2407 cftf082(&a[144], &w[nw - 8]);
2408 cftf081(&a[160], &w[nw - 8]);
2409 cftf081(&a[176], &w[nw - 8]);
2411 cftmdl1(64, &a[192], &w[nw - 32]);
2412 cftf081(&a[240], &w[nw - 8]);
2414 cftmdl2(64, &a[192], &w[nw - 64]);
2415 cftf082(&a[240], &w[nw - 8]);
2417 cftf081(&a[192], &w[nw - 8]);
2418 cftf082(&a[208], &w[nw - 8]);
2419 cftf081(&a[224], &w[nw - 8]);
2424 extern "C" void cftmdl1(
int n,
double *a,
double *w)
2426 int j, j0, j1, j2, j3, k, m, mh;
2427 double wn4r, wk1r, wk1i, wk3r, wk3i;
2428 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2436 x0i = a[1] + a[j2 + 1];
2438 x1i = a[1] - a[j2 + 1];
2439 x2r = a[j1] + a[j3];
2440 x2i = a[j1 + 1] + a[j3 + 1];
2441 x3r = a[j1] - a[j3];
2442 x3i = a[j1 + 1] - a[j3 + 1];
2446 a[j1 + 1] = x0i - x2i;
2448 a[j2 + 1] = x1i + x3r;
2450 a[j3 + 1] = x1i - x3r;
2453 for (j = 2; j < mh; j += 2) {
2463 x0i = a[j + 1] + a[j2 + 1];
2465 x1i = a[j + 1] - a[j2 + 1];
2466 x2r = a[j1] + a[j3];
2467 x2i = a[j1 + 1] + a[j3 + 1];
2468 x3r = a[j1] - a[j3];
2469 x3i = a[j1 + 1] - a[j3 + 1];
2471 a[j + 1] = x0i + x2i;
2473 a[j1 + 1] = x0i - x2i;
2476 a[j2] = wk1r * x0r - wk1i * x0i;
2477 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2480 a[j3] = wk3r * x0r + wk3i * x0i;
2481 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2486 x0r = a[j0] + a[j2];
2487 x0i = a[j0 + 1] + a[j2 + 1];
2488 x1r = a[j0] - a[j2];
2489 x1i = a[j0 + 1] - a[j2 + 1];
2490 x2r = a[j1] + a[j3];
2491 x2i = a[j1 + 1] + a[j3 + 1];
2492 x3r = a[j1] - a[j3];
2493 x3i = a[j1 + 1] - a[j3 + 1];
2495 a[j0 + 1] = x0i + x2i;
2497 a[j1 + 1] = x0i - x2i;
2500 a[j2] = wk1i * x0r - wk1r * x0i;
2501 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2504 a[j3] = wk3i * x0r + wk3r * x0i;
2505 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2511 x0r = a[j0] + a[j2];
2512 x0i = a[j0 + 1] + a[j2 + 1];
2513 x1r = a[j0] - a[j2];
2514 x1i = a[j0 + 1] - a[j2 + 1];
2515 x2r = a[j1] + a[j3];
2516 x2i = a[j1 + 1] + a[j3 + 1];
2517 x3r = a[j1] - a[j3];
2518 x3i = a[j1 + 1] - a[j3 + 1];
2520 a[j0 + 1] = x0i + x2i;
2522 a[j1 + 1] = x0i - x2i;
2525 a[j2] = wn4r * (x0r - x0i);
2526 a[j2 + 1] = wn4r * (x0i + x0r);
2529 a[j3] = -wn4r * (x0r + x0i);
2530 a[j3 + 1] = -wn4r * (x0i - x0r);
2534 extern "C" void cftmdl2(
int n,
double *a,
double *w)
2536 int j, j0, j1, j2, j3, k, kr, m, mh;
2537 double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2538 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2546 x0r = a[0] - a[j2 + 1];
2548 x1r = a[0] + a[j2 + 1];
2550 x2r = a[j1] - a[j3 + 1];
2551 x2i = a[j1 + 1] + a[j3];
2552 x3r = a[j1] + a[j3 + 1];
2553 x3i = a[j1 + 1] - a[j3];
2554 y0r = wn4r * (x2r - x2i);
2555 y0i = wn4r * (x2i + x2r);
2559 a[j1 + 1] = x0i - y0i;
2560 y0r = wn4r * (x3r - x3i);
2561 y0i = wn4r * (x3i + x3r);
2563 a[j2 + 1] = x1i + y0r;
2565 a[j3 + 1] = x1i - y0r;
2568 for (j = 2; j < mh; j += 2) {
2582 x0r = a[j] - a[j2 + 1];
2583 x0i = a[j + 1] + a[j2];
2584 x1r = a[j] + a[j2 + 1];
2585 x1i = a[j + 1] - a[j2];
2586 x2r = a[j1] - a[j3 + 1];
2587 x2i = a[j1 + 1] + a[j3];
2588 x3r = a[j1] + a[j3 + 1];
2589 x3i = a[j1 + 1] - a[j3];
2590 y0r = wk1r * x0r - wk1i * x0i;
2591 y0i = wk1r * x0i + wk1i * x0r;
2592 y2r = wd1r * x2r - wd1i * x2i;
2593 y2i = wd1r * x2i + wd1i * x2r;
2595 a[j + 1] = y0i + y2i;
2597 a[j1 + 1] = y0i - y2i;
2598 y0r = wk3r * x1r + wk3i * x1i;
2599 y0i = wk3r * x1i - wk3i * x1r;
2600 y2r = wd3r * x3r + wd3i * x3i;
2601 y2i = wd3r * x3i - wd3i * x3r;
2603 a[j2 + 1] = y0i + y2i;
2605 a[j3 + 1] = y0i - y2i;
2610 x0r = a[j0] - a[j2 + 1];
2611 x0i = a[j0 + 1] + a[j2];
2612 x1r = a[j0] + a[j2 + 1];
2613 x1i = a[j0 + 1] - a[j2];
2614 x2r = a[j1] - a[j3 + 1];
2615 x2i = a[j1 + 1] + a[j3];
2616 x3r = a[j1] + a[j3 + 1];
2617 x3i = a[j1 + 1] - a[j3];
2618 y0r = wd1i * x0r - wd1r * x0i;
2619 y0i = wd1i * x0i + wd1r * x0r;
2620 y2r = wk1i * x2r - wk1r * x2i;
2621 y2i = wk1i * x2i + wk1r * x2r;
2623 a[j0 + 1] = y0i + y2i;
2625 a[j1 + 1] = y0i - y2i;
2626 y0r = wd3i * x1r + wd3r * x1i;
2627 y0i = wd3i * x1i - wd3r * x1r;
2628 y2r = wk3i * x3r + wk3r * x3i;
2629 y2i = wk3i * x3i - wk3r * x3r;
2631 a[j2 + 1] = y0i + y2i;
2633 a[j3 + 1] = y0i - y2i;
2641 x0r = a[j0] - a[j2 + 1];
2642 x0i = a[j0 + 1] + a[j2];
2643 x1r = a[j0] + a[j2 + 1];
2644 x1i = a[j0 + 1] - a[j2];
2645 x2r = a[j1] - a[j3 + 1];
2646 x2i = a[j1 + 1] + a[j3];
2647 x3r = a[j1] + a[j3 + 1];
2648 x3i = a[j1 + 1] - a[j3];
2649 y0r = wk1r * x0r - wk1i * x0i;
2650 y0i = wk1r * x0i + wk1i * x0r;
2651 y2r = wk1i * x2r - wk1r * x2i;
2652 y2i = wk1i * x2i + wk1r * x2r;
2654 a[j0 + 1] = y0i + y2i;
2656 a[j1 + 1] = y0i - y2i;
2657 y0r = wk1i * x1r - wk1r * x1i;
2658 y0i = wk1i * x1i + wk1r * x1r;
2659 y2r = wk1r * x3r - wk1i * x3i;
2660 y2i = wk1r * x3i + wk1i * x3r;
2662 a[j2 + 1] = y0i - y2i;
2664 a[j3 + 1] = y0i + y2i;
2668 extern "C" void cftfx41(
int n,
double *a,
int nw,
double *w)
2670 void cftf161(
double *a,
double *w);
2671 void cftf162(
double *a,
double *w);
2672 void cftf081(
double *a,
double *w);
2673 void cftf082(
double *a,
double *w);
2676 cftf161(a, &w[nw - 8]);
2677 cftf162(&a[32], &w[nw - 32]);
2678 cftf161(&a[64], &w[nw - 8]);
2679 cftf161(&a[96], &w[nw - 8]);
2681 cftf081(a, &w[nw - 8]);
2682 cftf082(&a[16], &w[nw - 8]);
2683 cftf081(&a[32], &w[nw - 8]);
2684 cftf081(&a[48], &w[nw - 8]);
2689 extern "C" void cftf161(
double *a,
double *w)
2691 double wn4r, wk1r, wk1i,
2692 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2693 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2694 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2695 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2696 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2721 x2r = a[10] + a[26];
2722 x2i = a[11] + a[27];
2723 x3r = a[10] - a[26];
2724 x3i = a[11] - a[27];
2731 y9r = wk1r * x0r - wk1i * x0i;
2732 y9i = wk1r * x0i + wk1i * x0r;
2735 y13r = wk1i * x0r - wk1r * x0i;
2736 y13i = wk1i * x0i + wk1r * x0r;
2741 x2r = a[12] + a[28];
2742 x2i = a[13] + a[29];
2743 x3r = a[12] - a[28];
2744 x3i = a[13] - a[29];
2751 y10r = wn4r * (x0r - x0i);
2752 y10i = wn4r * (x0i + x0r);
2755 y14r = wn4r * (x0r + x0i);
2756 y14i = wn4r * (x0i - x0r);
2761 x2r = a[14] + a[30];
2762 x2i = a[15] + a[31];
2763 x3r = a[14] - a[30];
2764 x3i = a[15] - a[31];
2771 y11r = wk1i * x0r - wk1r * x0i;
2772 y11i = wk1i * x0i + wk1r * x0r;
2775 y15r = wk1r * x0r - wk1i * x0i;
2776 y15i = wk1r * x0i + wk1i * x0r;
2811 x2r = wn4r * (x0r - x0i);
2812 x2i = wn4r * (x0i + x0r);
2815 x3r = wn4r * (x0r - x0i);
2816 x3i = wn4r * (x0i + x0r);
2848 extern "C" void cftf162(
double *a,
double *w)
2850 double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
2851 x0r, x0i, x1r, x1i, x2r, x2i,
2852 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2853 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2854 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2855 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2868 x2r = wn4r * (x0r - x0i);
2869 x2i = wn4r * (x0i + x0r);
2878 x2r = wn4r * (x0r - x0i);
2879 x2i = wn4r * (x0i + x0r);
2886 x1r = wk1r * x0r - wk1i * x0i;
2887 x1i = wk1r * x0i + wk1i * x0r;
2888 x0r = a[10] - a[27];
2889 x0i = a[11] + a[26];
2890 x2r = wk3i * x0r - wk3r * x0i;
2891 x2i = wk3i * x0i + wk3r * x0r;
2898 x1r = wk3r * x0r - wk3i * x0i;
2899 x1i = wk3r * x0i + wk3i * x0r;
2900 x0r = a[10] + a[27];
2901 x0i = a[11] - a[26];
2902 x2r = wk1r * x0r + wk1i * x0i;
2903 x2i = wk1r * x0i - wk1i * x0r;
2910 x1r = wk2r * x0r - wk2i * x0i;
2911 x1i = wk2r * x0i + wk2i * x0r;
2912 x0r = a[12] - a[29];
2913 x0i = a[13] + a[28];
2914 x2r = wk2i * x0r - wk2r * x0i;
2915 x2i = wk2i * x0i + wk2r * x0r;
2922 x1r = wk2i * x0r - wk2r * x0i;
2923 x1i = wk2i * x0i + wk2r * x0r;
2924 x0r = a[12] + a[29];
2925 x0i = a[13] - a[28];
2926 x2r = wk2r * x0r - wk2i * x0i;
2927 x2i = wk2r * x0i + wk2i * x0r;
2934 x1r = wk3r * x0r - wk3i * x0i;
2935 x1i = wk3r * x0i + wk3i * x0r;
2936 x0r = a[14] - a[31];
2937 x0i = a[15] + a[30];
2938 x2r = wk1i * x0r - wk1r * x0i;
2939 x2i = wk1i * x0i + wk1r * x0r;
2946 x1r = wk1i * x0r + wk1r * x0i;
2947 x1i = wk1i * x0i - wk1r * x0r;
2948 x0r = a[14] + a[31];
2949 x0i = a[15] - a[30];
2950 x2r = wk3i * x0r - wk3r * x0i;
2951 x2i = wk3i * x0i + wk3r * x0r;
2976 x2r = wn4r * (x0r - x0i);
2977 x2i = wn4r * (x0i + x0r);
2986 x2r = wn4r * (x0r - x0i);
2987 x2i = wn4r * (x0i + x0r);
3012 x2r = wn4r * (x0r - x0i);
3013 x2i = wn4r * (x0i + x0r);
3022 x2r = wn4r * (x0r - x0i);
3023 x2i = wn4r * (x0i + x0r);
3031 extern "C" void cftf081(
double *a,
double *w)
3033 double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
3034 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3035 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3070 y5r = wn4r * (x0r - x0i);
3071 y5i = wn4r * (x0r + x0i);
3072 y7r = wn4r * (x2r - x2i);
3073 y7i = wn4r * (x2r + x2i);
3093 extern "C" void cftf082(
double *a,
double *w)
3095 double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
3096 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3097 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3108 y2r = wn4r * (x0r - x0i);
3109 y2i = wn4r * (x0i + x0r);
3112 y3r = wn4r * (x0r - x0i);
3113 y3i = wn4r * (x0i + x0r);
3116 y4r = wk1r * x0r - wk1i * x0i;
3117 y4i = wk1r * x0i + wk1i * x0r;
3120 y5r = wk1i * x0r - wk1r * x0i;
3121 y5i = wk1i * x0i + wk1r * x0r;
3124 y6r = wk1i * x0r - wk1r * x0i;
3125 y6i = wk1i * x0i + wk1r * x0r;
3128 y7r = wk1r * x0r - wk1i * x0i;
3129 y7i = wk1r * x0i + wk1i * x0r;
3165 extern "C" void cftf040(
double *a)
3167 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3188 extern "C" void cftb040(
double *a)
3190 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3211 extern "C" void cftx020(
double *a)
3224 extern "C" void rftfsub(
int n,
double *a,
int nc,
double *c)
3226 int j, k, kk, ks, m;
3227 double wkr, wki, xr, xi, yr, yi;
3232 for (j = 2; j < m; j += 2) {
3235 wkr = 0.5 - c[nc - kk];
3238 xi = a[j + 1] + a[k + 1];
3239 yr = wkr * xr - wki * xi;
3240 yi = wkr * xi + wki * xr;
3249 extern "C" void rftbsub(
int n,
double *a,
int nc,
double *c)
3251 int j, k, kk, ks, m;
3252 double wkr, wki, xr, xi, yr, yi;
3257 for (j = 2; j < m; j += 2) {
3260 wkr = 0.5 - c[nc - kk];
3263 xi = a[j + 1] + a[k + 1];
3264 yr = wkr * xr + wki * xi;
3265 yi = wkr * xi - wki * xr;
3274 extern "C" void dctsub(
int n,
double *a,
int nc,
double *c)
3276 int j, k, kk, ks, m;
3277 double wkr, wki, xr;
3282 for (j = 1; j < m; j++) {
3285 wkr = c[kk] - c[nc - kk];
3286 wki = c[kk] + c[nc - kk];
3287 xr = wki * a[j] - wkr * a[k];
3288 a[j] = wkr * a[j] + wki * a[k];
3295 extern "C" void dstsub(
int n,
double *a,
int nc,
double *c)
3297 int j, k, kk, ks, m;
3298 double wkr, wki, xr;
3303 for (j = 1; j < m; j++) {
3306 wkr = c[kk] - c[nc - kk];
3307 wki = c[kk] + c[nc - kk];
3308 xr = wki * a[k] - wkr * a[j];
3309 a[k] = wkr * a[k] + wki * a[j];