CLAM-Development  1.4.0
FFT_ooura.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 #include "FFT_ooura.hxx"
24 
25 #include "Assert.hxx"
26 #include "Audio.hxx"
27 #include "Spectrum.hxx"
28 #include "SpectrumConfig.hxx"
29 #include "CLAM_Math.hxx"
30 #include "ProcessingFactory.hxx"
31 
32 namespace CLAM
33 {
34 
35 namespace Hidden
36 {
37  static const char * metadata[] = {
38  "key", "FFT_ooura",
39  "category", "Analysis",
40  "description", "FFT_ooura",
41  0
42  };
43  static FactoryRegistrator<ProcessingFactory, FFT_ooura> reg = metadata;
44 }
45 
46  bool FFT_ooura::ConcreteConfigure(const ProcessingConfig& c) {
47  int oldSize=mSize;
49  if ( !isPowerOfTwo( mSize ) )
50  {
51  AddConfigErrorMessage("Configure failed: FFT Ooura algorithm only works for input buffers that are a power of two!");
52  return false;
53  }
54  if (mSize>0) {
55  if (mSize != oldSize) {
56  ReleaseMemory();
57  SetupMemory();
58  }
59  return true;
60  }
61  ReleaseMemory();
62  return false;
63  }
64 
65  void FFT_ooura::ReleaseMemory() {
66  if (ip) { delete[] ip; ip = 0; }
67  if (w) { delete[] w; w = 0; }
68  }
69 
70  void FFT_ooura::SetupMemory() {
71  int ipSize = (int)(2+(1<<(int)(log(mSize/2+0.5)/log(2.0))/2));
72  ip = new int[ipSize];
73  for (int i=0; i<ipSize; i++) ip[i] = 0;
74 
75  int wSize = (int)(mSize*5/8-1);
76  w = new TData[wSize];
77  for (int i=0; i<wSize; i++) w[i] = 0;
78  }
79 
81  : ip(0), w(0)
82  {
84  }
85 
87  : ip(0), w(0)
88  {
89  Configure(c);
90  };
91 
93  {
94  ReleaseMemory();
95  }
96 
97  bool FFT_ooura::Do()
98  {
100  bool toReturn = Do(mInput.GetAudio(), mOutput.GetData());
101  mInput.Consume();
102  mOutput.Produce();
103  return toReturn;
104  }
105 
106  bool FFT_ooura::Do(const Audio& in, Spectrum &out){
107  TData *inbuffer;
108 
110  "FFT_ooura: Do(): Not in execution mode");
112  "FFT_ooura: Do(): Not a power of two");
113 
114  out.SetSpectralRange(in.GetSampleRate()/2);
115 
116  switch(mState) {
117  case sComplex:
118  inbuffer = in.GetBuffer().GetPtr();
119  // Buffer dump. This is a kludge; the right way to do this
120  // is using a non-inplace version of rdft (which would
121  // not reduce performance).
122  for (int i=0; i<mSize; i++)
123  fftbuffer[i] = inbuffer[i];
124  rdft(mSize, 1, fftbuffer, ip, w);
125  ToComplex(out);
126  break;
127  case sComplexSync:
128  inbuffer = in.GetBuffer().GetPtr();
129  // Buffer dump. This is a kludge; the right way to do this
130  // is using a non-inplace version of rdft (which would
131  // not reduce performance).
132  for (int i=0; i<mSize; i++)
133  fftbuffer[i]=inbuffer[i];
134  rdft(mSize, 1, fftbuffer, ip, w);
135  ToComplex(out);
137  break;
138  case sOther:
139  CheckTypes(in,out);
140  inbuffer = in.GetBuffer().GetPtr();
141  // Buffer dump. This is a kludge; the right way to do this
142  // is using a non-inplace version of rdft (which would
143  // not reduce performance).
144  for (int i=0; i<mSize; i++)
145  fftbuffer[i]=inbuffer[i];
146  rdft(mSize, 1, fftbuffer, ip, w);
147 
148  ToOther(out);
149  break;
150  default:
151  CLAM_ASSERT(false, "FFT_ooura: Do(): Inconsistent state");
152  }
153 
154  return true;
155  }
156 
157 
158  void FFT_ooura::ToComplex(Spectrum &out) {
159  Array<Complex>& outbuffer = out.GetComplexArray();
160 
161  outbuffer[0].SetReal(fftbuffer[0]); // Real Values
162  outbuffer[0].SetImag(0); // Real Values
163  outbuffer[mSize/2].SetReal(fftbuffer[1]);
164  outbuffer[mSize/2].SetImag(0);
165 
166  for (int i=1; i<mSize/2; i++) {
167  outbuffer[i].SetReal(fftbuffer[2*i]);
168  outbuffer[i].SetImag(-fftbuffer[2*i+1]);
169  }
170 
171  outbuffer.SetSize(mSize/2+1);
172  }
173 
174  void FFT_ooura::ToOther(Spectrum &out)
175  {
176  if(out.HasComplexArray()) {
177  ToComplex(out);
178  out.SynchronizeTo(mComplexflags);
179  }
180  else {
181  SpecTypeFlags flags;
182  out.GetType(flags);
183 
184 
185  ToComplex(mComplexSpectrum);
186  out.SynchronizeTo(mComplexSpectrum);
187  }
188  }
189 
190 
192  // Original Ooura FFT functions, modified to accept TData instead of double
193 
194 void FFT_ooura::rdft(int n, int isgn, TData *a, int *ip, TData *w)
195 {
196  int nw, nc;
197  TData xi;
198 
199  nw = ip[0];
200  if (n > (nw << 2)) {
201  nw = n >> 2;
202  makewt(nw, ip, w);
203  }
204  nc = ip[1];
205  if (n > (nc << 2)) {
206  nc = n >> 2;
207  makect(nc, ip, w + nw);
208  }
209  if (isgn >= 0) {
210  if (n > 4) {
211  bitrv2(n, ip + 2, a);
212  cftfsub(n, a, w);
213  rftfsub(n, a, nc, w + nw);
214  } else if (n == 4) {
215  cftfsub(n, a, w);
216  }
217  xi = a[0] - a[1];
218  a[0] += a[1];
219  a[1] = xi;
220  } else {
221  a[1] = 0.5 * (a[0] - a[1]);
222  a[0] -= a[1];
223  if (n > 4) {
224  rftbsub(n, a, nc, w + nw);
225  bitrv2(n, ip + 2, a);
226  cftbsub(n, a, w);
227  } else if (n == 4) {
228  cftfsub(n, a, w);
229  }
230  }
231 }
232 
233  // workspace setup functions...
234 
235 void FFT_ooura::makewt(int nw, int *ip, TData *w)
236 {
237  int j, nwh;
238  TData delta, x, y;
239 
240  ip[0] = nw;
241  ip[1] = 1;
242  if (nw > 2) {
243  nwh = nw >> 1;
244  delta = CLAM_atan(1.0) / nwh;
245  w[0] = 1;
246  w[1] = 0;
247  w[nwh] = CLAM_cos(delta * nwh);
248  w[nwh + 1] = w[nwh];
249  if (nwh > 2) {
250  for (j = 2; j < nwh; j += 2) {
251  x = CLAM_cos(delta * j);
252  y = CLAM_sin(delta * j);
253  w[j] = x;
254  w[j + 1] = y;
255  w[nw - j] = y;
256  w[nw - j + 1] = x;
257  }
258  bitrv2(nw, ip + 2, w);
259  }
260  }
261 }
262 
263 
264 void FFT_ooura::makect(int nc, int *ip, TData *c)
265 {
266  int j, nch;
267  TData delta;
268 
269  ip[1] = nc;
270  if (nc > 1) {
271  nch = nc >> 1;
272  delta = CLAM_atan(1.0) / nch;
273  c[0] = CLAM_cos(delta * nch);
274  c[nch] = 0.5 * c[0];
275  for (j = 1; j < nch; j++) {
276  c[j] = 0.5 * CLAM_cos(delta * j);
277  c[nc - j] = 0.5 * CLAM_sin(delta * j);
278  }
279  }
280 }
281 
282 void FFT_ooura::bitrv2(int n, int *ip, TData *a)
283 {
284  int j, j1, k, k1, l, m, m2;
285  TData xr, xi, yr, yi;
286 
287  ip[0] = 0;
288  l = n;
289  m = 1;
290  while ((m << 3) < l) {
291  l >>= 1;
292  for (j = 0; j < m; j++) {
293  ip[m + j] = ip[j] + l;
294  }
295  m <<= 1;
296  }
297  m2 = 2 * m;
298  if ((m << 3) == l) {
299  for (k = 0; k < m; k++) {
300  for (j = 0; j < k; j++) {
301  j1 = 2 * j + ip[k];
302  k1 = 2 * k + ip[j];
303  xr = a[j1];
304  xi = a[j1 + 1];
305  yr = a[k1];
306  yi = a[k1 + 1];
307  a[j1] = yr;
308  a[j1 + 1] = yi;
309  a[k1] = xr;
310  a[k1 + 1] = xi;
311  j1 += m2;
312  k1 += 2 * m2;
313  xr = a[j1];
314  xi = a[j1 + 1];
315  yr = a[k1];
316  yi = a[k1 + 1];
317  a[j1] = yr;
318  a[j1 + 1] = yi;
319  a[k1] = xr;
320  a[k1 + 1] = xi;
321  j1 += m2;
322  k1 -= m2;
323  xr = a[j1];
324  xi = a[j1 + 1];
325  yr = a[k1];
326  yi = a[k1 + 1];
327  a[j1] = yr;
328  a[j1 + 1] = yi;
329  a[k1] = xr;
330  a[k1 + 1] = xi;
331  j1 += m2;
332  k1 += 2 * m2;
333  xr = a[j1];
334  xi = a[j1 + 1];
335  yr = a[k1];
336  yi = a[k1 + 1];
337  a[j1] = yr;
338  a[j1 + 1] = yi;
339  a[k1] = xr;
340  a[k1 + 1] = xi;
341  }
342  j1 = 2 * k + m2 + ip[k];
343  k1 = j1 + m2;
344  xr = a[j1];
345  xi = a[j1 + 1];
346  yr = a[k1];
347  yi = a[k1 + 1];
348  a[j1] = yr;
349  a[j1 + 1] = yi;
350  a[k1] = xr;
351  a[k1 + 1] = xi;
352  }
353  } else {
354  for (k = 1; k < m; k++) {
355  for (j = 0; j < k; j++) {
356  j1 = 2 * j + ip[k];
357  k1 = 2 * k + ip[j];
358  xr = a[j1];
359  xi = a[j1 + 1];
360  yr = a[k1];
361  yi = a[k1 + 1];
362  a[j1] = yr;
363  a[j1 + 1] = yi;
364  a[k1] = xr;
365  a[k1 + 1] = xi;
366  j1 += m2;
367  k1 += m2;
368  xr = a[j1];
369  xi = a[j1 + 1];
370  yr = a[k1];
371  yi = a[k1 + 1];
372  a[j1] = yr;
373  a[j1 + 1] = yi;
374  a[k1] = xr;
375  a[k1 + 1] = xi;
376  }
377  }
378  }
379 }
380 
381 void FFT_ooura::cftfsub(int n, TData *a, TData *w)
382 {
383  int j, j1, j2, j3, l;
384  TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
385 
386  l = 2;
387  if (n > 8) {
388  cft1st(n, a, w);
389  l = 8;
390  while ((l << 2) < n) {
391  cftmdl(n, l, a, w);
392  l <<= 2;
393  }
394  }
395  if ((l << 2) == n) {
396  for (j = 0; j < l; j += 2) {
397  j1 = j + l;
398  j2 = j1 + l;
399  j3 = j2 + l;
400  x0r = a[j] + a[j1];
401  x0i = a[j + 1] + a[j1 + 1];
402  x1r = a[j] - a[j1];
403  x1i = a[j + 1] - a[j1 + 1];
404  x2r = a[j2] + a[j3];
405  x2i = a[j2 + 1] + a[j3 + 1];
406  x3r = a[j2] - a[j3];
407  x3i = a[j2 + 1] - a[j3 + 1];
408  a[j] = x0r + x2r;
409  a[j + 1] = x0i + x2i;
410  a[j2] = x0r - x2r;
411  a[j2 + 1] = x0i - x2i;
412  a[j1] = x1r - x3i;
413  a[j1 + 1] = x1i + x3r;
414  a[j3] = x1r + x3i;
415  a[j3 + 1] = x1i - x3r;
416  }
417  } else {
418  for (j = 0; j < l; j += 2) {
419  j1 = j + l;
420  x0r = a[j] - a[j1];
421  x0i = a[j + 1] - a[j1 + 1];
422  a[j] += a[j1];
423  a[j + 1] += a[j1 + 1];
424  a[j1] = x0r;
425  a[j1 + 1] = x0i;
426  }
427  }
428 }
429 
430 void FFT_ooura::cftbsub(int n, TData *a, TData *w)
431 {
432  int j, j1, j2, j3, l;
433  TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
434 
435  l = 2;
436  if (n > 8) {
437  cft1st(n, a, w);
438  l = 8;
439  while ((l << 2) < n) {
440  cftmdl(n, l, a, w);
441  l <<= 2;
442  }
443  }
444  if ((l << 2) == n) {
445  for (j = 0; j < l; j += 2) {
446  j1 = j + l;
447  j2 = j1 + l;
448  j3 = j2 + l;
449  x0r = a[j] + a[j1];
450  x0i = -a[j + 1] - a[j1 + 1];
451  x1r = a[j] - a[j1];
452  x1i = -a[j + 1] + a[j1 + 1];
453  x2r = a[j2] + a[j3];
454  x2i = a[j2 + 1] + a[j3 + 1];
455  x3r = a[j2] - a[j3];
456  x3i = a[j2 + 1] - a[j3 + 1];
457  a[j] = x0r + x2r;
458  a[j + 1] = x0i - x2i;
459  a[j2] = x0r - x2r;
460  a[j2 + 1] = x0i + x2i;
461  a[j1] = x1r - x3i;
462  a[j1 + 1] = x1i - x3r;
463  a[j3] = x1r + x3i;
464  a[j3 + 1] = x1i + x3r;
465  }
466  } else {
467  for (j = 0; j < l; j += 2) {
468  j1 = j + l;
469  x0r = a[j] - a[j1];
470  x0i = -a[j + 1] + a[j1 + 1];
471  a[j] += a[j1];
472  a[j + 1] = -a[j + 1] - a[j1 + 1];
473  a[j1] = x0r;
474  a[j1 + 1] = x0i;
475  }
476  }
477 }
478 
479 void FFT_ooura::cft1st(int n, TData *a, TData *w)
480 {
481  int j, k1, k2;
482  TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
483  TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
484 
485  x0r = a[0] + a[2];
486  x0i = a[1] + a[3];
487  x1r = a[0] - a[2];
488  x1i = a[1] - a[3];
489  x2r = a[4] + a[6];
490  x2i = a[5] + a[7];
491  x3r = a[4] - a[6];
492  x3i = a[5] - a[7];
493  a[0] = x0r + x2r;
494  a[1] = x0i + x2i;
495  a[4] = x0r - x2r;
496  a[5] = x0i - x2i;
497  a[2] = x1r - x3i;
498  a[3] = x1i + x3r;
499  a[6] = x1r + x3i;
500  a[7] = x1i - x3r;
501  wk1r = w[2];
502  x0r = a[8] + a[10];
503  x0i = a[9] + a[11];
504  x1r = a[8] - a[10];
505  x1i = a[9] - a[11];
506  x2r = a[12] + a[14];
507  x2i = a[13] + a[15];
508  x3r = a[12] - a[14];
509  x3i = a[13] - a[15];
510  a[8] = x0r + x2r;
511  a[9] = x0i + x2i;
512  a[12] = x2i - x0i;
513  a[13] = x0r - x2r;
514  x0r = x1r - x3i;
515  x0i = x1i + x3r;
516  a[10] = wk1r * (x0r - x0i);
517  a[11] = wk1r * (x0r + x0i);
518  x0r = x3i + x1r;
519  x0i = x3r - x1i;
520  a[14] = wk1r * (x0i - x0r);
521  a[15] = wk1r * (x0i + x0r);
522  k1 = 0;
523  for (j = 16; j < n; j += 16) {
524  k1 += 2;
525  k2 = 2 * k1;
526  wk2r = w[k1];
527  wk2i = w[k1 + 1];
528  wk1r = w[k2];
529  wk1i = w[k2 + 1];
530  wk3r = wk1r - 2 * wk2i * wk1i;
531  wk3i = 2 * wk2i * wk1r - wk1i;
532  x0r = a[j] + a[j + 2];
533  x0i = a[j + 1] + a[j + 3];
534  x1r = a[j] - a[j + 2];
535  x1i = a[j + 1] - a[j + 3];
536  x2r = a[j + 4] + a[j + 6];
537  x2i = a[j + 5] + a[j + 7];
538  x3r = a[j + 4] - a[j + 6];
539  x3i = a[j + 5] - a[j + 7];
540  a[j] = x0r + x2r;
541  a[j + 1] = x0i + x2i;
542  x0r -= x2r;
543  x0i -= x2i;
544  a[j + 4] = wk2r * x0r - wk2i * x0i;
545  a[j + 5] = wk2r * x0i + wk2i * x0r;
546  x0r = x1r - x3i;
547  x0i = x1i + x3r;
548  a[j + 2] = wk1r * x0r - wk1i * x0i;
549  a[j + 3] = wk1r * x0i + wk1i * x0r;
550  x0r = x1r + x3i;
551  x0i = x1i - x3r;
552  a[j + 6] = wk3r * x0r - wk3i * x0i;
553  a[j + 7] = wk3r * x0i + wk3i * x0r;
554  wk1r = w[k2 + 2];
555  wk1i = w[k2 + 3];
556  wk3r = wk1r - 2 * wk2r * wk1i;
557  wk3i = 2 * wk2r * wk1r - wk1i;
558  x0r = a[j + 8] + a[j + 10];
559  x0i = a[j + 9] + a[j + 11];
560  x1r = a[j + 8] - a[j + 10];
561  x1i = a[j + 9] - a[j + 11];
562  x2r = a[j + 12] + a[j + 14];
563  x2i = a[j + 13] + a[j + 15];
564  x3r = a[j + 12] - a[j + 14];
565  x3i = a[j + 13] - a[j + 15];
566  a[j + 8] = x0r + x2r;
567  a[j + 9] = x0i + x2i;
568  x0r -= x2r;
569  x0i -= x2i;
570  a[j + 12] = -wk2i * x0r - wk2r * x0i;
571  a[j + 13] = -wk2i * x0i + wk2r * x0r;
572  x0r = x1r - x3i;
573  x0i = x1i + x3r;
574  a[j + 10] = wk1r * x0r - wk1i * x0i;
575  a[j + 11] = wk1r * x0i + wk1i * x0r;
576  x0r = x1r + x3i;
577  x0i = x1i - x3r;
578  a[j + 14] = wk3r * x0r - wk3i * x0i;
579  a[j + 15] = wk3r * x0i + wk3i * x0r;
580  }
581 }
582 
583 
584 void FFT_ooura::cftmdl(int n, int l, TData *a, TData *w)
585 {
586  int j, j1, j2, j3, k, k1, k2, m, m2;
587  TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
588  TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
589 
590  m = l << 2;
591  for (j = 0; j < l; j += 2) {
592  j1 = j + l;
593  j2 = j1 + l;
594  j3 = j2 + l;
595  x0r = a[j] + a[j1];
596  x0i = a[j + 1] + a[j1 + 1];
597  x1r = a[j] - a[j1];
598  x1i = a[j + 1] - a[j1 + 1];
599  x2r = a[j2] + a[j3];
600  x2i = a[j2 + 1] + a[j3 + 1];
601  x3r = a[j2] - a[j3];
602  x3i = a[j2 + 1] - a[j3 + 1];
603  a[j] = x0r + x2r;
604  a[j + 1] = x0i + x2i;
605  a[j2] = x0r - x2r;
606  a[j2 + 1] = x0i - x2i;
607  a[j1] = x1r - x3i;
608  a[j1 + 1] = x1i + x3r;
609  a[j3] = x1r + x3i;
610  a[j3 + 1] = x1i - x3r;
611  }
612  wk1r = w[2];
613  for (j = m; j < l + m; j += 2) {
614  j1 = j + l;
615  j2 = j1 + l;
616  j3 = j2 + l;
617  x0r = a[j] + a[j1];
618  x0i = a[j + 1] + a[j1 + 1];
619  x1r = a[j] - a[j1];
620  x1i = a[j + 1] - a[j1 + 1];
621  x2r = a[j2] + a[j3];
622  x2i = a[j2 + 1] + a[j3 + 1];
623  x3r = a[j2] - a[j3];
624  x3i = a[j2 + 1] - a[j3 + 1];
625  a[j] = x0r + x2r;
626  a[j + 1] = x0i + x2i;
627  a[j2] = x2i - x0i;
628  a[j2 + 1] = x0r - x2r;
629  x0r = x1r - x3i;
630  x0i = x1i + x3r;
631  a[j1] = wk1r * (x0r - x0i);
632  a[j1 + 1] = wk1r * (x0r + x0i);
633  x0r = x3i + x1r;
634  x0i = x3r - x1i;
635  a[j3] = wk1r * (x0i - x0r);
636  a[j3 + 1] = wk1r * (x0i + x0r);
637  }
638  k1 = 0;
639  m2 = 2 * m;
640  for (k = m2; k < n; k += m2) {
641  k1 += 2;
642  k2 = 2 * k1;
643  wk2r = w[k1];
644  wk2i = w[k1 + 1];
645  wk1r = w[k2];
646  wk1i = w[k2 + 1];
647  wk3r = wk1r - 2 * wk2i * wk1i;
648  wk3i = 2 * wk2i * wk1r - wk1i;
649  for (j = k; j < l + k; j += 2) {
650  j1 = j + l;
651  j2 = j1 + l;
652  j3 = j2 + l;
653  x0r = a[j] + a[j1];
654  x0i = a[j + 1] + a[j1 + 1];
655  x1r = a[j] - a[j1];
656  x1i = a[j + 1] - a[j1 + 1];
657  x2r = a[j2] + a[j3];
658  x2i = a[j2 + 1] + a[j3 + 1];
659  x3r = a[j2] - a[j3];
660  x3i = a[j2 + 1] - a[j3 + 1];
661  a[j] = x0r + x2r;
662  a[j + 1] = x0i + x2i;
663  x0r -= x2r;
664  x0i -= x2i;
665  a[j2] = wk2r * x0r - wk2i * x0i;
666  a[j2 + 1] = wk2r * x0i + wk2i * x0r;
667  x0r = x1r - x3i;
668  x0i = x1i + x3r;
669  a[j1] = wk1r * x0r - wk1i * x0i;
670  a[j1 + 1] = wk1r * x0i + wk1i * x0r;
671  x0r = x1r + x3i;
672  x0i = x1i - x3r;
673  a[j3] = wk3r * x0r - wk3i * x0i;
674  a[j3 + 1] = wk3r * x0i + wk3i * x0r;
675  }
676  wk1r = w[k2 + 2];
677  wk1i = w[k2 + 3];
678  wk3r = wk1r - 2 * wk2r * wk1i;
679  wk3i = 2 * wk2r * wk1r - wk1i;
680  for (j = k + m; j < l + (k + m); j += 2) {
681  j1 = j + l;
682  j2 = j1 + l;
683  j3 = j2 + l;
684  x0r = a[j] + a[j1];
685  x0i = a[j + 1] + a[j1 + 1];
686  x1r = a[j] - a[j1];
687  x1i = a[j + 1] - a[j1 + 1];
688  x2r = a[j2] + a[j3];
689  x2i = a[j2 + 1] + a[j3 + 1];
690  x3r = a[j2] - a[j3];
691  x3i = a[j2 + 1] - a[j3 + 1];
692  a[j] = x0r + x2r;
693  a[j + 1] = x0i + x2i;
694  x0r -= x2r;
695  x0i -= x2i;
696  a[j2] = -wk2i * x0r - wk2r * x0i;
697  a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
698  x0r = x1r - x3i;
699  x0i = x1i + x3r;
700  a[j1] = wk1r * x0r - wk1i * x0i;
701  a[j1 + 1] = wk1r * x0i + wk1i * x0r;
702  x0r = x1r + x3i;
703  x0i = x1i - x3r;
704  a[j3] = wk3r * x0r - wk3i * x0i;
705  a[j3 + 1] = wk3r * x0i + wk3i * x0r;
706  }
707  }
708 }
709 
710 
711 void FFT_ooura::rftfsub(int n, TData *a, int nc, TData *c)
712 {
713  int j, k, kk, ks, m;
714  TData wkr, wki, xr, xi, yr, yi;
715 
716  m = n >> 1;
717  ks = 2 * nc / m;
718  kk = 0;
719  for (j = 2; j < m; j += 2) {
720  k = n - j;
721  kk += ks;
722  wkr = 0.5 - c[nc - kk];
723  wki = c[kk];
724  xr = a[j] - a[k];
725  xi = a[j + 1] + a[k + 1];
726  yr = wkr * xr - wki * xi;
727  yi = wkr * xi + wki * xr;
728  a[j] -= yr;
729  a[j + 1] -= yi;
730  a[k] += yr;
731  a[k + 1] -= yi;
732  }
733 }
734 
735 
736 void FFT_ooura::rftbsub(int n, TData *a, int nc, TData *c)
737 {
738  int j, k, kk, ks, m;
739  TData wkr, wki, xr, xi, yr, yi;
740 
741  a[1] = -a[1];
742  m = n >> 1;
743  ks = 2 * nc / m;
744  kk = 0;
745  for (j = 2; j < m; j += 2) {
746  k = n - j;
747  kk += ks;
748  wkr = 0.5 - c[nc - kk];
749  wki = c[kk];
750  xr = a[j] - a[k];
751  xi = a[j + 1] + a[k + 1];
752  yr = wkr * xr + wki * xi;
753  yi = wkr * xi - wki * xr;
754  a[j] -= yr;
755  a[j + 1] = yi - a[j + 1];
756  a[k] += yr;
757  a[k + 1] = yi - a[k + 1];
758  }
759  a[m + 1] = -a[m + 1];
760 }
761 
762 
763 }; // CLAM
764