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