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