Lines Matching +full:n +full:- +full:1
1 /* Copyright Takuya OOURA, 1996-2001.
7 Package home: http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html
22 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
32 -------- Complex DFT (Discrete Fourier Transform) --------
35 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
37 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
38 (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
42 cdft(2*n, 1, a, ip, w);
45 cdft(2*n, -1, a, ip, w);
47 2*n :data length (int)
48 n >= 1, n = power of 2
49 a[0...2*n-1] :input/output data (double *)
52 a[2*j+1] = Im(x[j]), 0<=j<n
55 a[2*k+1] = Im(X[k]), 0<=k<n
57 length of ip >= 2+sqrt(n)
60 2+(1<<(int)(log(n+0.5)/log(2))/2).
61 ip[0],ip[1] are pointers of the cos/sin table.
62 w[0...n/2-1] :cos/sin table (double *)
66 cdft(2*n, -1, a, ip, w);
68 cdft(2*n, 1, a, ip, w);
69 for (j = 0; j <= 2 * n - 1; j++) {
70 a[j] *= 1.0 / n;
75 -------- Real DFT / Inverse of Real DFT --------
78 R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
79 I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
81 a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
82 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
83 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
87 rdft(n, 1, a, ip, w);
90 rdft(n, -1, a, ip, w);
92 n :data length (int)
93 n >= 2, n = power of 2
94 a[0...n-1] :input/output data (double *)
97 a[2*k] = R[k], 0<=k<n/2
98 a[2*k+1] = I[k], 0<k<n/2
99 a[1] = R[n/2]
102 a[2*j] = R[j], 0<=j<n/2
103 a[2*j+1] = I[j], 0<j<n/2
104 a[1] = R[n/2]
106 length of ip >= 2+sqrt(n/2)
109 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
110 ip[0],ip[1] are pointers of the cos/sin table.
111 w[0...n/2-1] :cos/sin table (double *)
115 rdft(n, 1, a, ip, w);
117 rdft(n, -1, a, ip, w);
118 for (j = 0; j <= n - 1; j++) {
119 a[j] *= 2.0 / n;
124 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
127 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
129 C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
133 ddct(n, 1, a, ip, w);
136 ddct(n, -1, a, ip, w);
138 n :data length (int)
139 n >= 2, n = power of 2
140 a[0...n-1] :input/output data (double *)
142 a[k] = C[k], 0<=k<n
144 length of ip >= 2+sqrt(n/2)
147 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
148 ip[0],ip[1] are pointers of the cos/sin table.
149 w[0...n*5/4-1] :cos/sin table (double *)
153 ddct(n, -1, a, ip, w);
156 ddct(n, 1, a, ip, w);
157 for (j = 0; j <= n - 1; j++) {
158 a[j] *= 2.0 / n;
163 -------- DST (Discrete Sine Transform) / Inverse of DST --------
166 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
168 S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
172 ddst(n, 1, a, ip, w);
175 ddst(n, -1, a, ip, w);
177 n :data length (int)
178 n >= 2, n = power of 2
179 a[0...n-1] :input/output data (double *)
182 a[j] = A[j], 0<j<n
183 a[0] = A[n]
185 a[k] = S[k], 0<=k<n
188 a[k] = S[k], 0<k<n
189 a[0] = S[n]
191 length of ip >= 2+sqrt(n/2)
194 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
195 ip[0],ip[1] are pointers of the cos/sin table.
196 w[0...n*5/4-1] :cos/sin table (double *)
200 ddst(n, -1, a, ip, w);
203 ddst(n, 1, a, ip, w);
204 for (j = 0; j <= n - 1; j++) {
205 a[j] *= 2.0 / n;
210 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
212 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
215 dfct(n, a, t, ip, w);
217 n :data length - 1 (int)
218 n >= 2, n = power of 2
219 a[0...n] :input/output data (double *)
221 a[k] = C[k], 0<=k<=n
222 t[0...n/2] :work area (double *)
224 length of ip >= 2+sqrt(n/4)
227 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
228 ip[0],ip[1] are pointers of the cos/sin table.
229 w[0...n*5/8-1] :cos/sin table (double *)
234 a[n] *= 0.5;
235 dfct(n, a, t, ip, w);
238 a[n] *= 0.5;
239 dfct(n, a, t, ip, w);
240 for (j = 0; j <= n; j++) {
241 a[j] *= 2.0 / n;
246 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
248 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
251 dfst(n, a, t, ip, w);
253 n :data length + 1 (int)
254 n >= 2, n = power of 2
255 a[0...n-1] :input/output data (double *)
257 a[k] = S[k], 0<k<n
259 t[0...n/2-1] :work area (double *)
261 length of ip >= 2+sqrt(n/4)
264 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
265 ip[0],ip[1] are pointers of the cos/sin table.
266 w[0...n*5/8-1] :cos/sin table (double *)
270 dfst(n, a, t, ip, w);
272 dfst(n, a, t, ip, w);
273 for (j = 1; j <= n - 1; j++) {
274 a[j] *= 2.0 / n;
309 static void bitrv2conj(int n, int *ip, double *a);
310 static void bitrv2(int n, int *ip, double *a);
311 static void cft1st(int n, double *a, double const *w);
312 static void cftbsub(int n, double *a, double const *w);
313 static void cftfsub(int n, double *a, double const *w);
314 static void cftmdl(int n, int l, double *a, double const *w);
315 static void dctsub(int n, double *a, int nc, double const *c);
316 static void dstsub(int n, double *a, int nc, double const *c);
319 static void rftbsub(int n, double *a, int nc, double const *c);
320 static void rftfsub(int n, double *a, int nc, double const *c);
323 void cdft(int n, int isgn, double *a, int *ip, double *w) in cdft() argument
325 if (n > FFT4G_MAX_SIZE) in cdft()
328 if (n > (ip[0] << 2)) { in cdft()
329 makewt(n >> 2, ip, w); in cdft()
331 if (n > 4) { in cdft()
333 bitrv2(n, ip + 2, a); in cdft()
334 cftfsub(n, a, w); in cdft()
336 bitrv2conj(n, ip + 2, a); in cdft()
337 cftbsub(n, a, w); in cdft()
339 } else if (n == 4) { in cdft()
340 cftfsub(n, a, w); in cdft()
345 void rdft(int n, int isgn, double *a, int *ip, double *w) in rdft() argument
350 if (n > FFT4G_MAX_SIZE) in rdft()
354 if (n > (nw << 2)) { in rdft()
355 nw = n >> 2; in rdft()
358 nc = ip[1]; in rdft()
359 if (n > (nc << 2)) { in rdft()
360 nc = n >> 2; in rdft()
364 if (n > 4) { in rdft()
365 bitrv2(n, ip + 2, a); in rdft()
366 cftfsub(n, a, w); in rdft()
367 rftfsub(n, a, nc, w + nw); in rdft()
368 } else if (n == 4) { in rdft()
369 cftfsub(n, a, w); in rdft()
371 xi = a[0] - a[1]; in rdft()
372 a[0] += a[1]; in rdft()
373 a[1] = xi; in rdft()
375 a[1] = 0.5 * (a[0] - a[1]); in rdft()
376 a[0] -= a[1]; in rdft()
377 if (n > 4) { in rdft()
378 rftbsub(n, a, nc, w + nw); in rdft()
379 bitrv2(n, ip + 2, a); in rdft()
380 cftbsub(n, a, w); in rdft()
381 } else if (n == 4) { in rdft()
382 cftfsub(n, a, w); in rdft()
388 void ddct(int n, int isgn, double *a, int *ip, double *w) in ddct() argument
393 if (n > FFT4G_MAX_SIZE) in ddct()
397 if (n > (nw << 2)) { in ddct()
398 nw = n >> 2; in ddct()
401 nc = ip[1]; in ddct()
402 if (n > nc) { in ddct()
403 nc = n; in ddct()
407 xr = a[n - 1]; in ddct()
408 for (j = n - 2; j >= 2; j -= 2) { in ddct()
409 a[j + 1] = a[j] - a[j - 1]; in ddct()
410 a[j] += a[j - 1]; in ddct()
412 a[1] = a[0] - xr; in ddct()
414 if (n > 4) { in ddct()
415 rftbsub(n, a, nc, w + nw); in ddct()
416 bitrv2(n, ip + 2, a); in ddct()
417 cftbsub(n, a, w); in ddct()
418 } else if (n == 4) { in ddct()
419 cftfsub(n, a, w); in ddct()
422 dctsub(n, a, nc, w + nw); in ddct()
424 if (n > 4) { in ddct()
425 bitrv2(n, ip + 2, a); in ddct()
426 cftfsub(n, a, w); in ddct()
427 rftfsub(n, a, nc, w + nw); in ddct()
428 } else if (n == 4) { in ddct()
429 cftfsub(n, a, w); in ddct()
431 xr = a[0] - a[1]; in ddct()
432 a[0] += a[1]; in ddct()
433 for (j = 2; j < n; j += 2) { in ddct()
434 a[j - 1] = a[j] - a[j + 1]; in ddct()
435 a[j] += a[j + 1]; in ddct()
437 a[n - 1] = xr; in ddct()
442 void ddst(int n, int isgn, double *a, int *ip, double *w) in ddst() argument
447 if (n > FFT4G_MAX_SIZE) in ddst()
451 if (n > (nw << 2)) { in ddst()
452 nw = n >> 2; in ddst()
455 nc = ip[1]; in ddst()
456 if (n > nc) { in ddst()
457 nc = n; in ddst()
461 xr = a[n - 1]; in ddst()
462 for (j = n - 2; j >= 2; j -= 2) { in ddst()
463 a[j + 1] = -a[j] - a[j - 1]; in ddst()
464 a[j] -= a[j - 1]; in ddst()
466 a[1] = a[0] + xr; in ddst()
467 a[0] -= xr; in ddst()
468 if (n > 4) { in ddst()
469 rftbsub(n, a, nc, w + nw); in ddst()
470 bitrv2(n, ip + 2, a); in ddst()
471 cftbsub(n, a, w); in ddst()
472 } else if (n == 4) { in ddst()
473 cftfsub(n, a, w); in ddst()
476 dstsub(n, a, nc, w + nw); in ddst()
478 if (n > 4) { in ddst()
479 bitrv2(n, ip + 2, a); in ddst()
480 cftfsub(n, a, w); in ddst()
481 rftfsub(n, a, nc, w + nw); in ddst()
482 } else if (n == 4) { in ddst()
483 cftfsub(n, a, w); in ddst()
485 xr = a[0] - a[1]; in ddst()
486 a[0] += a[1]; in ddst()
487 for (j = 2; j < n; j += 2) { in ddst()
488 a[j - 1] = -a[j] - a[j + 1]; in ddst()
489 a[j] -= a[j + 1]; in ddst()
491 a[n - 1] = -xr; in ddst()
496 void dfct(int n, double *a, double *t, int *ip, double *w) in dfct() argument
501 if (n > FFT4G_MAX_SIZE) in dfct()
505 if (n > (nw << 3)) { in dfct()
506 nw = n >> 3; in dfct()
509 nc = ip[1]; in dfct()
510 if (n > (nc << 1)) { in dfct()
511 nc = n >> 1; in dfct()
514 m = n >> 1; in dfct()
516 xi = a[0] + a[n]; in dfct()
517 a[0] -= a[n]; in dfct()
518 t[0] = xi - yi; in dfct()
520 if (n > 2) { in dfct()
521 mh = m >> 1; in dfct()
522 for (j = 1; j < mh; j++) { in dfct()
523 k = m - j; in dfct()
524 xr = a[j] - a[n - j]; in dfct()
525 xi = a[j] + a[n - j]; in dfct()
526 yr = a[k] - a[n - k]; in dfct()
527 yi = a[k] + a[n - k]; in dfct()
530 t[j] = xi - yi; in dfct()
533 t[mh] = a[mh] + a[n - mh]; in dfct()
534 a[mh] -= a[n - mh]; in dfct()
543 a[n - 1] = a[0] - a[1]; in dfct()
544 a[1] = a[0] + a[1]; in dfct()
545 for (j = m - 2; j >= 2; j -= 2) { in dfct()
546 a[2 * j + 1] = a[j] + a[j + 1]; in dfct()
547 a[2 * j - 1] = a[j] - a[j + 1]; in dfct()
560 a[n - l] = t[0] - t[1]; in dfct()
561 a[l] = t[0] + t[1]; in dfct()
565 a[k - l] = t[j] - t[j + 1]; in dfct()
566 a[k + l] = t[j] + t[j + 1]; in dfct()
568 l <<= 1; in dfct()
569 mh = m >> 1; in dfct()
571 k = m - j; in dfct()
572 t[j] = t[m + k] - t[m + j]; in dfct()
579 a[n] = t[2] - t[1]; in dfct()
580 a[0] = t[2] + t[1]; in dfct()
582 a[1] = a[0]; in dfct()
584 a[0] = t[1]; in dfct()
589 void dfst(int n, double *a, double *t, int *ip, double *w) in dfst() argument
594 if (n > FFT4G_MAX_SIZE) in dfst()
598 if (n > (nw << 3)) { in dfst()
599 nw = n >> 3; in dfst()
602 nc = ip[1]; in dfst()
603 if (n > (nc << 1)) { in dfst()
604 nc = n >> 1; in dfst()
607 if (n > 2) { in dfst()
608 m = n >> 1; in dfst()
609 mh = m >> 1; in dfst()
610 for (j = 1; j < mh; j++) { in dfst()
611 k = m - j; in dfst()
612 xr = a[j] + a[n - j]; in dfst()
613 xi = a[j] - a[n - j]; in dfst()
614 yr = a[k] + a[n - k]; in dfst()
615 yi = a[k] - a[n - k]; in dfst()
619 t[k] = xi - yi; in dfst()
621 t[0] = a[mh] - a[n - mh]; in dfst()
622 a[mh] += a[n - mh]; in dfst()
632 a[n - 1] = a[1] - a[0]; in dfst()
633 a[1] = a[0] + a[1]; in dfst()
634 for (j = m - 2; j >= 2; j -= 2) { in dfst()
635 a[2 * j + 1] = a[j] - a[j + 1]; in dfst()
636 a[2 * j - 1] = -a[j] - a[j + 1]; in dfst()
649 a[n - l] = t[1] - t[0]; in dfst()
650 a[l] = t[0] + t[1]; in dfst()
654 a[k - l] = -t[j] - t[j + 1]; in dfst()
655 a[k + l] = t[j] - t[j + 1]; in dfst()
657 l <<= 1; in dfst()
658 mh = m >> 1; in dfst()
659 for (j = 1; j < mh; j++) { in dfst()
660 k = m - j; in dfst()
662 t[k] = t[m + k] - t[m + j]; in dfst()
673 /* -------- initializing routines -------- */
682 ip[1] = 1; in makewt()
684 nwh = nw >> 1; in makewt()
686 w[0] = 1; in makewt()
687 w[1] = 0; in makewt()
689 w[nwh + 1] = w[nwh]; in makewt()
695 w[j + 1] = y; in makewt()
696 w[nw - j] = y; in makewt()
697 w[nw - j + 1] = x; in makewt()
710 ip[1] = nc; in makect()
711 if (nc > 1) { in makect()
712 nch = nc >> 1; in makect()
716 for (j = 1; j < nch; j++) { in makect()
718 c[nc - j] = 0.5 * sin(delta * j); in makect()
724 /* -------- child routines -------- */
727 static void bitrv2(int n, int *ip0, double *a) in bitrv2() argument
734 l = n; in bitrv2()
735 m = 1; in bitrv2()
737 l >>= 1; in bitrv2()
741 m <<= 1; in bitrv2()
750 xi = a[j1 + 1]; in bitrv2()
752 yi = a[k1 + 1]; in bitrv2()
754 a[j1 + 1] = yi; in bitrv2()
756 a[k1 + 1] = xi; in bitrv2()
760 xi = a[j1 + 1]; in bitrv2()
762 yi = a[k1 + 1]; in bitrv2()
764 a[j1 + 1] = yi; in bitrv2()
766 a[k1 + 1] = xi; in bitrv2()
768 k1 -= m2; in bitrv2()
770 xi = a[j1 + 1]; in bitrv2()
772 yi = a[k1 + 1]; in bitrv2()
774 a[j1 + 1] = yi; in bitrv2()
776 a[k1 + 1] = xi; in bitrv2()
780 xi = a[j1 + 1]; in bitrv2()
782 yi = a[k1 + 1]; in bitrv2()
784 a[j1 + 1] = yi; in bitrv2()
786 a[k1 + 1] = xi; in bitrv2()
791 xi = a[j1 + 1]; in bitrv2()
793 yi = a[k1 + 1]; in bitrv2()
795 a[j1 + 1] = yi; in bitrv2()
797 a[k1 + 1] = xi; in bitrv2()
800 for (k = 1; k < m; k++) { in bitrv2()
805 xi = a[j1 + 1]; in bitrv2()
807 yi = a[k1 + 1]; in bitrv2()
809 a[j1 + 1] = yi; in bitrv2()
811 a[k1 + 1] = xi; in bitrv2()
815 xi = a[j1 + 1]; in bitrv2()
817 yi = a[k1 + 1]; in bitrv2()
819 a[j1 + 1] = yi; in bitrv2()
821 a[k1 + 1] = xi; in bitrv2()
828 static void bitrv2conj(int n, int *ip0, double *a) in bitrv2conj() argument
835 l = n; in bitrv2conj()
836 m = 1; in bitrv2conj()
838 l >>= 1; in bitrv2conj()
842 m <<= 1; in bitrv2conj()
851 xi = -a[j1 + 1]; in bitrv2conj()
853 yi = -a[k1 + 1]; in bitrv2conj()
855 a[j1 + 1] = yi; in bitrv2conj()
857 a[k1 + 1] = xi; in bitrv2conj()
861 xi = -a[j1 + 1]; in bitrv2conj()
863 yi = -a[k1 + 1]; in bitrv2conj()
865 a[j1 + 1] = yi; in bitrv2conj()
867 a[k1 + 1] = xi; in bitrv2conj()
869 k1 -= m2; in bitrv2conj()
871 xi = -a[j1 + 1]; in bitrv2conj()
873 yi = -a[k1 + 1]; in bitrv2conj()
875 a[j1 + 1] = yi; in bitrv2conj()
877 a[k1 + 1] = xi; in bitrv2conj()
881 xi = -a[j1 + 1]; in bitrv2conj()
883 yi = -a[k1 + 1]; in bitrv2conj()
885 a[j1 + 1] = yi; in bitrv2conj()
887 a[k1 + 1] = xi; in bitrv2conj()
890 a[k1 + 1] = -a[k1 + 1]; in bitrv2conj()
894 xi = -a[j1 + 1]; in bitrv2conj()
896 yi = -a[k1 + 1]; in bitrv2conj()
898 a[j1 + 1] = yi; in bitrv2conj()
900 a[k1 + 1] = xi; in bitrv2conj()
902 a[k1 + 1] = -a[k1 + 1]; in bitrv2conj()
905 a[1] = -a[1]; in bitrv2conj()
906 a[m2 + 1] = -a[m2 + 1]; in bitrv2conj()
907 for (k = 1; k < m; k++) { in bitrv2conj()
912 xi = -a[j1 + 1]; in bitrv2conj()
914 yi = -a[k1 + 1]; in bitrv2conj()
916 a[j1 + 1] = yi; in bitrv2conj()
918 a[k1 + 1] = xi; in bitrv2conj()
922 xi = -a[j1 + 1]; in bitrv2conj()
924 yi = -a[k1 + 1]; in bitrv2conj()
926 a[j1 + 1] = yi; in bitrv2conj()
928 a[k1 + 1] = xi; in bitrv2conj()
931 a[k1 + 1] = -a[k1 + 1]; in bitrv2conj()
932 a[k1 + m2 + 1] = -a[k1 + m2 + 1]; in bitrv2conj()
938 static void cftfsub(int n, double *a, double const *w) in cftfsub() argument
944 if (n > 8) { in cftfsub()
945 cft1st(n, a, w); in cftfsub()
947 while ((l << 2) < n) { in cftfsub()
948 cftmdl(n, l, a, w); in cftfsub()
952 if ((l << 2) == n) { in cftfsub()
958 x0i = a[j + 1] + a[j1 + 1]; in cftfsub()
959 x1r = a[j] - a[j1]; in cftfsub()
960 x1i = a[j + 1] - a[j1 + 1]; in cftfsub()
962 x2i = a[j2 + 1] + a[j3 + 1]; in cftfsub()
963 x3r = a[j2] - a[j3]; in cftfsub()
964 x3i = a[j2 + 1] - a[j3 + 1]; in cftfsub()
966 a[j + 1] = x0i + x2i; in cftfsub()
967 a[j2] = x0r - x2r; in cftfsub()
968 a[j2 + 1] = x0i - x2i; in cftfsub()
969 a[j1] = x1r - x3i; in cftfsub()
970 a[j1 + 1] = x1i + x3r; in cftfsub()
972 a[j3 + 1] = x1i - x3r; in cftfsub()
977 x0r = a[j] - a[j1]; in cftfsub()
978 x0i = a[j + 1] - a[j1 + 1]; in cftfsub()
980 a[j + 1] += a[j1 + 1]; in cftfsub()
982 a[j1 + 1] = x0i; in cftfsub()
988 static void cftbsub(int n, double *a, double const *w) in cftbsub() argument
994 if (n > 8) { in cftbsub()
995 cft1st(n, a, w); in cftbsub()
997 while ((l << 2) < n) { in cftbsub()
998 cftmdl(n, l, a, w); in cftbsub()
1002 if ((l << 2) == n) { in cftbsub()
1008 x0i = -a[j + 1] - a[j1 + 1]; in cftbsub()
1009 x1r = a[j] - a[j1]; in cftbsub()
1010 x1i = -a[j + 1] + a[j1 + 1]; in cftbsub()
1012 x2i = a[j2 + 1] + a[j3 + 1]; in cftbsub()
1013 x3r = a[j2] - a[j3]; in cftbsub()
1014 x3i = a[j2 + 1] - a[j3 + 1]; in cftbsub()
1016 a[j + 1] = x0i - x2i; in cftbsub()
1017 a[j2] = x0r - x2r; in cftbsub()
1018 a[j2 + 1] = x0i + x2i; in cftbsub()
1019 a[j1] = x1r - x3i; in cftbsub()
1020 a[j1 + 1] = x1i - x3r; in cftbsub()
1022 a[j3 + 1] = x1i + x3r; in cftbsub()
1027 x0r = a[j] - a[j1]; in cftbsub()
1028 x0i = -a[j + 1] + a[j1 + 1]; in cftbsub()
1030 a[j + 1] = -a[j + 1] - a[j1 + 1]; in cftbsub()
1032 a[j1 + 1] = x0i; in cftbsub()
1038 static void cft1st(int n, double *a, double const *w) in cft1st() argument
1045 x0i = a[1] + a[3]; in cft1st()
1046 x1r = a[0] - a[2]; in cft1st()
1047 x1i = a[1] - a[3]; in cft1st()
1050 x3r = a[4] - a[6]; in cft1st()
1051 x3i = a[5] - a[7]; in cft1st()
1053 a[1] = x0i + x2i; in cft1st()
1054 a[4] = x0r - x2r; in cft1st()
1055 a[5] = x0i - x2i; in cft1st()
1056 a[2] = x1r - x3i; in cft1st()
1059 a[7] = x1i - x3r; in cft1st()
1063 x1r = a[8] - a[10]; in cft1st()
1064 x1i = a[9] - a[11]; in cft1st()
1067 x3r = a[12] - a[14]; in cft1st()
1068 x3i = a[13] - a[15]; in cft1st()
1071 a[12] = x2i - x0i; in cft1st()
1072 a[13] = x0r - x2r; in cft1st()
1073 x0r = x1r - x3i; in cft1st()
1075 a[10] = wk1r * (x0r - x0i); in cft1st()
1078 x0i = x3r - x1i; in cft1st()
1079 a[14] = wk1r * (x0i - x0r); in cft1st()
1082 for (j = 16; j < n; j += 16) { in cft1st()
1086 wk2i = w[k1 + 1]; in cft1st()
1088 wk1i = w[k2 + 1]; in cft1st()
1089 wk3r = wk1r - 2 * wk2i * wk1i; in cft1st()
1090 wk3i = 2 * wk2i * wk1r - wk1i; in cft1st()
1092 x0i = a[j + 1] + a[j + 3]; in cft1st()
1093 x1r = a[j] - a[j + 2]; in cft1st()
1094 x1i = a[j + 1] - a[j + 3]; in cft1st()
1097 x3r = a[j + 4] - a[j + 6]; in cft1st()
1098 x3i = a[j + 5] - a[j + 7]; in cft1st()
1100 a[j + 1] = x0i + x2i; in cft1st()
1101 x0r -= x2r; in cft1st()
1102 x0i -= x2i; in cft1st()
1103 a[j + 4] = wk2r * x0r - wk2i * x0i; in cft1st()
1105 x0r = x1r - x3i; in cft1st()
1107 a[j + 2] = wk1r * x0r - wk1i * x0i; in cft1st()
1110 x0i = x1i - x3r; in cft1st()
1111 a[j + 6] = wk3r * x0r - wk3i * x0i; in cft1st()
1115 wk3r = wk1r - 2 * wk2r * wk1i; in cft1st()
1116 wk3i = 2 * wk2r * wk1r - wk1i; in cft1st()
1119 x1r = a[j + 8] - a[j + 10]; in cft1st()
1120 x1i = a[j + 9] - a[j + 11]; in cft1st()
1123 x3r = a[j + 12] - a[j + 14]; in cft1st()
1124 x3i = a[j + 13] - a[j + 15]; in cft1st()
1127 x0r -= x2r; in cft1st()
1128 x0i -= x2i; in cft1st()
1129 a[j + 12] = -wk2i * x0r - wk2r * x0i; in cft1st()
1130 a[j + 13] = -wk2i * x0i + wk2r * x0r; in cft1st()
1131 x0r = x1r - x3i; in cft1st()
1133 a[j + 10] = wk1r * x0r - wk1i * x0i; in cft1st()
1136 x0i = x1i - x3r; in cft1st()
1137 a[j + 14] = wk3r * x0r - wk3i * x0i; in cft1st()
1143 static void cftmdl(int n, int l, double *a, double const *w) in cftmdl() argument
1155 x0i = a[j + 1] + a[j1 + 1]; in cftmdl()
1156 x1r = a[j] - a[j1]; in cftmdl()
1157 x1i = a[j + 1] - a[j1 + 1]; in cftmdl()
1159 x2i = a[j2 + 1] + a[j3 + 1]; in cftmdl()
1160 x3r = a[j2] - a[j3]; in cftmdl()
1161 x3i = a[j2 + 1] - a[j3 + 1]; in cftmdl()
1163 a[j + 1] = x0i + x2i; in cftmdl()
1164 a[j2] = x0r - x2r; in cftmdl()
1165 a[j2 + 1] = x0i - x2i; in cftmdl()
1166 a[j1] = x1r - x3i; in cftmdl()
1167 a[j1 + 1] = x1i + x3r; in cftmdl()
1169 a[j3 + 1] = x1i - x3r; in cftmdl()
1177 x0i = a[j + 1] + a[j1 + 1]; in cftmdl()
1178 x1r = a[j] - a[j1]; in cftmdl()
1179 x1i = a[j + 1] - a[j1 + 1]; in cftmdl()
1181 x2i = a[j2 + 1] + a[j3 + 1]; in cftmdl()
1182 x3r = a[j2] - a[j3]; in cftmdl()
1183 x3i = a[j2 + 1] - a[j3 + 1]; in cftmdl()
1185 a[j + 1] = x0i + x2i; in cftmdl()
1186 a[j2] = x2i - x0i; in cftmdl()
1187 a[j2 + 1] = x0r - x2r; in cftmdl()
1188 x0r = x1r - x3i; in cftmdl()
1190 a[j1] = wk1r * (x0r - x0i); in cftmdl()
1191 a[j1 + 1] = wk1r * (x0r + x0i); in cftmdl()
1193 x0i = x3r - x1i; in cftmdl()
1194 a[j3] = wk1r * (x0i - x0r); in cftmdl()
1195 a[j3 + 1] = wk1r * (x0i + x0r); in cftmdl()
1199 for (k = m2; k < n; k += m2) { in cftmdl()
1203 wk2i = w[k1 + 1]; in cftmdl()
1205 wk1i = w[k2 + 1]; in cftmdl()
1206 wk3r = wk1r - 2 * wk2i * wk1i; in cftmdl()
1207 wk3i = 2 * wk2i * wk1r - wk1i; in cftmdl()
1213 x0i = a[j + 1] + a[j1 + 1]; in cftmdl()
1214 x1r = a[j] - a[j1]; in cftmdl()
1215 x1i = a[j + 1] - a[j1 + 1]; in cftmdl()
1217 x2i = a[j2 + 1] + a[j3 + 1]; in cftmdl()
1218 x3r = a[j2] - a[j3]; in cftmdl()
1219 x3i = a[j2 + 1] - a[j3 + 1]; in cftmdl()
1221 a[j + 1] = x0i + x2i; in cftmdl()
1222 x0r -= x2r; in cftmdl()
1223 x0i -= x2i; in cftmdl()
1224 a[j2] = wk2r * x0r - wk2i * x0i; in cftmdl()
1225 a[j2 + 1] = wk2r * x0i + wk2i * x0r; in cftmdl()
1226 x0r = x1r - x3i; in cftmdl()
1228 a[j1] = wk1r * x0r - wk1i * x0i; in cftmdl()
1229 a[j1 + 1] = wk1r * x0i + wk1i * x0r; in cftmdl()
1231 x0i = x1i - x3r; in cftmdl()
1232 a[j3] = wk3r * x0r - wk3i * x0i; in cftmdl()
1233 a[j3 + 1] = wk3r * x0i + wk3i * x0r; in cftmdl()
1237 wk3r = wk1r - 2 * wk2r * wk1i; in cftmdl()
1238 wk3i = 2 * wk2r * wk1r - wk1i; in cftmdl()
1244 x0i = a[j + 1] + a[j1 + 1]; in cftmdl()
1245 x1r = a[j] - a[j1]; in cftmdl()
1246 x1i = a[j + 1] - a[j1 + 1]; in cftmdl()
1248 x2i = a[j2 + 1] + a[j3 + 1]; in cftmdl()
1249 x3r = a[j2] - a[j3]; in cftmdl()
1250 x3i = a[j2 + 1] - a[j3 + 1]; in cftmdl()
1252 a[j + 1] = x0i + x2i; in cftmdl()
1253 x0r -= x2r; in cftmdl()
1254 x0i -= x2i; in cftmdl()
1255 a[j2] = -wk2i * x0r - wk2r * x0i; in cftmdl()
1256 a[j2 + 1] = -wk2i * x0i + wk2r * x0r; in cftmdl()
1257 x0r = x1r - x3i; in cftmdl()
1259 a[j1] = wk1r * x0r - wk1i * x0i; in cftmdl()
1260 a[j1 + 1] = wk1r * x0i + wk1i * x0r; in cftmdl()
1262 x0i = x1i - x3r; in cftmdl()
1263 a[j3] = wk3r * x0r - wk3i * x0i; in cftmdl()
1264 a[j3 + 1] = wk3r * x0i + wk3i * x0r; in cftmdl()
1270 static void rftfsub(int n, double *a, int nc, double const *c) in rftfsub() argument
1275 m = n >> 1; in rftfsub()
1279 k = n - j; in rftfsub()
1281 wkr = 0.5 - c[nc - kk]; in rftfsub()
1283 xr = a[j] - a[k]; in rftfsub()
1284 xi = a[j + 1] + a[k + 1]; in rftfsub()
1285 yr = wkr * xr - wki * xi; in rftfsub()
1287 a[j] -= yr; in rftfsub()
1288 a[j + 1] -= yi; in rftfsub()
1290 a[k + 1] -= yi; in rftfsub()
1295 static void rftbsub(int n, double *a, int nc, double const *c) in rftbsub() argument
1300 a[1] = -a[1]; in rftbsub()
1301 m = n >> 1; in rftbsub()
1305 k = n - j; in rftbsub()
1307 wkr = 0.5 - c[nc - kk]; in rftbsub()
1309 xr = a[j] - a[k]; in rftbsub()
1310 xi = a[j + 1] + a[k + 1]; in rftbsub()
1312 yi = wkr * xi - wki * xr; in rftbsub()
1313 a[j] -= yr; in rftbsub()
1314 a[j + 1] = yi - a[j + 1]; in rftbsub()
1316 a[k + 1] = yi - a[k + 1]; in rftbsub()
1318 a[m + 1] = -a[m + 1]; in rftbsub()
1322 static void dctsub(int n, double *a, int nc, double const *c) in dctsub() argument
1327 m = n >> 1; in dctsub()
1328 ks = nc / n; in dctsub()
1330 for (j = 1; j < m; j++) { in dctsub()
1331 k = n - j; in dctsub()
1333 wkr = c[kk] - c[nc - kk]; in dctsub()
1334 wki = c[kk] + c[nc - kk]; in dctsub()
1335 xr = wki * a[j] - wkr * a[k]; in dctsub()
1343 static void dstsub(int n, double *a, int nc, double const *c) in dstsub() argument
1348 m = n >> 1; in dstsub()
1349 ks = nc / n; in dstsub()
1351 for (j = 1; j < m; j++) { in dstsub()
1352 k = n - j; in dstsub()
1354 wkr = c[kk] - c[nc - kk]; in dstsub()
1355 wki = c[kk] + c[nc - kk]; in dstsub()
1356 xr = wki * a[k] - wkr * a[j]; in dstsub()