1/*=============================================================================
2|
3| NAME
4|
5| FFT.cpp
6|
7| DESCRIPTION
8|
9| This program does an in-place one dimensional complex discrete
10| Fourier transform. It uses an FFT algorithm for a number of
11| input points which is a power of two. It's implemented in C++
12| as a derived class of the STL vector type.
13|
14| AUTHOR
15|
16| Sean O'Connor 19 November 2005
17|
18| LEGAL
19|
20| FFT Version 1.6 - An FFT utility library in C++.
21| Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved.
22|
23| This program is free software: you can redistribute it and/or modify
24| it under the terms of the GNU General Public License as published by
25| the Free Software Foundation, either version 3 of the License, or
26| (at your option) any later version.
27|
28| This program is distributed in the hope that it will be useful,
29| but WITHOUT ANY WARRANTY; without even the implied warranty of
30| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31| GNU General Public License for more details.
32|
33| You should have received a copy of the GNU General Public License
34| along with this program. If not, see <http://www.gnu.org/licenses/>.
35|
36| The author's address is seanerikoconnor!AT!gmail!DOT!com
37| with !DOT! replaced by . and the !AT! replaced by @
38|
39+============================================================================*/
40
41/*------------------------------------------------------------------------------
42| Include Files |
43------------------------------------------------------------------------------*/
44
45#include <iostream> // Basic stream I/O.
46#include <fstream> // File stream I/O.
47#include <sstream> // String stream I/O.
48#include <cmath> // Basic math functions.
49#include <complex> // Complex data type and operations.
50#include <vector> // STL vector class.
51#include <algorithm> // Iterators.
52#include <string> // STL string class.
53#include <stdexcept> // Exceptions.
54#include <numbers> // Constants like pi_v<float type>
55
56using namespace std ; // Let us omit the std:: prefix for every STL class.
57
58#include "FFT.hpp" // Our own FFT stuff.
59
60
61/*=============================================================================
62|
63| NAME
64|
65| FastFourierTransform
66|
67| DESCRIPTION
68|
69| Default constructor. Just call the base class constructor for
70| the vector type, then initialize the derived class field.
71|
72+============================================================================*/
73
74template <typename FloatType>
75FastFourierTransform<FloatType>::FastFourierTransform()
76 : vector< complex<FloatType> >()
77 , direction( false )
78{
79 // Get the sine and cosine tables initialized.
80 /// MAX_TABLE_SIZE ) ;
81 starting_multiplier.clear();
82}
83
84
85/*=============================================================================
86|
87| NAME
88|
89| ~FastFourierTransform
90|
91| DESCRIPTION
92|
93| Destructor. Doesn't need to do much.
94|
95+============================================================================*/
96
97template <typename FloatType>
98FastFourierTransform<FloatType>::~FastFourierTransform()
99{
100}
101
102
103
104/*=============================================================================
105|
106| NAME
107|
108| FastFourierTransform
109|
110| DESCRIPTION
111|
112| Copy constructor.
113|
114+============================================================================*/
115
116template <typename FloatType>
117FastFourierTransform<FloatType>::FastFourierTransform( const FastFourierTransform & transform )
118 : vector< complex<FloatType> >( transform ) // Initialize base class first.
119{
120 // Copy derived class stuff.
121 direction = transform.direction ;
122}
123
124
125/*=============================================================================
126|
127| NAME
128|
129| =
130|
131| DESCRIPTION
132|
133| Assignment operator.
134|
135+============================================================================*/
136
137template <typename FloatType>
138FastFourierTransform<FloatType> & FastFourierTransform<FloatType>::operator=( const FastFourierTransform & transform )
139{
140 // Check for initializing to oneself: pass back a reference to the unchanged object.
141 if (this == &transform)
142 {
143 // Pass back the object so we can concatenate assignments.
144 return *this ;
145 }
146
147 // Call base class assignment operator first.
148 FastFourierTransform<FloatType>::operator=( transform ) ;
149
150 // Assignment for derived class fields.
151 direction = transform.direction ;
152
153
154 // Pass back the object so we can concatenate assignments.
155 return *this ;
156}
157
158
159/*=============================================================================
160|
161| NAME
162|
163| fft
164|
165| DESCRIPTION
166|
167|
168| This member function does an in-place one dimensional complex
169| discrete Fourier transform. It uses an FFT algorithm for a number
170| of input points which is a power of two.
171|
172| PARAMETERS
173|
174| direction The direction of the FFT, true for a forward transform and
175| false for an inverse transform.
176|
177| EXCEPTIONS
178|
179| Throws FastFourierTransform type exceptions for the cases where
180|
181| -Number of points in the transform is less than 2.
182| -Number of points is greater than the sine/cosine table can handle.
183| -Number of points is not a power of 2.
184|
185| EXAMPLE
186|
187| Suppose the vector contains the data
188|
189| X ( 0 ) = 1 + i
190| X ( 1 ) = 2 + 2 i
191| X ( 2 ) = 3 + 3 i
192| X ( 3 ) = 4 + 4 i
193|
194| then the forward FFT returns
195|
196| X ( 0 ) = 5 + 5 i
197| X ( 1 ) = -2
198| X ( 2 ) = -1 - i
199| X ( 3 ) = - 2 i
200|
201| and the inverse FFT returns the original data to within roundoff.
202|
203| ^
204| The Fourier coefficients X( k ) are defined for k = 0, . . ., N - 1
205| by the formula
206|
207| N - 1 2 Pi i
208| ^ 1 ____ - ------ j k
209| X( k ) = --------- \ N
210| ----- / X( j ) e
211| \ / N ----
212| v j = 0
213|
214| where e is 2.1718 . . ., Pi = 3.14159. . . and i is the square root
215| of -1.
216|
217| This is the formula for the forward transform. The inverse transform
218| is the same, except that the minus sign in the exponent above is a
219| plus sign instead.
220| ---
221| The normalizing factor 1 / \/ N in front was picked so that the
222| forward and inverse transforms look the same.
223|
224| The particular type of FFT algorithm we use is a
225|
226| Complex,
227| Two to the N,
228| Decimation in frequency,
229| Input in order, output bit reversed
230|
231| type of FFT. It is described in the book
232|
233| THE FAST FOURIER TRANSFORM,
234| F. Oran Brigham, Prentice Hall, 1974
235|
236| Also see my notes for an in-depth explanation.
237|
238| BUGS
239|
240| The bit-reversal routine could be made faster by using bit shifting
241| operations.
242|
243| You might want to use long double precision complex numbers to get more
244| accuracy.
245|
246| Check for a power of 2 probably ought to use the machine's floating point
247| epsilon (computed automatically).
248|
249+============================================================================*/
250
251template <typename FloatType>
252void FastFourierTransform<FloatType>::fft( bool direction )
253{
254// Number of complex points in the FFT from the vector's size.
255unsigned int N = this->size() ;
256
257/*
258 * Compute log base 2 of N. Perturb slightly to avoid roundoff error.
259 * For example, (int)( 4.99999999999 ) = 4 (wrong!), but
260 * (int)( 4.999999999 + .01 ) = 5 (correct!)
261 */
262int log_of_n = (int) ( log( (FloatType) N ) / ( log( 2.0 ) ) + .01 ) ;
263
264
265/*
266 * Error checking. Print a message to the terminal if any one of these
267 * conditions is true:
268 *
269 * N < 2.
270 * log2( N ) > MAX_TABLE_SIZE + 1.
271 * N is not a power of 2.
272 */
273if (N < 2)
274{
275 ostringstream os ;
276 os << "FFT has less than two points. N = " << N << " file: " << __FILE__ << " line: " << __LINE__ ;;
277 throw FastFourierTransformException( os.str() ) ;
278}
279else if ( (log_of_n - 1) > MAX_TABLE_SIZE)
280{
281 ostringstream os ;
282 os << "FFT has too many points for its sin/cos table. Table size = " << MAX_TABLE_SIZE << " file: " << __FILE__ << " line: " << __LINE__ ;
283 throw FastFourierTransformException( os.str() ) ;
284
285}
286else if ( (int) pow( 2.0, log_of_n ) != N )
287{
288 ostringstream os ;
289 os << "Number of points in the FFT is not a power of 2. N = " << N << " file: " << __FILE__ << " line: " << __LINE__ ;;
290 throw FastFourierTransformException( os.str() ) ;
291}
292
293
294/*
295 * If called_already is false, generate the sine and cosine table for
296 * the first time, then set called_already to true. Now later calls to
297 * this function don't recreate the sine and cosine table.
298 */
299if (!called_already)
300{
301 for (int i = 0 ; i <= MAX_TABLE_SIZE - 1 ; ++i)
302 {
303 // Current angle in sine and cosine table.
304 FloatType angle = -PI / pow( (FloatType) 2, (FloatType) i ) ;
305
306 // e^t = cos t + i sin t
307 starting_multiplier[ i ] = complex<FloatType>(cos(angle), sin(angle)) ;
308 }
309
310 called_already = true ;
311}
312
313
314
315/*
316 * Compute the FFT according to its signal flow graph with
317 * computations arranged into butterflies, groups of butterflies
318 * with the same multiplier, and columns.
319 */
320int butterfly_height = N / 2 ;
321int number_of_groups = 1 ;
322
323
324/*
325 * A column refers to a column of the FFT signal flow graph. A group is
326 * a collection of butterflies which have the same multiplier W^p. A
327 * butterfly is the basic computational unit of the FFT.
328 */
329for (int column = 1 ; column <= log_of_n ; ++column)
330{
331 // The multiplier W^p for a group of butterflies.
332 complex<FloatType> multiplier = starting_multiplier[ column - 1 ] ;
333
334 // Modify starting multiplier W^p when computing the inverse transform.
335 if (direction == false)
336 multiplier = conj( multiplier ) ;
337
338 // Difference between current multiplier and next.
339 complex<FloatType> increment = multiplier ;
340
341 // Compute the first group.
342 for (int i = 0 ; i <= butterfly_height - 1 ; ++i)
343 {
344 // Lower branch of a butterfly.
345 int offset = i + butterfly_height ;
346
347 complex<FloatType> temp = (*this)[ i ] + (*this)[ offset ] ;
348 (*this)[ offset ] = (*this)[ i ] - (*this)[ offset ] ;
349 (*this)[ i ] = temp ;
350 }
351
352 // Now do the remaining groups.
353 for (int group = 2 ; group <= number_of_groups ; ++group)
354 {
355 // Array index of first butterfly in group.
356 int start_of_group = reverse_bits( group - 1, log_of_n ) ;
357
358 // Array index of last butterfly in group.
359 int end_of_group = start_of_group + butterfly_height - 1 ;
360
361 // Compute all butterflies in a group.
362 for (int i = start_of_group ; i <= end_of_group ; ++i)
363 {
364 int offset = i + butterfly_height ;
365
366 // Temporary storage in a butterfly calculation.
367 complex<FloatType> product = (*this)[ offset ] * multiplier ;
368
369 complex<FloatType> temp = (*this)[ i ] + product ;
370 (*this)[ offset ] = (*this)[ i ] - product ;
371 (*this)[ i ] = temp ;
372
373 }
374
375 multiplier *= increment ;
376 } // end group
377
378 butterfly_height = butterfly_height / 2 ;
379 number_of_groups = number_of_groups * 2 ;
380
381} // end column
382
383
384
385// Normalize the results by \/ N and permute them into order, two at a time.
386FloatType normalizer = (FloatType) 1.0 / sqrt( (FloatType) N ) ;
387
388
389for (int i = 0 ; i <= N - 1 ; ++i)
390{
391 // Bit-reversed array index i.
392 int rev_i = reverse_bits( i, log_of_n ) ;
393
394 if ( rev_i > i ) // Exchange is not yet done.
395 {
396 complex<FloatType> temp = (*this)[ i ] ;
397 (*this)[ i ] = (*this)[ rev_i ] * normalizer ;
398 (*this)[ rev_i ] = temp * normalizer ;
399 }
400 else if ( rev_i == i ) // exchange with itself
401 {
402 (*this)[ i ] = (*this)[ i ] * normalizer ;
403 }
404}
405
406return ;
407
408} /* ====================== end of function fft ============================= */
409
410
411
412/*=============================================================================
413|
414| NAME
415|
416| reverse_bits
417|
418| DESCRIPTION
419|
420| This function reverses the order of the bits in a number.
421|
422| INPUT
423|
424| n The number whose bits are to be reversed.
425| num_bits The number of bits in N, counting leading zeros.
426|
427| OUTPUT
428|
429| The number N with all its num_bits bits reversed.
430|
431| EXAMPLE CALLING SEQUENCE
432|
433| reverse_bits( (10) , 2 ) = (01) = 1, but
434| 2 2
435|
436| reverse_bits( (10) , 3 ) = reverse_bits( (010) , 3 ) = (010) = (10)
437| 2 2 2 2
438| = 10
439|
440|
441| METHOD
442|
443| The algorithm works by a method similar to base conversion.
444| As the low order bits are broken off, one by one, they are
445| multiplied by 2 and accumulated to generate the bit-reversed
446| number.
447|
448| BUGS
449|
450| The bit-reversal routine could be made faster by using C shift operators.
451|
452+============================================================================*/
453
454template <typename FloatType>
455int FastFourierTransform<FloatType>::reverse_bits( int n, int num_bits )
456{
457 int
458 remainder = n, // Bits remaining to be reversed.
459 reversed_n = 0 ; // Bit-reversed value of n.
460
461 for (int bit_num = 1 ; bit_num <= num_bits ; ++bit_num)
462 {
463 int bit = remainder % 2 ; /* Next least significant bit. */
464 reversed_n = bit + 2 * reversed_n ; /* Shift left and add to buffer. */
465 remainder = (int)( remainder / 2 ) ; /* Remaining bits. */
466 }
467
468 return( reversed_n ) ;
469
470} /* ===================== end of function reverse_bits ===================== */
471
472
473/*==============================================================================
474| Forced Template Instantiations |
475==============================================================================*/
476
477// C++ doesn't automatically generate templated classes or functions until they
478// are first used. So depending on the order of compilation, the linker may say
479// the templated functions are undefined.
480//
481// Therefore, explicitly instantiate ALL templates here:
482
483template FastFourierTransform<float>::FastFourierTransform() ;
484template FastFourierTransform<float>::FastFourierTransform( const FastFourierTransform<float> & ) ;
485template FastFourierTransform<float>::~FastFourierTransform() ;
486template void FastFourierTransform<float>::fft( bool ) ;
487
488template FastFourierTransform<double>::FastFourierTransform() ;
489template FastFourierTransform<double>::FastFourierTransform( const FastFourierTransform<double> & ) ;
490template FastFourierTransform<double>::~FastFourierTransform() ;
491template void FastFourierTransform<double>::fft( bool ) ;
