1 /* EAW - May 11th, 2001: Changed all doubles to floats */
2
3 /*
4 Copyright:
5 Copyright(C) 1996-1999 Takuya OOURA
6 email: ooura@mmm.t.u-tokyo.ac.jp
7 download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
8 You may use, copy, modify this code for any purpose and
9 without fee. You may distribute this ORIGINAL package.
10 */
11
12 /*
13 Fast Fourier/Cosine/Sine Transform
14 dimension :one
15 data length :power of 2
16 decimation :frequency
17 radix :4, 2
18 data :inplace
19 table :use
20 functions
21 cdft: Complex Discrete Fourier Transform
22 rdft: Real Discrete Fourier Transform
23 ddct: Discrete Cosine Transform
24 ddst: Discrete Sine Transform
25 dfct: Cosine Transform of RDFT (Real Symmetric DFT)
26 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
27 function prototypes
28 void cdft(int, int, float *, int *, float *);
29 void rdft(int, int, float *, int *, float *);
30 void ddct(int, int, float *, int *, float *);
31 void ddst(int, int, float *, int *, float *);
32 void dfct(int, float *, float *, int *, float *);
33 void dfst(int, float *, float *, int *, float *);
34
35
36 -------- Complex DFT (Discrete Fourier Transform) --------
37 [definition]
38 <case1>
39 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
40 <case2>
41 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
42 (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
43 [usage]
44 <case1>
45 ip[0] = 0; // first time only
46 cdft(2*n, 1, a, ip, w);
47 <case2>
48 ip[0] = 0; // first time only
49 cdft(2*n, -1, a, ip, w);
50 [parameters]
51 2*n :data length (int)
52 n >= 1, n = power of 2
53 a[0...2*n-1] :input/output data (float *)
54 input data
55 a[2*j] = Re(x[j]),
56 a[2*j+1] = Im(x[j]), 0<=j<n
57 output data
58 a[2*k] = Re(X[k]),
59 a[2*k+1] = Im(X[k]), 0<=k<n
60 ip[0...*] :work area for bit reversal (int *)
61 length of ip >= 2+sqrt(n)
62 strictly,
63 length of ip >=
64 2+(1<<(int)(log(n+0.5)/log(2))/2).
65 ip[0],ip[1] are pointers of the cos/sin table.
66 w[0...n/2-1] :cos/sin table (float *)
67 w[],ip[] are initialized if ip[0] == 0.
68 [remark]
69 Inverse of
70 cdft(2*n, -1, a, ip, w);
71 is
72 cdft(2*n, 1, a, ip, w);
73 for (j = 0; j <= 2 * n - 1; j++) {
74 a[j] *= 1.0 / n;
75 }
76 .
77
78
79 -------- Real DFT / Inverse of Real DFT --------
80 [definition]
81 <case1> RDFT
82 R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
83 I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
84 <case2> IRDFT (excluding scale)
85 a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
86 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
87 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
88 [usage]
89 <case1>
90 ip[0] = 0; // first time only
91 rdft(n, 1, a, ip, w);
92 <case2>
93 ip[0] = 0; // first time only
94 rdft(n, -1, a, ip, w);
95 [parameters]
96 n :data length (int)
97 n >= 2, n = power of 2
98 a[0...n-1] :input/output data (float *)
99 <case1>
100 output data
101 a[2*k] = R[k], 0<=k<n/2
102 a[2*k+1] = I[k], 0<k<n/2
103 a[1] = R[n/2]
104 <case2>
105 input data
106 a[2*j] = R[j], 0<=j<n/2
107 a[2*j+1] = I[j], 0<j<n/2
108 a[1] = R[n/2]
109 ip[0...*] :work area for bit reversal (int *)
110 length of ip >= 2+sqrt(n/2)
111 strictly,
112 length of ip >=
113 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
114 ip[0],ip[1] are pointers of the cos/sin table.
115 w[0...n/2-1] :cos/sin table (float *)
116 w[],ip[] are initialized if ip[0] == 0.
117 [remark]
118 Inverse of
119 rdft(n, 1, a, ip, w);
120 is
121 rdft(n, -1, a, ip, w);
122 for (j = 0; j <= n - 1; j++) {
123 a[j] *= 2.0 / n;
124 }
125 .
126
127
128 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
129 [definition]
130 <case1> IDCT (excluding scale)
131 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
132 <case2> DCT
133 C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
134 [usage]
135 <case1>
136 ip[0] = 0; // first time only
137 ddct(n, 1, a, ip, w);
138 <case2>
139 ip[0] = 0; // first time only
140 ddct(n, -1, a, ip, w);
141 [parameters]
142 n :data length (int)
143 n >= 2, n = power of 2
144 a[0...n-1] :input/output data (float *)
145 output data
146 a[k] = C[k], 0<=k<n
147 ip[0...*] :work area for bit reversal (int *)
148 length of ip >= 2+sqrt(n/2)
149 strictly,
150 length of ip >=
151 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
152 ip[0],ip[1] are pointers of the cos/sin table.
153 w[0...n*5/4-1] :cos/sin table (float *)
154 w[],ip[] are initialized if ip[0] == 0.
155 [remark]
156 Inverse of
157 ddct(n, -1, a, ip, w);
158 is
159 a[0] *= 0.5;
160 ddct(n, 1, a, ip, w);
161 for (j = 0; j <= n - 1; j++) {
162 a[j] *= 2.0 / n;
163 }
164 .
165
166
167 -------- DST (Discrete Sine Transform) / Inverse of DST --------
168 [definition]
169 <case1> IDST (excluding scale)
170 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
171 <case2> DST
172 S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
173 [usage]
174 <case1>
175 ip[0] = 0; // first time only
176 ddst(n, 1, a, ip, w);
177 <case2>
178 ip[0] = 0; // first time only
179 ddst(n, -1, a, ip, w);
180 [parameters]
181 n :data length (int)
182 n >= 2, n = power of 2
183 a[0...n-1] :input/output data (float *)
184 <case1>
185 input data
186 a[j] = A[j], 0<j<n
187 a[0] = A[n]
188 output data
189 a[k] = S[k], 0<=k<n
190 <case2>
191 output data
192 a[k] = S[k], 0<k<n
193 a[0] = S[n]
194 ip[0...*] :work area for bit reversal (int *)
195 length of ip >= 2+sqrt(n/2)
196 strictly,
197 length of ip >=
198 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
199 ip[0],ip[1] are pointers of the cos/sin table.
200 w[0...n*5/4-1] :cos/sin table (float *)
201 w[],ip[] are initialized if ip[0] == 0.
202 [remark]
203 Inverse of
204 ddst(n, -1, a, ip, w);
205 is
206 a[0] *= 0.5;
207 ddst(n, 1, a, ip, w);
208 for (j = 0; j <= n - 1; j++) {
209 a[j] *= 2.0 / n;
210 }
211 .
212
213
214 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
215 [definition]
216 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
217 [usage]
218 ip[0] = 0; // first time only
219 dfct(n, a, t, ip, w);
220 [parameters]
221 n :data length - 1 (int)
222 n >= 2, n = power of 2
223 a[0...n] :input/output data (float *)
224 output data
225 a[k] = C[k], 0<=k<=n
226 t[0...n/2] :work area (float *)
227 ip[0...*] :work area for bit reversal (int *)
228 length of ip >= 2+sqrt(n/4)
229 strictly,
230 length of ip >=
231 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
232 ip[0],ip[1] are pointers of the cos/sin table.
233 w[0...n*5/8-1] :cos/sin table (float *)
234 w[],ip[] are initialized if ip[0] == 0.
235 [remark]
236 Inverse of
237 a[0] *= 0.5;
238 a[n] *= 0.5;
239 dfct(n, a, t, ip, w);
240 is
241 a[0] *= 0.5;
242 a[n] *= 0.5;
243 dfct(n, a, t, ip, w);
244 for (j = 0; j <= n; j++) {
245 a[j] *= 2.0 / n;
246 }
247 .
248
249
250 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
251 [definition]
252 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
253 [usage]
254 ip[0] = 0; // first time only
255 dfst(n, a, t, ip, w);
256 [parameters]
257 n :data length + 1 (int)
258 n >= 2, n = power of 2
259 a[0...n-1] :input/output data (float *)
260 output data
261 a[k] = S[k], 0<k<n
262 (a[0] is used for work area)
263 t[0...n/2-1] :work area (float *)
264 ip[0...*] :work area for bit reversal (int *)
265 length of ip >= 2+sqrt(n/4)
266 strictly,
267 length of ip >=
268 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
269 ip[0],ip[1] are pointers of the cos/sin table.
270 w[0...n*5/8-1] :cos/sin table (float *)
271 w[],ip[] are initialized if ip[0] == 0.
272 [remark]
273 Inverse of
274 dfst(n, a, t, ip, w);
275 is
276 dfst(n, a, t, ip, w);
277 for (j = 1; j <= n - 1; j++) {
278 a[j] *= 2.0 / n;
279 }
280 .
281
282
283 Appendix :
284 The cos/sin table is recalculated when the larger table required.
285 w[] and ip[] are compatible with all routines.
286 */
287
288 #include <math.h>
289
290
291 namespace TimidityPlus
292 {
293 void makewt(int nw, int *ip, float *w);
294 void bitrv2(int n, int *ip, float *a);
295 void bitrv2conj(int n, int *ip, float *a);
296 void cftfsub(int n, float *a, float *w);
297 void cftbsub(int n, float *a, float *w);
298 void makect(int nc, int *ip, float *c);
299 void rftfsub(int n, float *a, int nc, float *c);
300 void rftbsub(int n, float *a, int nc, float *c);
301 void dctsub(int n, float *a, int nc, float *c);
302 void dstsub(int n, float *a, int nc, float *c);
303 void cft1st(int n, float *a, float *w);
304 void cftmdl(int n, int l, float *a, float *w);
305
306
307
308 #ifdef _MSC_VER
309 #pragma warning(disable:4244) // conversion from 'double' to 'float', possible loss of data
310 #endif
311
312
cdft(int n,int isgn,float * a,int * ip,float * w)313 void cdft(int n, int isgn, float *a, int *ip, float *w)
314 {
315 if (n > (ip[0] << 2)) {
316 makewt(n >> 2, ip, w);
317 }
318 if (n > 4) {
319 if (isgn >= 0) {
320 bitrv2(n, ip + 2, a);
321 cftfsub(n, a, w);
322 }
323 else {
324 bitrv2conj(n, ip + 2, a);
325 cftbsub(n, a, w);
326 }
327 }
328 else if (n == 4) {
329 cftfsub(n, a, w);
330 }
331 }
332
333
rdft(int n,int isgn,float * a,int * ip,float * w)334 void rdft(int n, int isgn, float *a, int *ip, float *w)
335 {
336 int nw, nc;
337 float xi;
338
339 nw = ip[0];
340 if (n > (nw << 2)) {
341 nw = n >> 2;
342 makewt(nw, ip, w);
343 }
344 nc = ip[1];
345 if (n > (nc << 2)) {
346 nc = n >> 2;
347 makect(nc, ip, w + nw);
348 }
349 if (isgn >= 0) {
350 if (n > 4) {
351 bitrv2(n, ip + 2, a);
352 cftfsub(n, a, w);
353 rftfsub(n, a, nc, w + nw);
354 }
355 else if (n == 4) {
356 cftfsub(n, a, w);
357 }
358 xi = a[0] - a[1];
359 a[0] += a[1];
360 a[1] = xi;
361 }
362 else {
363 a[1] = 0.5 * (a[0] - a[1]);
364 a[0] -= a[1];
365 if (n > 4) {
366 rftbsub(n, a, nc, w + nw);
367 bitrv2(n, ip + 2, a);
368 cftbsub(n, a, w);
369 }
370 else if (n == 4) {
371 cftfsub(n, a, w);
372 }
373 }
374 }
375
376
ddct(int n,int isgn,float * a,int * ip,float * w)377 void ddct(int n, int isgn, float *a, int *ip, float *w)
378 {
379 int j, nw, nc;
380 float xr;
381
382 nw = ip[0];
383 if (n > (nw << 2)) {
384 nw = n >> 2;
385 makewt(nw, ip, w);
386 }
387 nc = ip[1];
388 if (n > nc) {
389 nc = n;
390 makect(nc, ip, w + nw);
391 }
392 if (isgn < 0) {
393 xr = a[n - 1];
394 for (j = n - 2; j >= 2; j -= 2) {
395 a[j + 1] = a[j] - a[j - 1];
396 a[j] += a[j - 1];
397 }
398 a[1] = a[0] - xr;
399 a[0] += xr;
400 if (n > 4) {
401 rftbsub(n, a, nc, w + nw);
402 bitrv2(n, ip + 2, a);
403 cftbsub(n, a, w);
404 }
405 else if (n == 4) {
406 cftfsub(n, a, w);
407 }
408 }
409 dctsub(n, a, nc, w + nw);
410 if (isgn >= 0) {
411 if (n > 4) {
412 bitrv2(n, ip + 2, a);
413 cftfsub(n, a, w);
414 rftfsub(n, a, nc, w + nw);
415 }
416 else if (n == 4) {
417 cftfsub(n, a, w);
418 }
419 xr = a[0] - a[1];
420 a[0] += a[1];
421 for (j = 2; j < n; j += 2) {
422 a[j - 1] = a[j] - a[j + 1];
423 a[j] += a[j + 1];
424 }
425 a[n - 1] = xr;
426 }
427 }
428
429
ddst(int n,int isgn,float * a,int * ip,float * w)430 void ddst(int n, int isgn, float *a, int *ip, float *w)
431 {
432 int j, nw, nc;
433 float xr;
434
435 nw = ip[0];
436 if (n > (nw << 2)) {
437 nw = n >> 2;
438 makewt(nw, ip, w);
439 }
440 nc = ip[1];
441 if (n > nc) {
442 nc = n;
443 makect(nc, ip, w + nw);
444 }
445 if (isgn < 0) {
446 xr = a[n - 1];
447 for (j = n - 2; j >= 2; j -= 2) {
448 a[j + 1] = -a[j] - a[j - 1];
449 a[j] -= a[j - 1];
450 }
451 a[1] = a[0] + xr;
452 a[0] -= xr;
453 if (n > 4) {
454 rftbsub(n, a, nc, w + nw);
455 bitrv2(n, ip + 2, a);
456 cftbsub(n, a, w);
457 }
458 else if (n == 4) {
459 cftfsub(n, a, w);
460 }
461 }
462 dstsub(n, a, nc, w + nw);
463 if (isgn >= 0) {
464 if (n > 4) {
465 bitrv2(n, ip + 2, a);
466 cftfsub(n, a, w);
467 rftfsub(n, a, nc, w + nw);
468 }
469 else if (n == 4) {
470 cftfsub(n, a, w);
471 }
472 xr = a[0] - a[1];
473 a[0] += a[1];
474 for (j = 2; j < n; j += 2) {
475 a[j - 1] = -a[j] - a[j + 1];
476 a[j] -= a[j + 1];
477 }
478 a[n - 1] = -xr;
479 }
480 }
481
482
dfct(int n,float * a,float * t,int * ip,float * w)483 void dfct(int n, float *a, float *t, int *ip, float *w)
484 {
485 int j, k, l, m, mh, nw, nc;
486 float xr, xi, yr, yi;
487
488 nw = ip[0];
489 if (n > (nw << 3)) {
490 nw = n >> 3;
491 makewt(nw, ip, w);
492 }
493 nc = ip[1];
494 if (n > (nc << 1)) {
495 nc = n >> 1;
496 makect(nc, ip, w + nw);
497 }
498 m = n >> 1;
499 yi = a[m];
500 xi = a[0] + a[n];
501 a[0] -= a[n];
502 t[0] = xi - yi;
503 t[m] = xi + yi;
504 if (n > 2) {
505 mh = m >> 1;
506 for (j = 1; j < mh; j++) {
507 k = m - j;
508 xr = a[j] - a[n - j];
509 xi = a[j] + a[n - j];
510 yr = a[k] - a[n - k];
511 yi = a[k] + a[n - k];
512 a[j] = xr;
513 a[k] = yr;
514 t[j] = xi - yi;
515 t[k] = xi + yi;
516 }
517 t[mh] = a[mh] + a[n - mh];
518 a[mh] -= a[n - mh];
519 dctsub(m, a, nc, w + nw);
520 if (m > 4) {
521 bitrv2(m, ip + 2, a);
522 cftfsub(m, a, w);
523 rftfsub(m, a, nc, w + nw);
524 }
525 else if (m == 4) {
526 cftfsub(m, a, w);
527 }
528 a[n - 1] = a[0] - a[1];
529 a[1] = a[0] + a[1];
530 for (j = m - 2; j >= 2; j -= 2) {
531 a[2 * j + 1] = a[j] + a[j + 1];
532 a[2 * j - 1] = a[j] - a[j + 1];
533 }
534 l = 2;
535 m = mh;
536 while (m >= 2) {
537 dctsub(m, t, nc, w + nw);
538 if (m > 4) {
539 bitrv2(m, ip + 2, t);
540 cftfsub(m, t, w);
541 rftfsub(m, t, nc, w + nw);
542 }
543 else if (m == 4) {
544 cftfsub(m, t, w);
545 }
546 a[n - l] = t[0] - t[1];
547 a[l] = t[0] + t[1];
548 k = 0;
549 for (j = 2; j < m; j += 2) {
550 k += l << 2;
551 a[k - l] = t[j] - t[j + 1];
552 a[k + l] = t[j] + t[j + 1];
553 }
554 l <<= 1;
555 mh = m >> 1;
556 for (j = 0; j < mh; j++) {
557 k = m - j;
558 t[j] = t[m + k] - t[m + j];
559 t[k] = t[m + k] + t[m + j];
560 }
561 t[mh] = t[m + mh];
562 m = mh;
563 }
564 a[l] = t[0];
565 a[n] = t[2] - t[1];
566 a[0] = t[2] + t[1];
567 }
568 else {
569 a[1] = a[0];
570 a[2] = t[0];
571 a[0] = t[1];
572 }
573 }
574
575
dfst(int n,float * a,float * t,int * ip,float * w)576 void dfst(int n, float *a, float *t, int *ip, float *w)
577 {
578 int j, k, l, m, mh, nw, nc;
579 float xr, xi, yr, yi;
580
581 nw = ip[0];
582 if (n > (nw << 3)) {
583 nw = n >> 3;
584 makewt(nw, ip, w);
585 }
586 nc = ip[1];
587 if (n > (nc << 1)) {
588 nc = n >> 1;
589 makect(nc, ip, w + nw);
590 }
591 if (n > 2) {
592 m = n >> 1;
593 mh = m >> 1;
594 for (j = 1; j < mh; j++) {
595 k = m - j;
596 xr = a[j] + a[n - j];
597 xi = a[j] - a[n - j];
598 yr = a[k] + a[n - k];
599 yi = a[k] - a[n - k];
600 a[j] = xr;
601 a[k] = yr;
602 t[j] = xi + yi;
603 t[k] = xi - yi;
604 }
605 t[0] = a[mh] - a[n - mh];
606 a[mh] += a[n - mh];
607 a[0] = a[m];
608 dstsub(m, a, nc, w + nw);
609 if (m > 4) {
610 bitrv2(m, ip + 2, a);
611 cftfsub(m, a, w);
612 rftfsub(m, a, nc, w + nw);
613 }
614 else if (m == 4) {
615 cftfsub(m, a, w);
616 }
617 a[n - 1] = a[1] - a[0];
618 a[1] = a[0] + a[1];
619 for (j = m - 2; j >= 2; j -= 2) {
620 a[2 * j + 1] = a[j] - a[j + 1];
621 a[2 * j - 1] = -a[j] - a[j + 1];
622 }
623 l = 2;
624 m = mh;
625 while (m >= 2) {
626 dstsub(m, t, nc, w + nw);
627 if (m > 4) {
628 bitrv2(m, ip + 2, t);
629 cftfsub(m, t, w);
630 rftfsub(m, t, nc, w + nw);
631 }
632 else if (m == 4) {
633 cftfsub(m, t, w);
634 }
635 a[n - l] = t[1] - t[0];
636 a[l] = t[0] + t[1];
637 k = 0;
638 for (j = 2; j < m; j += 2) {
639 k += l << 2;
640 a[k - l] = -t[j] - t[j + 1];
641 a[k + l] = t[j] - t[j + 1];
642 }
643 l <<= 1;
644 mh = m >> 1;
645 for (j = 1; j < mh; j++) {
646 k = m - j;
647 t[j] = t[m + k] + t[m + j];
648 t[k] = t[m + k] - t[m + j];
649 }
650 t[0] = t[m + mh];
651 m = mh;
652 }
653 a[l] = t[0];
654 }
655 a[0] = 0;
656 }
657
658
659 /* -------- initializing routines -------- */
660
661
662
makewt(int nw,int * ip,float * w)663 void makewt(int nw, int *ip, float *w)
664 {
665 void bitrv2(int n, int *ip, float *a);
666 int j, nwh;
667 float delta, x, y;
668
669 ip[0] = nw;
670 ip[1] = 1;
671 if (nw > 2) {
672 nwh = nw >> 1;
673 delta = atan(1.0) / nwh;
674 w[0] = 1;
675 w[1] = 0;
676 w[nwh] = cos(delta * nwh);
677 w[nwh + 1] = w[nwh];
678 if (nwh > 2) {
679 for (j = 2; j < nwh; j += 2) {
680 x = cos(delta * j);
681 y = sin(delta * j);
682 w[j] = x;
683 w[j + 1] = y;
684 w[nw - j] = y;
685 w[nw - j + 1] = x;
686 }
687 bitrv2(nw, ip + 2, w);
688 }
689 }
690 }
691
692
makect(int nc,int * ip,float * c)693 void makect(int nc, int *ip, float *c)
694 {
695 int j, nch;
696 float delta;
697
698 ip[1] = nc;
699 if (nc > 1) {
700 nch = nc >> 1;
701 delta = atan(1.0) / nch;
702 c[0] = cos(delta * nch);
703 c[nch] = 0.5 * c[0];
704 for (j = 1; j < nch; j++) {
705 c[j] = 0.5 * cos(delta * j);
706 c[nc - j] = 0.5 * sin(delta * j);
707 }
708 }
709 }
710
711
712 /* -------- child routines -------- */
713
714
bitrv2(int n,int * ip,float * a)715 void bitrv2(int n, int *ip, float *a)
716 {
717 int j, j1, k, k1, l, m, m2;
718 float xr, xi, yr, yi;
719
720 ip[0] = 0;
721 l = n;
722 m = 1;
723 while ((m << 3) < l) {
724 l >>= 1;
725 for (j = 0; j < m; j++) {
726 ip[m + j] = ip[j] + l;
727 }
728 m <<= 1;
729 }
730 m2 = 2 * m;
731 if ((m << 3) == l) {
732 for (k = 0; k < m; k++) {
733 for (j = 0; j < k; j++) {
734 j1 = 2 * j + ip[k];
735 k1 = 2 * k + ip[j];
736 xr = a[j1];
737 xi = a[j1 + 1];
738 yr = a[k1];
739 yi = a[k1 + 1];
740 a[j1] = yr;
741 a[j1 + 1] = yi;
742 a[k1] = xr;
743 a[k1 + 1] = xi;
744 j1 += m2;
745 k1 += 2 * m2;
746 xr = a[j1];
747 xi = a[j1 + 1];
748 yr = a[k1];
749 yi = a[k1 + 1];
750 a[j1] = yr;
751 a[j1 + 1] = yi;
752 a[k1] = xr;
753 a[k1 + 1] = xi;
754 j1 += m2;
755 k1 -= m2;
756 xr = a[j1];
757 xi = a[j1 + 1];
758 yr = a[k1];
759 yi = a[k1 + 1];
760 a[j1] = yr;
761 a[j1 + 1] = yi;
762 a[k1] = xr;
763 a[k1 + 1] = xi;
764 j1 += m2;
765 k1 += 2 * m2;
766 xr = a[j1];
767 xi = a[j1 + 1];
768 yr = a[k1];
769 yi = a[k1 + 1];
770 a[j1] = yr;
771 a[j1 + 1] = yi;
772 a[k1] = xr;
773 a[k1 + 1] = xi;
774 }
775 j1 = 2 * k + m2 + ip[k];
776 k1 = j1 + m2;
777 xr = a[j1];
778 xi = a[j1 + 1];
779 yr = a[k1];
780 yi = a[k1 + 1];
781 a[j1] = yr;
782 a[j1 + 1] = yi;
783 a[k1] = xr;
784 a[k1 + 1] = xi;
785 }
786 }
787 else {
788 for (k = 1; k < m; k++) {
789 for (j = 0; j < k; j++) {
790 j1 = 2 * j + ip[k];
791 k1 = 2 * k + ip[j];
792 xr = a[j1];
793 xi = a[j1 + 1];
794 yr = a[k1];
795 yi = a[k1 + 1];
796 a[j1] = yr;
797 a[j1 + 1] = yi;
798 a[k1] = xr;
799 a[k1 + 1] = xi;
800 j1 += m2;
801 k1 += m2;
802 xr = a[j1];
803 xi = a[j1 + 1];
804 yr = a[k1];
805 yi = a[k1 + 1];
806 a[j1] = yr;
807 a[j1 + 1] = yi;
808 a[k1] = xr;
809 a[k1 + 1] = xi;
810 }
811 }
812 }
813 }
814
815
bitrv2conj(int n,int * ip,float * a)816 void bitrv2conj(int n, int *ip, float *a)
817 {
818 int j, j1, k, k1, l, m, m2;
819 float xr, xi, yr, yi;
820
821 ip[0] = 0;
822 l = n;
823 m = 1;
824 while ((m << 3) < l) {
825 l >>= 1;
826 for (j = 0; j < m; j++) {
827 ip[m + j] = ip[j] + l;
828 }
829 m <<= 1;
830 }
831 m2 = 2 * m;
832 if ((m << 3) == l) {
833 for (k = 0; k < m; k++) {
834 for (j = 0; j < k; j++) {
835 j1 = 2 * j + ip[k];
836 k1 = 2 * k + ip[j];
837 xr = a[j1];
838 xi = -a[j1 + 1];
839 yr = a[k1];
840 yi = -a[k1 + 1];
841 a[j1] = yr;
842 a[j1 + 1] = yi;
843 a[k1] = xr;
844 a[k1 + 1] = xi;
845 j1 += m2;
846 k1 += 2 * m2;
847 xr = a[j1];
848 xi = -a[j1 + 1];
849 yr = a[k1];
850 yi = -a[k1 + 1];
851 a[j1] = yr;
852 a[j1 + 1] = yi;
853 a[k1] = xr;
854 a[k1 + 1] = xi;
855 j1 += m2;
856 k1 -= m2;
857 xr = a[j1];
858 xi = -a[j1 + 1];
859 yr = a[k1];
860 yi = -a[k1 + 1];
861 a[j1] = yr;
862 a[j1 + 1] = yi;
863 a[k1] = xr;
864 a[k1 + 1] = xi;
865 j1 += m2;
866 k1 += 2 * m2;
867 xr = a[j1];
868 xi = -a[j1 + 1];
869 yr = a[k1];
870 yi = -a[k1 + 1];
871 a[j1] = yr;
872 a[j1 + 1] = yi;
873 a[k1] = xr;
874 a[k1 + 1] = xi;
875 }
876 k1 = 2 * k + ip[k];
877 a[k1 + 1] = -a[k1 + 1];
878 j1 = k1 + m2;
879 k1 = j1 + m2;
880 xr = a[j1];
881 xi = -a[j1 + 1];
882 yr = a[k1];
883 yi = -a[k1 + 1];
884 a[j1] = yr;
885 a[j1 + 1] = yi;
886 a[k1] = xr;
887 a[k1 + 1] = xi;
888 k1 += m2;
889 a[k1 + 1] = -a[k1 + 1];
890 }
891 }
892 else {
893 a[1] = -a[1];
894 a[m2 + 1] = -a[m2 + 1];
895 for (k = 1; k < m; k++) {
896 for (j = 0; j < k; j++) {
897 j1 = 2 * j + ip[k];
898 k1 = 2 * k + ip[j];
899 xr = a[j1];
900 xi = -a[j1 + 1];
901 yr = a[k1];
902 yi = -a[k1 + 1];
903 a[j1] = yr;
904 a[j1 + 1] = yi;
905 a[k1] = xr;
906 a[k1 + 1] = xi;
907 j1 += m2;
908 k1 += m2;
909 xr = a[j1];
910 xi = -a[j1 + 1];
911 yr = a[k1];
912 yi = -a[k1 + 1];
913 a[j1] = yr;
914 a[j1 + 1] = yi;
915 a[k1] = xr;
916 a[k1 + 1] = xi;
917 }
918 k1 = 2 * k + ip[k];
919 a[k1 + 1] = -a[k1 + 1];
920 a[k1 + m2 + 1] = -a[k1 + m2 + 1];
921 }
922 }
923 }
924
925
cftfsub(int n,float * a,float * w)926 void cftfsub(int n, float *a, float *w)
927 {
928 int j, j1, j2, j3, l;
929 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
930
931 l = 2;
932 if (n > 8) {
933 cft1st(n, a, w);
934 l = 8;
935 while ((l << 2) < n) {
936 cftmdl(n, l, a, w);
937 l <<= 2;
938 }
939 }
940 if ((l << 2) == n) {
941 for (j = 0; j < l; j += 2) {
942 j1 = j + l;
943 j2 = j1 + l;
944 j3 = j2 + l;
945 x0r = a[j] + a[j1];
946 x0i = a[j + 1] + a[j1 + 1];
947 x1r = a[j] - a[j1];
948 x1i = a[j + 1] - a[j1 + 1];
949 x2r = a[j2] + a[j3];
950 x2i = a[j2 + 1] + a[j3 + 1];
951 x3r = a[j2] - a[j3];
952 x3i = a[j2 + 1] - a[j3 + 1];
953 a[j] = x0r + x2r;
954 a[j + 1] = x0i + x2i;
955 a[j2] = x0r - x2r;
956 a[j2 + 1] = x0i - x2i;
957 a[j1] = x1r - x3i;
958 a[j1 + 1] = x1i + x3r;
959 a[j3] = x1r + x3i;
960 a[j3 + 1] = x1i - x3r;
961 }
962 }
963 else {
964 for (j = 0; j < l; j += 2) {
965 j1 = j + l;
966 x0r = a[j] - a[j1];
967 x0i = a[j + 1] - a[j1 + 1];
968 a[j] += a[j1];
969 a[j + 1] += a[j1 + 1];
970 a[j1] = x0r;
971 a[j1 + 1] = x0i;
972 }
973 }
974 }
975
976
cftbsub(int n,float * a,float * w)977 void cftbsub(int n, float *a, float *w)
978 {
979 int j, j1, j2, j3, l;
980 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
981
982 l = 2;
983 if (n > 8) {
984 cft1st(n, a, w);
985 l = 8;
986 while ((l << 2) < n) {
987 cftmdl(n, l, a, w);
988 l <<= 2;
989 }
990 }
991 if ((l << 2) == n) {
992 for (j = 0; j < l; j += 2) {
993 j1 = j + l;
994 j2 = j1 + l;
995 j3 = j2 + l;
996 x0r = a[j] + a[j1];
997 x0i = -a[j + 1] - a[j1 + 1];
998 x1r = a[j] - a[j1];
999 x1i = -a[j + 1] + a[j1 + 1];
1000 x2r = a[j2] + a[j3];
1001 x2i = a[j2 + 1] + a[j3 + 1];
1002 x3r = a[j2] - a[j3];
1003 x3i = a[j2 + 1] - a[j3 + 1];
1004 a[j] = x0r + x2r;
1005 a[j + 1] = x0i - x2i;
1006 a[j2] = x0r - x2r;
1007 a[j2 + 1] = x0i + x2i;
1008 a[j1] = x1r - x3i;
1009 a[j1 + 1] = x1i - x3r;
1010 a[j3] = x1r + x3i;
1011 a[j3 + 1] = x1i + x3r;
1012 }
1013 }
1014 else {
1015 for (j = 0; j < l; j += 2) {
1016 j1 = j + l;
1017 x0r = a[j] - a[j1];
1018 x0i = -a[j + 1] + a[j1 + 1];
1019 a[j] += a[j1];
1020 a[j + 1] = -a[j + 1] - a[j1 + 1];
1021 a[j1] = x0r;
1022 a[j1 + 1] = x0i;
1023 }
1024 }
1025 }
1026
1027
cft1st(int n,float * a,float * w)1028 void cft1st(int n, float *a, float *w)
1029 {
1030 int j, k1, k2;
1031 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1032 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1033
1034 x0r = a[0] + a[2];
1035 x0i = a[1] + a[3];
1036 x1r = a[0] - a[2];
1037 x1i = a[1] - a[3];
1038 x2r = a[4] + a[6];
1039 x2i = a[5] + a[7];
1040 x3r = a[4] - a[6];
1041 x3i = a[5] - a[7];
1042 a[0] = x0r + x2r;
1043 a[1] = x0i + x2i;
1044 a[4] = x0r - x2r;
1045 a[5] = x0i - x2i;
1046 a[2] = x1r - x3i;
1047 a[3] = x1i + x3r;
1048 a[6] = x1r + x3i;
1049 a[7] = x1i - x3r;
1050 wk1r = w[2];
1051 x0r = a[8] + a[10];
1052 x0i = a[9] + a[11];
1053 x1r = a[8] - a[10];
1054 x1i = a[9] - a[11];
1055 x2r = a[12] + a[14];
1056 x2i = a[13] + a[15];
1057 x3r = a[12] - a[14];
1058 x3i = a[13] - a[15];
1059 a[8] = x0r + x2r;
1060 a[9] = x0i + x2i;
1061 a[12] = x2i - x0i;
1062 a[13] = x0r - x2r;
1063 x0r = x1r - x3i;
1064 x0i = x1i + x3r;
1065 a[10] = wk1r * (x0r - x0i);
1066 a[11] = wk1r * (x0r + x0i);
1067 x0r = x3i + x1r;
1068 x0i = x3r - x1i;
1069 a[14] = wk1r * (x0i - x0r);
1070 a[15] = wk1r * (x0i + x0r);
1071 k1 = 0;
1072 for (j = 16; j < n; j += 16) {
1073 k1 += 2;
1074 k2 = 2 * k1;
1075 wk2r = w[k1];
1076 wk2i = w[k1 + 1];
1077 wk1r = w[k2];
1078 wk1i = w[k2 + 1];
1079 wk3r = wk1r - 2 * wk2i * wk1i;
1080 wk3i = 2 * wk2i * wk1r - wk1i;
1081 x0r = a[j] + a[j + 2];
1082 x0i = a[j + 1] + a[j + 3];
1083 x1r = a[j] - a[j + 2];
1084 x1i = a[j + 1] - a[j + 3];
1085 x2r = a[j + 4] + a[j + 6];
1086 x2i = a[j + 5] + a[j + 7];
1087 x3r = a[j + 4] - a[j + 6];
1088 x3i = a[j + 5] - a[j + 7];
1089 a[j] = x0r + x2r;
1090 a[j + 1] = x0i + x2i;
1091 x0r -= x2r;
1092 x0i -= x2i;
1093 a[j + 4] = wk2r * x0r - wk2i * x0i;
1094 a[j + 5] = wk2r * x0i + wk2i * x0r;
1095 x0r = x1r - x3i;
1096 x0i = x1i + x3r;
1097 a[j + 2] = wk1r * x0r - wk1i * x0i;
1098 a[j + 3] = wk1r * x0i + wk1i * x0r;
1099 x0r = x1r + x3i;
1100 x0i = x1i - x3r;
1101 a[j + 6] = wk3r * x0r - wk3i * x0i;
1102 a[j + 7] = wk3r * x0i + wk3i * x0r;
1103 wk1r = w[k2 + 2];
1104 wk1i = w[k2 + 3];
1105 wk3r = wk1r - 2 * wk2r * wk1i;
1106 wk3i = 2 * wk2r * wk1r - wk1i;
1107 x0r = a[j + 8] + a[j + 10];
1108 x0i = a[j + 9] + a[j + 11];
1109 x1r = a[j + 8] - a[j + 10];
1110 x1i = a[j + 9] - a[j + 11];
1111 x2r = a[j + 12] + a[j + 14];
1112 x2i = a[j + 13] + a[j + 15];
1113 x3r = a[j + 12] - a[j + 14];
1114 x3i = a[j + 13] - a[j + 15];
1115 a[j + 8] = x0r + x2r;
1116 a[j + 9] = x0i + x2i;
1117 x0r -= x2r;
1118 x0i -= x2i;
1119 a[j + 12] = -wk2i * x0r - wk2r * x0i;
1120 a[j + 13] = -wk2i * x0i + wk2r * x0r;
1121 x0r = x1r - x3i;
1122 x0i = x1i + x3r;
1123 a[j + 10] = wk1r * x0r - wk1i * x0i;
1124 a[j + 11] = wk1r * x0i + wk1i * x0r;
1125 x0r = x1r + x3i;
1126 x0i = x1i - x3r;
1127 a[j + 14] = wk3r * x0r - wk3i * x0i;
1128 a[j + 15] = wk3r * x0i + wk3i * x0r;
1129 }
1130 }
1131
1132
cftmdl(int n,int l,float * a,float * w)1133 void cftmdl(int n, int l, float *a, float *w)
1134 {
1135 int j, j1, j2, j3, k, k1, k2, m, m2;
1136 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1137 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1138
1139 m = l << 2;
1140 for (j = 0; j < l; j += 2) {
1141 j1 = j + l;
1142 j2 = j1 + l;
1143 j3 = j2 + l;
1144 x0r = a[j] + a[j1];
1145 x0i = a[j + 1] + a[j1 + 1];
1146 x1r = a[j] - a[j1];
1147 x1i = a[j + 1] - a[j1 + 1];
1148 x2r = a[j2] + a[j3];
1149 x2i = a[j2 + 1] + a[j3 + 1];
1150 x3r = a[j2] - a[j3];
1151 x3i = a[j2 + 1] - a[j3 + 1];
1152 a[j] = x0r + x2r;
1153 a[j + 1] = x0i + x2i;
1154 a[j2] = x0r - x2r;
1155 a[j2 + 1] = x0i - x2i;
1156 a[j1] = x1r - x3i;
1157 a[j1 + 1] = x1i + x3r;
1158 a[j3] = x1r + x3i;
1159 a[j3 + 1] = x1i - x3r;
1160 }
1161 wk1r = w[2];
1162 for (j = m; j < l + m; j += 2) {
1163 j1 = j + l;
1164 j2 = j1 + l;
1165 j3 = j2 + l;
1166 x0r = a[j] + a[j1];
1167 x0i = a[j + 1] + a[j1 + 1];
1168 x1r = a[j] - a[j1];
1169 x1i = a[j + 1] - a[j1 + 1];
1170 x2r = a[j2] + a[j3];
1171 x2i = a[j2 + 1] + a[j3 + 1];
1172 x3r = a[j2] - a[j3];
1173 x3i = a[j2 + 1] - a[j3 + 1];
1174 a[j] = x0r + x2r;
1175 a[j + 1] = x0i + x2i;
1176 a[j2] = x2i - x0i;
1177 a[j2 + 1] = x0r - x2r;
1178 x0r = x1r - x3i;
1179 x0i = x1i + x3r;
1180 a[j1] = wk1r * (x0r - x0i);
1181 a[j1 + 1] = wk1r * (x0r + x0i);
1182 x0r = x3i + x1r;
1183 x0i = x3r - x1i;
1184 a[j3] = wk1r * (x0i - x0r);
1185 a[j3 + 1] = wk1r * (x0i + x0r);
1186 }
1187 k1 = 0;
1188 m2 = 2 * m;
1189 for (k = m2; k < n; k += m2) {
1190 k1 += 2;
1191 k2 = 2 * k1;
1192 wk2r = w[k1];
1193 wk2i = w[k1 + 1];
1194 wk1r = w[k2];
1195 wk1i = w[k2 + 1];
1196 wk3r = wk1r - 2 * wk2i * wk1i;
1197 wk3i = 2 * wk2i * wk1r - wk1i;
1198 for (j = k; j < l + k; j += 2) {
1199 j1 = j + l;
1200 j2 = j1 + l;
1201 j3 = j2 + l;
1202 x0r = a[j] + a[j1];
1203 x0i = a[j + 1] + a[j1 + 1];
1204 x1r = a[j] - a[j1];
1205 x1i = a[j + 1] - a[j1 + 1];
1206 x2r = a[j2] + a[j3];
1207 x2i = a[j2 + 1] + a[j3 + 1];
1208 x3r = a[j2] - a[j3];
1209 x3i = a[j2 + 1] - a[j3 + 1];
1210 a[j] = x0r + x2r;
1211 a[j + 1] = x0i + x2i;
1212 x0r -= x2r;
1213 x0i -= x2i;
1214 a[j2] = wk2r * x0r - wk2i * x0i;
1215 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
1216 x0r = x1r - x3i;
1217 x0i = x1i + x3r;
1218 a[j1] = wk1r * x0r - wk1i * x0i;
1219 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1220 x0r = x1r + x3i;
1221 x0i = x1i - x3r;
1222 a[j3] = wk3r * x0r - wk3i * x0i;
1223 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1224 }
1225 wk1r = w[k2 + 2];
1226 wk1i = w[k2 + 3];
1227 wk3r = wk1r - 2 * wk2r * wk1i;
1228 wk3i = 2 * wk2r * wk1r - wk1i;
1229 for (j = k + m; j < l + (k + m); j += 2) {
1230 j1 = j + l;
1231 j2 = j1 + l;
1232 j3 = j2 + l;
1233 x0r = a[j] + a[j1];
1234 x0i = a[j + 1] + a[j1 + 1];
1235 x1r = a[j] - a[j1];
1236 x1i = a[j + 1] - a[j1 + 1];
1237 x2r = a[j2] + a[j3];
1238 x2i = a[j2 + 1] + a[j3 + 1];
1239 x3r = a[j2] - a[j3];
1240 x3i = a[j2 + 1] - a[j3 + 1];
1241 a[j] = x0r + x2r;
1242 a[j + 1] = x0i + x2i;
1243 x0r -= x2r;
1244 x0i -= x2i;
1245 a[j2] = -wk2i * x0r - wk2r * x0i;
1246 a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
1247 x0r = x1r - x3i;
1248 x0i = x1i + x3r;
1249 a[j1] = wk1r * x0r - wk1i * x0i;
1250 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1251 x0r = x1r + x3i;
1252 x0i = x1i - x3r;
1253 a[j3] = wk3r * x0r - wk3i * x0i;
1254 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1255 }
1256 }
1257 }
1258
1259
rftfsub(int n,float * a,int nc,float * c)1260 void rftfsub(int n, float *a, int nc, float *c)
1261 {
1262 int j, k, kk, ks, m;
1263 float wkr, wki, xr, xi, yr, yi;
1264
1265 m = n >> 1;
1266 ks = 2 * nc / m;
1267 kk = 0;
1268 for (j = 2; j < m; j += 2) {
1269 k = n - j;
1270 kk += ks;
1271 wkr = 0.5 - c[nc - kk];
1272 wki = c[kk];
1273 xr = a[j] - a[k];
1274 xi = a[j + 1] + a[k + 1];
1275 yr = wkr * xr - wki * xi;
1276 yi = wkr * xi + wki * xr;
1277 a[j] -= yr;
1278 a[j + 1] -= yi;
1279 a[k] += yr;
1280 a[k + 1] -= yi;
1281 }
1282 }
1283
1284
rftbsub(int n,float * a,int nc,float * c)1285 void rftbsub(int n, float *a, int nc, float *c)
1286 {
1287 int j, k, kk, ks, m;
1288 float wkr, wki, xr, xi, yr, yi;
1289
1290 a[1] = -a[1];
1291 m = n >> 1;
1292 ks = 2 * nc / m;
1293 kk = 0;
1294 for (j = 2; j < m; j += 2) {
1295 k = n - j;
1296 kk += ks;
1297 wkr = 0.5 - c[nc - kk];
1298 wki = c[kk];
1299 xr = a[j] - a[k];
1300 xi = a[j + 1] + a[k + 1];
1301 yr = wkr * xr + wki * xi;
1302 yi = wkr * xi - wki * xr;
1303 a[j] -= yr;
1304 a[j + 1] = yi - a[j + 1];
1305 a[k] += yr;
1306 a[k + 1] = yi - a[k + 1];
1307 }
1308 a[m + 1] = -a[m + 1];
1309 }
1310
1311
dctsub(int n,float * a,int nc,float * c)1312 void dctsub(int n, float *a, int nc, float *c)
1313 {
1314 int j, k, kk, ks, m;
1315 float wkr, wki, xr;
1316
1317 m = n >> 1;
1318 ks = nc / n;
1319 kk = 0;
1320 for (j = 1; j < m; j++) {
1321 k = n - j;
1322 kk += ks;
1323 wkr = c[kk] - c[nc - kk];
1324 wki = c[kk] + c[nc - kk];
1325 xr = wki * a[j] - wkr * a[k];
1326 a[j] = wkr * a[j] + wki * a[k];
1327 a[k] = xr;
1328 }
1329 a[m] *= c[0];
1330 }
1331
1332
dstsub(int n,float * a,int nc,float * c)1333 void dstsub(int n, float *a, int nc, float *c)
1334 {
1335 int j, k, kk, ks, m;
1336 float wkr, wki, xr;
1337
1338 m = n >> 1;
1339 ks = nc / n;
1340 kk = 0;
1341 for (j = 1; j < m; j++) {
1342 k = n - j;
1343 kk += ks;
1344 wkr = c[kk] - c[nc - kk];
1345 wki = c[kk] + c[nc - kk];
1346 xr = wki * a[k] - wkr * a[j];
1347 a[k] = wkr * a[k] + wki * a[j];
1348 a[j] = xr;
1349 }
1350 a[m] *= c[0];
1351 }
1352
1353 }
1354