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