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