CSL  5.2
FFTReal.cpp
Go to the documentation of this file.
1 /*****************************************************************************
2 * *
3 * DIGITAL SIGNAL PROCESSING TOOLS *
4 * Version 1.01, 1999/11/07 *
5 * (c) 1999 - Laurent de Soras *
6 * *
7 * FFTReal.cpp *
8 * Fourier transformation of real number arrays. *
9 * Portable ISO C++ *
10 * *
11 * Tab = 3 *
12 *****************************************************************************/
13 
14 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
15 
16 #include "FFTReal.h"
17 #include <cassert>
18 #include <cmath>
19 
20 //#include <stdio.h>
21 
22 #if defined (_MSC_VER)
23 #pragma pack (push, 8)
24 #endif // _MSC_VER
25 
26 //using namespace csl;
27 
28 /*\\\ PUBLIC MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
29 
30 /*==========================================================================*/
31 /* Name: Constructor */
32 /* Input parameters: */
33 /* - length: length of the array on which we want to do a FFT. */
34 /* Range: power of 2 only, > 0. */
35 /* Throws: std::bad_alloc, anything */
36 /*==========================================================================*/
37 
38 FFTReal::FFTReal (const long length)
39 : _bit_rev_lut (int (floor (logf (length) / logf (2) + 0.5f)))
40 , _trigo_lut (int (floor (logf (length) / logf (2) + 0.5f)))
41 , _sqrt2_2 (flt_t (sqrtf (2) * 0.5f))
42 , _length (length)
43 , _nbr_bits (int (floor (logf (length) / logf (2) + 0.5f)))
44 {
45  assert ((1L << _nbr_bits) == length);
46 
47  _buffer_ptr = 0;
48  if (_nbr_bits > 2)
49  {
50  _buffer_ptr = new flt_t [_length];
51  }
52 }
53 
54 /*==========================================================================*/
55 /* Name: Destructor */
56 /*==========================================================================*/
57 
58 FFTReal::~FFTReal (void)
59 {
60 // printf("\t\t~FFTReal\n");
61  delete [] _buffer_ptr;
62  _buffer_ptr = 0;
63 }
64 
65 /*==========================================================================*/
66 /* Name: do_fft */
67 /* Description: Compute the FFT of the array. */
68 /* Input parameters: */
69 /* - x: pointer on the source array (time). */
70 /* Output parameters: */
71 /* - f: pointer on the destination array (frequencies). */
72 /* f [0...length(x)/2] = real values, */
73 /* f [length(x)/2+1...length(x)-1] = imaginary values of */
74 /* coefficents 1...length(x)/2-1. */
75 /* Throws: Nothing */
76 /*==========================================================================*/
77 
78 void FFTReal::do_fft (flt_t f [], const flt_t x []) const
79 {
80 
81 /*______________________________________________
82  *
83  * General case
84  *______________________________________________
85  */
86 
87  if (_nbr_bits > 2)
88  {
89  flt_t * sf;
90  flt_t * df;
91 
92  if (_nbr_bits & 1)
93  {
94  df = _buffer_ptr;
95  sf = f;
96  }
97  else
98  {
99  df = f;
100  sf = _buffer_ptr;
101  }
102 
103  /* Do the transformation in several pass */
104  {
105  int pass;
106  long nbr_coef;
107  long h_nbr_coef;
108  long d_nbr_coef;
109  long coef_index;
110 
111  /* First and second pass at once */
112  {
113  const long * const bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
114  coef_index = 0;
115  do
116  {
117  const long rev_index_0 = bit_rev_lut_ptr [coef_index];
118  const long rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
119  const long rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
120  const long rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
121 
122  flt_t * const df2 = df + coef_index;
123  df2 [1] = x [rev_index_0] - x [rev_index_1];
124  df2 [3] = x [rev_index_2] - x [rev_index_3];
125 
126  const flt_t sf_0 = x [rev_index_0] + x [rev_index_1];
127  const flt_t sf_2 = x [rev_index_2] + x [rev_index_3];
128 
129  df2 [0] = sf_0 + sf_2;
130  df2 [2] = sf_0 - sf_2;
131 
132  coef_index += 4;
133  }
134  while (coef_index < _length);
135  }
136 
137  /* Third pass */
138  {
139  coef_index = 0;
140  const flt_t sqrt2_2 = _sqrt2_2;
141  do
142  {
143  flt_t v;
144 
145  sf [coef_index] = df [coef_index] + df [coef_index + 4];
146  sf [coef_index + 4] = df [coef_index] - df [coef_index + 4];
147  sf [coef_index + 2] = df [coef_index + 2];
148  sf [coef_index + 6] = df [coef_index + 6];
149 
150  v = (df [coef_index + 5] - df [coef_index + 7]) * sqrt2_2;
151  sf [coef_index + 1] = df [coef_index + 1] + v;
152  sf [coef_index + 3] = df [coef_index + 1] - v;
153 
154  v = (df [coef_index + 5] + df [coef_index + 7]) * sqrt2_2;
155  sf [coef_index + 5] = v + df [coef_index + 3];
156  sf [coef_index + 7] = v - df [coef_index + 3];
157 
158  coef_index += 8;
159  }
160  while (coef_index < _length);
161  }
162 
163  /* Next pass */
164  for (pass = 3; pass < _nbr_bits; ++pass)
165  {
166  coef_index = 0;
167  nbr_coef = 1 << pass;
168  h_nbr_coef = nbr_coef >> 1;
169  d_nbr_coef = nbr_coef << 1;
170  const flt_t * const cos_ptr = _trigo_lut.get_ptr (pass);
171  do
172  {
173  long i;
174  const flt_t * const sf1r = sf + coef_index;
175  const flt_t * const sf2r = sf1r + nbr_coef;
176  flt_t * const dfr = df + coef_index;
177  flt_t * const dfi = dfr + nbr_coef;
178 
179  /* Extreme coefficients are always real */
180  dfr [0] = sf1r [0] + sf2r [0];
181  dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
182  dfr [h_nbr_coef] = sf1r [h_nbr_coef];
183  dfi [h_nbr_coef] = sf2r [h_nbr_coef];
184 
185  /* Others are conjugate complex numbers */
186  const flt_t * const sf1i = sf1r + h_nbr_coef;
187  const flt_t * const sf2i = sf1i + nbr_coef;
188  for (i = 1; i < h_nbr_coef; ++ i)
189  {
190  const flt_t c = cos_ptr [i]; // cos (i*PI/nbr_coef);
191  const flt_t s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
192  flt_t v;
193 
194  v = sf2r [i] * c - sf2i [i] * s;
195  dfr [i] = sf1r [i] + v;
196  dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
197 
198  v = sf2r [i] * s + sf2i [i] * c;
199  dfi [i] = v + sf1i [i];
200  dfi [nbr_coef - i] = v - sf1i [i];
201  }
202 
203  coef_index += d_nbr_coef;
204  }
205  while (coef_index < _length);
206 
207  /* Prepare to the next pass */
208  {
209  flt_t * const temp_ptr = df;
210  df = sf;
211  sf = temp_ptr;
212  }
213  }
214  }
215  }
216 
217 /*______________________________________________
218  *
219  * Special cases
220  *______________________________________________
221  */
222 
223  /* 4-point FFT */
224  else if (_nbr_bits == 2)
225  {
226  f [1] = x [0] - x [2];
227  f [3] = x [1] - x [3];
228 
229  const flt_t b_0 = x [0] + x [2];
230  const flt_t b_2 = x [1] + x [3];
231 
232  f [0] = b_0 + b_2;
233  f [2] = b_0 - b_2;
234  }
235 
236  /* 2-point FFT */
237  else if (_nbr_bits == 1)
238  {
239  f [0] = x [0] + x [1];
240  f [1] = x [0] - x [1];
241  }
242 
243  /* 1-point FFT */
244  else
245  {
246  f [0] = x [0];
247  }
248 }
249 
250 /*==========================================================================*/
251 /* Name: do_ifft */
252 /* Description: Compute the inverse FFT of the array. Notice that */
253 /* IFFT (FFT (x)) = x * length (x). Data must be */
254 /* post-scaled. */
255 /* Input parameters: */
256 /* - f: pointer on the source array (frequencies). */
257 /* f [0...length(x)/2] = real values, */
258 /* f [length(x)/2+1...length(x)] = imaginary values of */
259 /* coefficents 1...length(x)-1. */
260 /* Output parameters: */
261 /* - x: pointer on the destination array (time). */
262 /* Throws: Nothing */
263 /*==========================================================================*/
264 
265 void FFTReal::do_ifft (const flt_t f [], flt_t x []) const
266 {
267 
268 /*______________________________________________
269  *
270  * General case
271  *______________________________________________
272  */
273 
274  if (_nbr_bits > 2)
275  {
276  flt_t * sf = const_cast <flt_t *> (f);
277  flt_t * df;
278  flt_t * df_temp;
279 
280  if (_nbr_bits & 1)
281  {
282  df = _buffer_ptr;
283  df_temp = x;
284  }
285  else
286  {
287  df = x;
288  df_temp = _buffer_ptr;
289  }
290 
291  /* Do the transformation in several pass */
292  {
293  int pass;
294  long nbr_coef;
295  long h_nbr_coef;
296  long d_nbr_coef;
297  long coef_index;
298 
299  /* First pass */
300  for (pass = _nbr_bits - 1; pass >= 3; --pass)
301  {
302  coef_index = 0;
303  nbr_coef = 1 << pass;
304  h_nbr_coef = nbr_coef >> 1;
305  d_nbr_coef = nbr_coef << 1;
306  const flt_t *const cos_ptr = _trigo_lut.get_ptr (pass);
307  do
308  {
309  long i;
310  const flt_t * const sfr = sf + coef_index;
311  const flt_t * const sfi = sfr + nbr_coef;
312  flt_t * const df1r = df + coef_index;
313  flt_t * const df2r = df1r + nbr_coef;
314 
315  /* Extreme coefficients are always real */
316  df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
317  df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
318  df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
319  df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
320 
321  /* Others are conjugate complex numbers */
322  flt_t * const df1i = df1r + h_nbr_coef;
323  flt_t * const df2i = df1i + nbr_coef;
324  for (i = 1; i < h_nbr_coef; ++ i)
325  {
326  df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
327  df1i [i] = sfi [i] - sfi [nbr_coef - i];
328 
329  const flt_t c = cos_ptr [i]; // cos (i*PI/nbr_coef);
330  const flt_t s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
331  const flt_t vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
332  const flt_t vi = sfi [i] + sfi [nbr_coef - i];
333 
334  df2r [i] = vr * c + vi * s;
335  df2i [i] = vi * c - vr * s;
336  }
337  coef_index += d_nbr_coef;
338  }
339  while (coef_index < _length);
340 
341  /* Prepare to the next pass */
342  if (pass < _nbr_bits - 1)
343  {
344  flt_t * const temp_ptr = df;
345  df = sf;
346  sf = temp_ptr;
347  }
348  else
349  {
350  sf = df;
351  df = df_temp;
352  }
353  }
354 
355  /* Antepenultimate pass */
356  {
357  const flt_t sqrt2_2 = _sqrt2_2;
358  coef_index = 0;
359  do
360  {
361  df [coef_index] = sf [coef_index] + sf [coef_index + 4];
362  df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
363  df [coef_index + 2] = sf [coef_index + 2] * 2;
364  df [coef_index + 6] = sf [coef_index + 6] * 2;
365 
366  df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
367  df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
368 
369  const flt_t vr = sf [coef_index + 1] - sf [coef_index + 3];
370  const flt_t vi = sf [coef_index + 5] + sf [coef_index + 7];
371 
372  df [coef_index + 5] = (vr + vi) * sqrt2_2;
373  df [coef_index + 7] = (vi - vr) * sqrt2_2;
374 
375  coef_index += 8;
376  }
377  while (coef_index < _length);
378  }
379 
380  /* Penultimate and last pass at once */
381  {
382  coef_index = 0;
383  const long * bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
384  const flt_t * sf2 = df;
385  do
386  {
387  {
388  const flt_t b_0 = sf2 [0] + sf2 [2];
389  const flt_t b_2 = sf2 [0] - sf2 [2];
390  const flt_t b_1 = sf2 [1] * 2;
391  const flt_t b_3 = sf2 [3] * 2;
392 
393  x [bit_rev_lut_ptr [0]] = b_0 + b_1;
394  x [bit_rev_lut_ptr [1]] = b_0 - b_1;
395  x [bit_rev_lut_ptr [2]] = b_2 + b_3;
396  x [bit_rev_lut_ptr [3]] = b_2 - b_3;
397  }
398  {
399  const flt_t b_0 = sf2 [4] + sf2 [6];
400  const flt_t b_2 = sf2 [4] - sf2 [6];
401  const flt_t b_1 = sf2 [5] * 2;
402  const flt_t b_3 = sf2 [7] * 2;
403 
404  x [bit_rev_lut_ptr [4]] = b_0 + b_1;
405  x [bit_rev_lut_ptr [5]] = b_0 - b_1;
406  x [bit_rev_lut_ptr [6]] = b_2 + b_3;
407  x [bit_rev_lut_ptr [7]] = b_2 - b_3;
408  }
409 
410  sf2 += 8;
411  coef_index += 8;
412  bit_rev_lut_ptr += 8;
413  }
414  while (coef_index < _length);
415  }
416  }
417  }
418 
419 /*______________________________________________
420  *
421  * Special cases
422  *______________________________________________
423  */
424 
425  /* 4-point IFFT */
426  else if (_nbr_bits == 2)
427  {
428  const flt_t b_0 = f [0] + f [2];
429  const flt_t b_2 = f [0] - f [2];
430 
431  x [0] = b_0 + f [1] * 2;
432  x [2] = b_0 - f [1] * 2;
433  x [1] = b_2 + f [3] * 2;
434  x [3] = b_2 - f [3] * 2;
435  }
436 
437  /* 2-point IFFT */
438  else if (_nbr_bits == 1)
439  {
440  x [0] = f [0] + f [1];
441  x [1] = f [0] - f [1];
442  }
443 
444  /* 1-point IFFT */
445  else
446  {
447  x [0] = f [0];
448  }
449 }
450 
451 /*==========================================================================*/
452 /* Name: rescale */
453 /* Description: Scale an array by divide each element by its length. */
454 /* This function should be called after FFT + IFFT. */
455 /* Input/Output parameters: */
456 /* - x: pointer on array to rescale (time or frequency). */
457 /* Throws: Nothing */
458 /*==========================================================================*/
459 
460 void FFTReal::rescale (flt_t x []) const
461 {
462  const flt_t mul = flt_t (1.0 / _length);
463  long i = _length - 1;
464 
465  do
466  {
467  x [i] *= mul;
468  --i;
469  }
470  while (i >= 0);
471 }
472 
473 /*\\\ NESTED CLASS MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
474 
475 /*==========================================================================*/
476 /* Name: Constructor */
477 /* Input parameters: */
478 /* - nbr_bits: number of bits of the array on which we want to do a */
479 /* FFT. Range: > 0 */
480 /* Throws: std::bad_alloc */
481 /*==========================================================================*/
482 
483 FFTReal::BitReversedLUT::BitReversedLUT (const int nbr_bits)
484 {
485  long length;
486  long cnt;
487  long br_index;
488  long bit;
489 
490  length = 1L << nbr_bits;
491  _ptr = new long [length];
492 
493  br_index = 0;
494  _ptr [0] = 0;
495  for (cnt = 1; cnt < length; ++cnt)
496  {
497  /* ++br_index (bit reversed) */
498  bit = length >> 1;
499  while (((br_index ^= bit) & bit) == 0)
500  {
501  bit >>= 1;
502  }
503 
504  _ptr [cnt] = br_index;
505  }
506 }
507 
508 /*==========================================================================*/
509 /* Name: Destructor */
510 /*==========================================================================*/
511 
512 FFTReal::BitReversedLUT::~BitReversedLUT (void)
513 {
514  delete [] _ptr;
515  _ptr = 0;
516 }
517 
518 /*==========================================================================*/
519 /* Name: Constructor */
520 /* Input parameters: */
521 /* - nbr_bits: number of bits of the array on which we want to do a */
522 /* FFT. Range: > 0 */
523 /* Throws: std::bad_alloc, anything */
524 /*==========================================================================*/
525 
526 FFTReal::TrigoLUT::TrigoLUT (const int nbr_bits)
527 {
528  long total_len;
529 
530  _ptr = 0;
531  if (nbr_bits > 3)
532  {
533  total_len = (1L << (nbr_bits - 1)) - 4;
534  _ptr = new flt_t [total_len];
535 
536  const double PI = atanf (1) * 4;
537  for (int level = 3; level < nbr_bits; ++level)
538  {
539  const long level_len = 1L << (level - 1);
540  flt_t * const level_ptr = const_cast<flt_t *> (get_ptr (level));
541  const double mul = PI / (level_len << 1);
542 
543  for (long i = 0; i < level_len; ++ i)
544  {
545  level_ptr [i] = (flt_t) cos (i * mul);
546  }
547  }
548  }
549 }
550 
551 /*==========================================================================*/
552 /* Name: Destructor */
553 /*==========================================================================*/
554 
555 FFTReal::TrigoLUT::~TrigoLUT (void)
556 {
557  delete [] _ptr;
558  _ptr = 0;
559 }
560 
561 #if defined (_MSC_VER)
562 #pragma pack (pop)
563 #endif // _MSC_VER
564 
565 /*****************************************************************************
566 
567  LEGAL
568 
569  Source code may be freely used for any purpose, including commercial
570  applications. Programs must display in their "About" dialog-box (or
571  documentation) a text telling they use these routines by Laurent de Soras.
572  Modified source code can be distributed, but modifications must be clearly
573  indicated.
574 
575  CONTACT
576 
577  Laurent de Soras
578  92 avenue Albert 1er
579  92500 Rueil-Malmaison
580  France
581 
582  ldesoras@club-internet.fr
583 
584 *****************************************************************************/
585 
586 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/