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
cdft(int n,int isgn,float * a,int * ip,float * w)289 void 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
rdft(int n,int isgn,float * a,int * ip,float * w)314 void 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.5 * (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
ddct(int n,int isgn,float * a,int * ip,float * w)361 void 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
ddst(int n,int isgn,float * a,int * ip,float * w)420 void 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
dfct(int n,float * a,float * t,int * ip,float * w)479 void 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
dfst(int n,float * a,float * t,int * ip,float * w)575 void 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
makewt(int nw,int * ip,float * w)667 void 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 = atan(1.0) / nwh;
678 w[0] = 1;
679 w[1] = 0;
680 w[nwh] = cos(delta * nwh);
681 w[nwh + 1] = w[nwh];
682 if (nwh > 2) {
683 for (j = 2; j < nwh; j += 2) {
684 x = cos(delta * j);
685 y = 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
makect(int nc,int * ip,float * c)697 void 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 = atan(1.0) / nch;
706 c[0] = cos(delta * nch);
707 c[nch] = 0.5 * c[0];
708 for (j = 1; j < nch; j++) {
709 c[j] = 0.5 * cos(delta * j);
710 c[nc - j] = 0.5 * sin(delta * j);
711 }
712 }
713 }
714
715
716 /* -------- child routines -------- */
717
718
bitrv2(int n,int * ip,float * a)719 void 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
bitrv2conj(int n,int * ip,float * a)819 void 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
cftfsub(int n,float * a,float * w)928 void 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
cftbsub(int n,float * a,float * w)980 void 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
cft1st(int n,float * a,float * w)1032 void 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
cftmdl(int n,int l,float * a,float * w)1137 void 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
rftfsub(int n,float * a,int nc,float * c)1264 void 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.5 - 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
rftbsub(int n,float * a,int nc,float * c)1289 void 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.5 - 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
dctsub(int n,float * a,int nc,float * c)1316 void 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
dstsub(int n,float * a,int nc,float * c)1337 void 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
1357