278 void cdft(
int n,
int isgn,
double *a,
int *ip,
double *w)
280 void makewt(
int nw,
int *ip,
double *w);
281 void bitrv2(
int n,
int *ip,
double *a);
282 void bitrv2conj(
int n,
int *ip,
double *a);
283 void cftfsub(
int n,
double *a,
double *w);
284 void cftbsub(
int n,
double *a,
double *w);
286 if (n > (ip[0] << 2)) {
287 makewt(n >> 2, ip, w);
291 bitrv2(n, ip + 2, a);
294 bitrv2conj(n, ip + 2, a);
303 void rdft(
int n,
int isgn,
double *a,
int *ip,
double *w)
305 void makewt(
int nw,
int *ip,
double *w);
306 void makect(
int nc,
int *ip,
double *c);
307 void bitrv2(
int n,
int *ip,
double *a);
308 void cftfsub(
int n,
double *a,
double *w);
309 void cftbsub(
int n,
double *a,
double *w);
310 void rftfsub(
int n,
double *a,
int nc,
double *c);
311 void rftbsub(
int n,
double *a,
int nc,
double *c);
323 makect(nc, ip, w + nw);
327 bitrv2(n, ip + 2, a);
329 rftfsub(n, a, nc, w + nw);
337 a[1] = 0.5 * (a[0] - a[1]);
340 rftbsub(n, a, nc, w + nw);
341 bitrv2(n, ip + 2, a);
350 void ddct(
int n,
int isgn,
double *a,
int *ip,
double *w)
352 void makewt(
int nw,
int *ip,
double *w);
353 void makect(
int nc,
int *ip,
double *c);
354 void bitrv2(
int n,
int *ip,
double *a);
355 void cftfsub(
int n,
double *a,
double *w);
356 void cftbsub(
int n,
double *a,
double *w);
357 void rftfsub(
int n,
double *a,
int nc,
double *c);
358 void rftbsub(
int n,
double *a,
int nc,
double *c);
359 void dctsub(
int n,
double *a,
int nc,
double *c);
371 makect(nc, ip, w + nw);
375 for (j = n - 2; j >= 2; j -= 2) {
376 a[j + 1] = a[j] - a[j - 1];
382 rftbsub(n, a, nc, w + nw);
383 bitrv2(n, ip + 2, a);
389 dctsub(n, a, nc, w + nw);
392 bitrv2(n, ip + 2, a);
394 rftfsub(n, a, nc, w + nw);
400 for (j = 2; j < n; j += 2) {
401 a[j - 1] = a[j] - a[j + 1];
409 void ddst(
int n,
int isgn,
double *a,
int *ip,
double *w)
411 void makewt(
int nw,
int *ip,
double *w);
412 void makect(
int nc,
int *ip,
double *c);
413 void bitrv2(
int n,
int *ip,
double *a);
414 void cftfsub(
int n,
double *a,
double *w);
415 void cftbsub(
int n,
double *a,
double *w);
416 void rftfsub(
int n,
double *a,
int nc,
double *c);
417 void rftbsub(
int n,
double *a,
int nc,
double *c);
418 void dstsub(
int n,
double *a,
int nc,
double *c);
430 makect(nc, ip, w + nw);
434 for (j = n - 2; j >= 2; j -= 2) {
435 a[j + 1] = -a[j] - a[j - 1];
441 rftbsub(n, a, nc, w + nw);
442 bitrv2(n, ip + 2, a);
448 dstsub(n, a, nc, w + nw);
451 bitrv2(n, ip + 2, a);
453 rftfsub(n, a, nc, w + nw);
459 for (j = 2; j < n; j += 2) {
460 a[j - 1] = -a[j] - a[j + 1];
468 void dfct(
int n,
double *a,
double *t,
int *ip,
double *w)
470 void makewt(
int nw,
int *ip,
double *w);
471 void makect(
int nc,
int *ip,
double *c);
472 void bitrv2(
int n,
int *ip,
double *a);
473 void cftfsub(
int n,
double *a,
double *w);
474 void rftfsub(
int n,
double *a,
int nc,
double *c);
475 void dctsub(
int n,
double *a,
int nc,
double *c);
476 int j, k, l, m, mh, nw, nc;
477 double xr, xi, yr, yi;
487 makect(nc, ip, w + nw);
497 for (j = 1; j < mh; j++) {
499 xr = a[j] - a[n - j];
500 xi = a[j] + a[n - j];
501 yr = a[k] - a[n - k];
502 yi = a[k] + a[n - k];
508 t[mh] = a[mh] + a[n - mh];
510 dctsub(m, a, nc, w + nw);
512 bitrv2(m, ip + 2, a);
514 rftfsub(m, a, nc, w + nw);
518 a[n - 1] = a[0] - a[1];
520 for (j = m - 2; j >= 2; j -= 2) {
521 a[2 * j + 1] = a[j] + a[j + 1];
522 a[2 * j - 1] = a[j] - a[j + 1];
527 dctsub(m, t, nc, w + nw);
529 bitrv2(m, ip + 2, t);
531 rftfsub(m, t, nc, w + nw);
535 a[n - l] = t[0] - t[1];
538 for (j = 2; j < m; j += 2) {
540 a[k - l] = t[j] - t[j + 1];
541 a[k + l] = t[j] + t[j + 1];
545 for (j = 0; j < mh; j++) {
547 t[j] = t[m + k] - t[m + j];
548 t[k] = t[m + k] + t[m + j];
564 void dfst(
int n,
double *a,
double *t,
int *ip,
double *w)
566 void makewt(
int nw,
int *ip,
double *w);
567 void makect(
int nc,
int *ip,
double *c);
568 void bitrv2(
int n,
int *ip,
double *a);
569 void cftfsub(
int n,
double *a,
double *w);
570 void rftfsub(
int n,
double *a,
int nc,
double *c);
571 void dstsub(
int n,
double *a,
int nc,
double *c);
572 int j, k, l, m, mh, nw, nc;
573 double xr, xi, yr, yi;
583 makect(nc, ip, w + nw);
588 for (j = 1; j < mh; j++) {
590 xr = a[j] + a[n - j];
591 xi = a[j] - a[n - j];
592 yr = a[k] + a[n - k];
593 yi = a[k] - a[n - k];
599 t[0] = a[mh] - a[n - mh];
602 dstsub(m, a, nc, w + nw);
604 bitrv2(m, ip + 2, a);
606 rftfsub(m, a, nc, w + nw);
610 a[n - 1] = a[1] - a[0];
612 for (j = m - 2; j >= 2; j -= 2) {
613 a[2 * j + 1] = a[j] - a[j + 1];
614 a[2 * j - 1] = -a[j] - a[j + 1];
619 dstsub(m, t, nc, w + nw);
621 bitrv2(m, ip + 2, t);
623 rftfsub(m, t, nc, w + nw);
627 a[n - l] = t[1] - t[0];
630 for (j = 2; j < m; j += 2) {
632 a[k - l] = -t[j] - t[j + 1];
633 a[k + l] = t[j] - t[j + 1];
637 for (j = 1; j < mh; j++) {
639 t[j] = t[m + k] + t[m + j];
640 t[k] = t[m + k] - t[m + j];
656 void makewt(
int nw,
int *ip,
double *w)
658 void bitrv2(
int n,
int *ip,
double *a);
666 delta = atan(1.0) / nwh;
669 w[nwh] = cos(delta * nwh);
672 for (j = 2; j < nwh; j += 2) {
680 for (j = nwh - 2; j >= 2; j -= 2) {
686 bitrv2(nw, ip + 2, w);
692 void makect(
int nc,
int *ip,
double *c)
700 delta = atan(1.0) / nch;
701 c[0] = cos(delta * nch);
703 for (j = 1; j < nch; j++) {
704 c[j] = 0.5 * cos(delta * j);
705 c[nc - j] = 0.5 * sin(delta * j);
714 void bitrv2(
int n,
int *ip,
double *a)
716 int j, j1, k, k1, l, m, m2;
717 double xr, xi, yr, yi;
722 while ((m << 3) < l) {
724 for (j = 0; j < m; j++) {
725 ip[m + j] = ip[j] + l;
731 for (k = 0; k < m; k++) {
732 for (j = 0; j < k; j++) {
774 j1 = 2 * k + m2 + ip[k];
786 for (k = 1; k < m; k++) {
787 for (j = 0; j < k; j++) {
814 void bitrv2conj(
int n,
int *ip,
double *a)
816 int j, j1, k, k1, l, m, m2;
817 double xr, xi, yr, yi;
822 while ((m << 3) < l) {
824 for (j = 0; j < m; j++) {
825 ip[m + j] = ip[j] + l;
831 for (k = 0; k < m; k++) {
832 for (j = 0; j < k; j++) {
875 a[k1 + 1] = -a[k1 + 1];
887 a[k1 + 1] = -a[k1 + 1];
891 a[m2 + 1] = -a[m2 + 1];
892 for (k = 1; k < m; k++) {
893 for (j = 0; j < k; j++) {
916 a[k1 + 1] = -a[k1 + 1];
917 a[k1 + m2 + 1] = -a[k1 + m2 + 1];
923 void cftfsub(
int n,
double *a,
double *w)
925 void cft1st(
int n,
double *a,
double *w);
926 void cftmdl(
int n,
int l,
double *a,
double *w);
927 int j, j1, j2, j3, l;
928 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
934 while ((l << 3) <= n) {
940 for (j = 0; j < l; j += 2) {
945 x0i = a[j + 1] + a[j1 + 1];
947 x1i = a[j + 1] - a[j1 + 1];
949 x2i = a[j2 + 1] + a[j3 + 1];
951 x3i = a[j2 + 1] - a[j3 + 1];
953 a[j + 1] = x0i + x2i;
955 a[j2 + 1] = x0i - x2i;
957 a[j1 + 1] = x1i + x3r;
959 a[j3 + 1] = x1i - x3r;
961 }
else if ((l << 1) == n) {
962 for (j = 0; j < l; j += 2) {
965 x0i = a[j + 1] - a[j1 + 1];
967 a[j + 1] += a[j1 + 1];
975 void cftbsub(
int n,
double *a,
double *w)
977 void cft1st(
int n,
double *a,
double *w);
978 void cftmdl(
int n,
int l,
double *a,
double *w);
979 int j, j1, j2, j3, j4, j5, j6, j7, l;
980 double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
981 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
982 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
988 while ((l << 3) < n) {
995 for (j = 0; j < l; j += 2) {
1004 x0i = -a[j + 1] - a[j1 + 1];
1006 x1i = -a[j + 1] + a[j1 + 1];
1007 x2r = a[j2] + a[j3];
1008 x2i = a[j2 + 1] + a[j3 + 1];
1009 x3r = a[j2] - a[j3];
1010 x3i = a[j2 + 1] - a[j3 + 1];
1019 x0r = a[j4] + a[j5];
1020 x0i = a[j4 + 1] + a[j5 + 1];
1021 x1r = a[j4] - a[j5];
1022 x1i = a[j4 + 1] - a[j5 + 1];
1023 x2r = a[j6] + a[j7];
1024 x2i = a[j6 + 1] + a[j7 + 1];
1025 x3r = a[j6] - a[j7];
1026 x3i = a[j6 + 1] - a[j7 + 1];
1035 y5r = wn4r * (x0r - x0i);
1036 y5i = wn4r * (x0r + x0i);
1037 y7r = wn4r * (x2r - x2i);
1038 y7i = wn4r * (x2r + x2i);
1040 a[j1 + 1] = y1i - y5i;
1042 a[j5 + 1] = y1i + y5i;
1044 a[j3 + 1] = y3i - y7r;
1046 a[j7 + 1] = y3i + y7r;
1048 a[j + 1] = y0i - y4i;
1050 a[j4 + 1] = y0i + y4i;
1052 a[j2 + 1] = y2i - y6r;
1054 a[j6 + 1] = y2i + y6r;
1056 }
else if ((l << 2) == n) {
1057 for (j = 0; j < l; j += 2) {
1062 x0i = -a[j + 1] - a[j1 + 1];
1064 x1i = -a[j + 1] + a[j1 + 1];
1065 x2r = a[j2] + a[j3];
1066 x2i = a[j2 + 1] + a[j3 + 1];
1067 x3r = a[j2] - a[j3];
1068 x3i = a[j2 + 1] - a[j3 + 1];
1070 a[j + 1] = x0i - x2i;
1072 a[j2 + 1] = x0i + x2i;
1074 a[j1 + 1] = x1i - x3r;
1076 a[j3 + 1] = x1i + x3r;
1079 for (j = 0; j < l; j += 2) {
1082 x0i = -a[j + 1] + a[j1 + 1];
1084 a[j + 1] = -a[j + 1] - a[j1 + 1];
1092 void cft1st(
int n,
double *a,
double *w)
1095 double wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
1096 wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1097 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1098 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
1099 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1122 x2r = a[12] + a[14];
1123 x2i = a[13] + a[15];
1124 x3r = a[12] - a[14];
1125 x3i = a[13] - a[15];
1134 y5r = wn4r * (x0r - x0i);
1135 y5i = wn4r * (x0r + x0i);
1136 y7r = wn4r * (x2r - x2i);
1137 y7i = wn4r * (x2r + x2i);
1157 x0r = a[16] + a[18];
1158 x0i = a[17] + a[19];
1159 x1r = a[16] - a[18];
1160 x1i = a[17] - a[19];
1161 x2r = a[20] + a[22];
1162 x2i = a[21] + a[23];
1163 x3r = a[20] - a[22];
1164 x3i = a[21] - a[23];
1173 x0r = a[24] + a[26];
1174 x0i = a[25] + a[27];
1175 x1r = a[24] - a[26];
1176 x1i = a[25] - a[27];
1177 x2r = a[28] + a[30];
1178 x2i = a[29] + a[31];
1179 x3r = a[28] - a[30];
1180 x3i = a[29] - a[31];
1189 y5r = wk1i * x0r - wk1r * x0i;
1190 y5i = wk1i * x0i + wk1r * x0r;
1191 y7r = wk1r * x2r + wk1i * x2i;
1192 y7i = wk1r * x2i - wk1i * x2r;
1193 x0r = wk1r * y1r - wk1i * y1i;
1194 x0i = wk1r * y1i + wk1i * y1r;
1199 x0r = wk1i * y3r - wk1r * y3i;
1200 x0i = wk1i * y3i + wk1r * y3r;
1211 a[20] = wn4r * (x0r - x0i);
1212 a[21] = wn4r * (x0i + x0r);
1215 a[28] = wn4r * (x0r - x0i);
1216 a[29] = wn4r * (x0i + x0r);
1218 for (j = 32; j < n; j += 16) {
1225 wk3r = wk1r - wtmp * wk1i;
1226 wk3i = wtmp * wk1r - wk1i;
1227 wk4r = 1 - wtmp * wk2i;
1230 wk5r = wk3r - wtmp * wk1i;
1231 wk5i = wtmp * wk1r - wk3i;
1232 wk6r = wk2r - wtmp * wk2i;
1233 wk6i = wtmp * wk2r - wk2i;
1234 wk7r = wk1r - wtmp * wk3i;
1235 wk7i = wtmp * wk3r - wk1i;
1236 x0r = a[j] + a[j + 2];
1237 x0i = a[j + 1] + a[j + 3];
1238 x1r = a[j] - a[j + 2];
1239 x1i = a[j + 1] - a[j + 3];
1240 x2r = a[j + 4] + a[j + 6];
1241 x2i = a[j + 5] + a[j + 7];
1242 x3r = a[j + 4] - a[j + 6];
1243 x3i = a[j + 5] - a[j + 7];
1252 x0r = a[j + 8] + a[j + 10];
1253 x0i = a[j + 9] + a[j + 11];
1254 x1r = a[j + 8] - a[j + 10];
1255 x1i = a[j + 9] - a[j + 11];
1256 x2r = a[j + 12] + a[j + 14];
1257 x2i = a[j + 13] + a[j + 15];
1258 x3r = a[j + 12] - a[j + 14];
1259 x3i = a[j + 13] - a[j + 15];
1268 y5r = wn4r * (x0r - x0i);
1269 y5i = wn4r * (x0r + x0i);
1270 y7r = wn4r * (x2r - x2i);
1271 y7i = wn4r * (x2r + x2i);
1274 a[j + 2] = wk1r * x0r - wk1i * x0i;
1275 a[j + 3] = wk1r * x0i + wk1i * x0r;
1278 a[j + 10] = wk5r * x0r - wk5i * x0i;
1279 a[j + 11] = wk5r * x0i + wk5i * x0r;
1282 a[j + 6] = wk3r * x0r - wk3i * x0i;
1283 a[j + 7] = wk3r * x0i + wk3i * x0r;
1286 a[j + 14] = wk7r * x0r - wk7i * x0i;
1287 a[j + 15] = wk7r * x0i + wk7i * x0r;
1289 a[j + 1] = y0i + y4i;
1292 a[j + 8] = wk4r * x0r - wk4i * x0i;
1293 a[j + 9] = wk4r * x0i + wk4i * x0r;
1296 a[j + 4] = wk2r * x0r - wk2i * x0i;
1297 a[j + 5] = wk2r * x0i + wk2i * x0r;
1300 a[j + 12] = wk6r * x0r - wk6i * x0i;
1301 a[j + 13] = wk6r * x0i + wk6i * x0r;
1307 void cftmdl(
int n,
int l,
double *a,
double *w)
1309 int j, j1, j2, j3, j4, j5, j6, j7, k, k1, m;
1310 double wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
1311 wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1312 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1313 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
1314 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1318 for (j = 0; j < l; j += 2) {
1327 x0i = a[j + 1] + a[j1 + 1];
1329 x1i = a[j + 1] - a[j1 + 1];
1330 x2r = a[j2] + a[j3];
1331 x2i = a[j2 + 1] + a[j3 + 1];
1332 x3r = a[j2] - a[j3];
1333 x3i = a[j2 + 1] - a[j3 + 1];
1342 x0r = a[j4] + a[j5];
1343 x0i = a[j4 + 1] + a[j5 + 1];
1344 x1r = a[j4] - a[j5];
1345 x1i = a[j4 + 1] - a[j5 + 1];
1346 x2r = a[j6] + a[j7];
1347 x2i = a[j6 + 1] + a[j7 + 1];
1348 x3r = a[j6] - a[j7];
1349 x3i = a[j6 + 1] - a[j7 + 1];
1358 y5r = wn4r * (x0r - x0i);
1359 y5i = wn4r * (x0r + x0i);
1360 y7r = wn4r * (x2r - x2i);
1361 y7i = wn4r * (x2r + x2i);
1363 a[j1 + 1] = y1i + y5i;
1365 a[j5 + 1] = y1i - y5i;
1367 a[j3 + 1] = y3i + y7r;
1369 a[j7 + 1] = y3i - y7r;
1371 a[j + 1] = y0i + y4i;
1373 a[j4 + 1] = y0i - y4i;
1375 a[j2 + 1] = y2i + y6r;
1377 a[j6 + 1] = y2i - y6r;
1382 for (j = m; j < l + m; j += 2) {
1391 x0i = a[j + 1] + a[j1 + 1];
1393 x1i = a[j + 1] - a[j1 + 1];
1394 x2r = a[j2] + a[j3];
1395 x2i = a[j2 + 1] + a[j3 + 1];
1396 x3r = a[j2] - a[j3];
1397 x3i = a[j2 + 1] - a[j3 + 1];
1406 x0r = a[j4] + a[j5];
1407 x0i = a[j4 + 1] + a[j5 + 1];
1408 x1r = a[j4] - a[j5];
1409 x1i = a[j4 + 1] - a[j5 + 1];
1410 x2r = a[j6] + a[j7];
1411 x2i = a[j6 + 1] + a[j7 + 1];
1412 x3r = a[j6] - a[j7];
1413 x3i = a[j6 + 1] - a[j7 + 1];
1422 y5r = wk1i * x0r - wk1r * x0i;
1423 y5i = wk1i * x0i + wk1r * x0r;
1424 y7r = wk1r * x2r + wk1i * x2i;
1425 y7i = wk1r * x2i - wk1i * x2r;
1426 x0r = wk1r * y1r - wk1i * y1i;
1427 x0i = wk1r * y1i + wk1i * y1r;
1429 a[j1 + 1] = x0i + y5i;
1431 a[j5 + 1] = x0r - y5r;
1432 x0r = wk1i * y3r - wk1r * y3i;
1433 x0i = wk1i * y3i + wk1r * y3r;
1435 a[j3 + 1] = x0i + y7i;
1437 a[j7 + 1] = x0r + y7r;
1439 a[j + 1] = y0i + y4i;
1441 a[j4 + 1] = y0r - y4r;
1444 a[j2] = wn4r * (x0r - x0i);
1445 a[j2 + 1] = wn4r * (x0i + x0r);
1448 a[j6] = wn4r * (x0r - x0i);
1449 a[j6 + 1] = wn4r * (x0i + x0r);
1452 for (k = 2 * m; k < n; k += m) {
1459 wk3r = wk1r - wtmp * wk1i;
1460 wk3i = wtmp * wk1r - wk1i;
1461 wk4r = 1 - wtmp * wk2i;
1464 wk5r = wk3r - wtmp * wk1i;
1465 wk5i = wtmp * wk1r - wk3i;
1466 wk6r = wk2r - wtmp * wk2i;
1467 wk6i = wtmp * wk2r - wk2i;
1468 wk7r = wk1r - wtmp * wk3i;
1469 wk7i = wtmp * wk3r - wk1i;
1470 for (j = k; j < l + k; j += 2) {
1479 x0i = a[j + 1] + a[j1 + 1];
1481 x1i = a[j + 1] - a[j1 + 1];
1482 x2r = a[j2] + a[j3];
1483 x2i = a[j2 + 1] + a[j3 + 1];
1484 x3r = a[j2] - a[j3];
1485 x3i = a[j2 + 1] - a[j3 + 1];
1494 x0r = a[j4] + a[j5];
1495 x0i = a[j4 + 1] + a[j5 + 1];
1496 x1r = a[j4] - a[j5];
1497 x1i = a[j4 + 1] - a[j5 + 1];
1498 x2r = a[j6] + a[j7];
1499 x2i = a[j6 + 1] + a[j7 + 1];
1500 x3r = a[j6] - a[j7];
1501 x3i = a[j6 + 1] - a[j7 + 1];
1510 y5r = wn4r * (x0r - x0i);
1511 y5i = wn4r * (x0r + x0i);
1512 y7r = wn4r * (x2r - x2i);
1513 y7i = wn4r * (x2r + x2i);
1516 a[j1] = wk1r * x0r - wk1i * x0i;
1517 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1520 a[j5] = wk5r * x0r - wk5i * x0i;
1521 a[j5 + 1] = wk5r * x0i + wk5i * x0r;
1524 a[j3] = wk3r * x0r - wk3i * x0i;
1525 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1528 a[j7] = wk7r * x0r - wk7i * x0i;
1529 a[j7 + 1] = wk7r * x0i + wk7i * x0r;
1531 a[j + 1] = y0i + y4i;
1534 a[j4] = wk4r * x0r - wk4i * x0i;
1535 a[j4 + 1] = wk4r * x0i + wk4i * x0r;
1538 a[j2] = wk2r * x0r - wk2i * x0i;
1539 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
1542 a[j6] = wk6r * x0r - wk6i * x0i;
1543 a[j6 + 1] = wk6r * x0i + wk6i * x0r;
1550 void rftfsub(
int n,
double *a,
int nc,
double *c)
1552 int j, k, kk, ks, m;
1553 double wkr, wki, xr, xi, yr, yi;
1558 for (j = 2; j < m; j += 2) {
1561 wkr = 0.5 - c[nc - kk];
1564 xi = a[j + 1] + a[k + 1];
1565 yr = wkr * xr - wki * xi;
1566 yi = wkr * xi + wki * xr;
1575 void rftbsub(
int n,
double *a,
int nc,
double *c)
1577 int j, k, kk, ks, m;
1578 double wkr, wki, xr, xi, yr, yi;
1584 for (j = 2; j < m; j += 2) {
1587 wkr = 0.5 - c[nc - kk];
1590 xi = a[j + 1] + a[k + 1];
1591 yr = wkr * xr + wki * xi;
1592 yi = wkr * xi - wki * xr;
1594 a[j + 1] = yi - a[j + 1];
1596 a[k + 1] = yi - a[k + 1];
1598 a[m + 1] = -a[m + 1];
1602 void dctsub(
int n,
double *a,
int nc,
double *c)
1604 int j, k, kk, ks, m;
1605 double wkr, wki, xr;
1610 for (j = 1; j < m; j++) {
1613 wkr = c[kk] - c[nc - kk];
1614 wki = c[kk] + c[nc - kk];
1615 xr = wki * a[j] - wkr * a[k];
1616 a[j] = wkr * a[j] + wki * a[k];
1623 void dstsub(
int n,
double *a,
int nc,
double *c)
1625 int j, k, kk, ks, m;
1626 double wkr, wki, xr;
1631 for (j = 1; j < m; j++) {
1634 wkr = c[kk] - c[nc - kk];
1635 wki = c[kk] + c[nc - kk];
1636 xr = wki * a[k] - wkr * a[j];
1637 a[k] = wkr * a[k] + wki * a[j];