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