FFTReal.hpp (21034B)
1 /***************************************************************************** 2 3 FFTReal.hpp 4 By Laurent de Soras 5 6 --- Legal stuff --- 7 8 This program is free software. It comes without any warranty, to 9 the extent permitted by applicable law. You can redistribute it 10 and/or modify it under the terms of the Do What The Fuck You Want 11 To Public License, Version 2, as published by Sam Hocevar. See 12 http://sam.zoy.org/wtfpl/COPYING for more details. 13 14 *Tab=3***********************************************************************/ 15 16 17 18 #if defined (ffft_FFTReal_CURRENT_CODEHEADER) 19 #error Recursive inclusion of FFTReal code header. 20 #endif 21 #define ffft_FFTReal_CURRENT_CODEHEADER 22 23 #if ! defined (ffft_FFTReal_CODEHEADER_INCLUDED) 24 #define ffft_FFTReal_CODEHEADER_INCLUDED 25 26 27 28 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ 29 30 #include <cassert> 31 #include <cmath> 32 33 34 35 namespace ffft 36 { 37 38 39 40 static inline bool FFTReal_is_pow2 (long x) 41 { 42 assert (x > 0); 43 44 return ((x & -x) == x); 45 } 46 47 48 49 static inline int FFTReal_get_next_pow2 (long x) 50 { 51 --x; 52 53 int p = 0; 54 while ((x & ~0xFFFFL) != 0) 55 { 56 p += 16; 57 x >>= 16; 58 } 59 while ((x & ~0xFL) != 0) 60 { 61 p += 4; 62 x >>= 4; 63 } 64 while (x > 0) 65 { 66 ++p; 67 x >>= 1; 68 } 69 70 return (p); 71 } 72 73 74 75 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ 76 77 78 79 /* 80 ============================================================================== 81 Name: ctor 82 Input parameters: 83 - length: length of the array on which we want to do a FFT. Range: power of 84 2 only, > 0. 85 Throws: std::bad_alloc 86 ============================================================================== 87 */ 88 89 template <class DT> 90 FFTReal <DT>::FFTReal (long length) 91 : _length (length) 92 , _nbr_bits (FFTReal_get_next_pow2 (length)) 93 , _br_lut () 94 , _trigo_lut () 95 , _buffer (length) 96 , _trigo_osc () 97 { 98 assert (FFTReal_is_pow2 (length)); 99 assert (_nbr_bits <= MAX_BIT_DEPTH); 100 101 init_br_lut (); 102 init_trigo_lut (); 103 init_trigo_osc (); 104 } 105 106 107 108 /* 109 ============================================================================== 110 Name: get_length 111 Description: 112 Returns the number of points processed by this FFT object. 113 Returns: The number of points, power of 2, > 0. 114 Throws: Nothing 115 ============================================================================== 116 */ 117 118 template <class DT> 119 long FFTReal <DT>::get_length () const 120 { 121 return (_length); 122 } 123 124 125 126 /* 127 ============================================================================== 128 Name: do_fft 129 Description: 130 Compute the FFT of the array. 131 Input parameters: 132 - x: pointer on the source array (time). 133 Output parameters: 134 - f: pointer on the destination array (frequencies). 135 f [0...length(x)/2] = real values, 136 f [length(x)/2+1...length(x)-1] = negative imaginary values of 137 coefficents 1...length(x)/2-1. 138 Throws: Nothing 139 ============================================================================== 140 */ 141 142 template <class DT> 143 void FFTReal <DT>::do_fft (DataType f [], const DataType x []) const 144 { 145 assert (f != 0); 146 assert (f != use_buffer ()); 147 assert (x != 0); 148 assert (x != use_buffer ()); 149 assert (x != f); 150 151 // General case 152 if (_nbr_bits > 2) 153 { 154 compute_fft_general (f, x); 155 } 156 157 // 4-point FFT 158 else if (_nbr_bits == 2) 159 { 160 f [1] = x [0] - x [2]; 161 f [3] = x [1] - x [3]; 162 163 const DataType b_0 = x [0] + x [2]; 164 const DataType b_2 = x [1] + x [3]; 165 166 f [0] = b_0 + b_2; 167 f [2] = b_0 - b_2; 168 } 169 170 // 2-point FFT 171 else if (_nbr_bits == 1) 172 { 173 f [0] = x [0] + x [1]; 174 f [1] = x [0] - x [1]; 175 } 176 177 // 1-point FFT 178 else 179 { 180 f [0] = x [0]; 181 } 182 } 183 184 185 186 /* 187 ============================================================================== 188 Name: do_ifft 189 Description: 190 Compute the inverse FFT of the array. Note that data must be post-scaled: 191 IFFT (FFT (x)) = x * length (x). 192 Input parameters: 193 - f: pointer on the source array (frequencies). 194 f [0...length(x)/2] = real values 195 f [length(x)/2+1...length(x)-1] = negative imaginary values of 196 coefficents 1...length(x)/2-1. 197 Output parameters: 198 - x: pointer on the destination array (time). 199 Throws: Nothing 200 ============================================================================== 201 */ 202 203 template <class DT> 204 void FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const 205 { 206 assert (f != 0); 207 assert (f != use_buffer ()); 208 assert (x != 0); 209 assert (x != use_buffer ()); 210 assert (x != f); 211 212 // General case 213 if (_nbr_bits > 2) 214 { 215 compute_ifft_general (f, x); 216 } 217 218 // 4-point IFFT 219 else if (_nbr_bits == 2) 220 { 221 const DataType b_0 = f [0] + f [2]; 222 const DataType b_2 = f [0] - f [2]; 223 224 x [0] = b_0 + f [1] * 2; 225 x [2] = b_0 - f [1] * 2; 226 x [1] = b_2 + f [3] * 2; 227 x [3] = b_2 - f [3] * 2; 228 } 229 230 // 2-point IFFT 231 else if (_nbr_bits == 1) 232 { 233 x [0] = f [0] + f [1]; 234 x [1] = f [0] - f [1]; 235 } 236 237 // 1-point IFFT 238 else 239 { 240 x [0] = f [0]; 241 } 242 } 243 244 245 246 /* 247 ============================================================================== 248 Name: rescale 249 Description: 250 Scale an array by divide each element by its length. This function should 251 be called after FFT + IFFT. 252 Input parameters: 253 - x: pointer on array to rescale (time or frequency). 254 Throws: Nothing 255 ============================================================================== 256 */ 257 258 template <class DT> 259 void FFTReal <DT>::rescale (DataType x []) const 260 { 261 const DataType mul = DataType (1.0 / _length); 262 263 if (_length < 4) 264 { 265 long i = _length - 1; 266 do 267 { 268 x [i] *= mul; 269 --i; 270 } 271 while (i >= 0); 272 } 273 274 else 275 { 276 assert ((_length & 3) == 0); 277 278 // Could be optimized with SIMD instruction sets (needs alignment check) 279 long i = _length - 4; 280 do 281 { 282 x [i + 0] *= mul; 283 x [i + 1] *= mul; 284 x [i + 2] *= mul; 285 x [i + 3] *= mul; 286 i -= 4; 287 } 288 while (i >= 0); 289 } 290 } 291 292 293 294 /* 295 ============================================================================== 296 Name: use_buffer 297 Description: 298 Access the internal buffer, whose length is the FFT one. 299 Buffer content will be erased at each do_fft() / do_ifft() call! 300 This buffer cannot be used as: 301 - source for FFT or IFFT done with this object 302 - destination for FFT or IFFT done with this object 303 Returns: 304 Buffer start address 305 Throws: Nothing 306 ============================================================================== 307 */ 308 309 template <class DT> 310 typename FFTReal <DT>::DataType * FFTReal <DT>::use_buffer () const 311 { 312 return (&_buffer [0]); 313 } 314 315 316 317 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ 318 319 320 321 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ 322 323 324 325 template <class DT> 326 void FFTReal <DT>::init_br_lut () 327 { 328 const long length = 1L << _nbr_bits; 329 _br_lut.resize (length); 330 331 _br_lut [0] = 0; 332 long br_index = 0; 333 for (long cnt = 1; cnt < length; ++cnt) 334 { 335 // ++br_index (bit reversed) 336 long bit = length >> 1; 337 while (((br_index ^= bit) & bit) == 0) 338 { 339 bit >>= 1; 340 } 341 342 _br_lut [cnt] = br_index; 343 } 344 } 345 346 347 348 template <class DT> 349 void FFTReal <DT>::init_trigo_lut () 350 { 351 using namespace std; 352 353 if (_nbr_bits > 3) 354 { 355 const long total_len = (1L << (_nbr_bits - 1)) - 4; 356 _trigo_lut.resize (total_len); 357 358 for (int level = 3; level < _nbr_bits; ++level) 359 { 360 const long level_len = 1L << (level - 1); 361 DataType * const level_ptr = 362 &_trigo_lut [get_trigo_level_index (level)]; 363 const double mul = PI / (level_len << 1); 364 365 for (long i = 0; i < level_len; ++ i) 366 { 367 level_ptr [i] = static_cast <DataType> (cos (i * mul)); 368 } 369 } 370 } 371 } 372 373 374 375 template <class DT> 376 void FFTReal <DT>::init_trigo_osc () 377 { 378 const int nbr_osc = _nbr_bits - TRIGO_BD_LIMIT; 379 if (nbr_osc > 0) 380 { 381 _trigo_osc.resize (nbr_osc); 382 383 for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt) 384 { 385 OscType & osc = _trigo_osc [osc_cnt]; 386 387 const long len = 1L << (TRIGO_BD_LIMIT + osc_cnt); 388 const double mul = (0.5 * PI) / len; 389 osc.set_step (mul); 390 } 391 } 392 } 393 394 395 396 template <class DT> 397 const long * FFTReal <DT>::get_br_ptr () const 398 { 399 return (&_br_lut [0]); 400 } 401 402 403 404 template <class DT> 405 const typename FFTReal <DT>::DataType * FFTReal <DT>::get_trigo_ptr (int level) const 406 { 407 assert (level >= 3); 408 409 return (&_trigo_lut [get_trigo_level_index (level)]); 410 } 411 412 413 414 template <class DT> 415 long FFTReal <DT>::get_trigo_level_index (int level) const 416 { 417 assert (level >= 3); 418 419 return ((1L << (level - 1)) - 4); 420 } 421 422 423 424 // Transform in several passes 425 template <class DT> 426 void FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const 427 { 428 assert (f != 0); 429 assert (f != use_buffer ()); 430 assert (x != 0); 431 assert (x != use_buffer ()); 432 assert (x != f); 433 434 DataType * sf; 435 DataType * df; 436 437 if ((_nbr_bits & 1) != 0) 438 { 439 df = use_buffer (); 440 sf = f; 441 } 442 else 443 { 444 df = f; 445 sf = use_buffer (); 446 } 447 448 compute_direct_pass_1_2 (df, x); 449 compute_direct_pass_3 (sf, df); 450 451 for (int pass = 3; pass < _nbr_bits; ++ pass) 452 { 453 compute_direct_pass_n (df, sf, pass); 454 455 DataType * const temp_ptr = df; 456 df = sf; 457 sf = temp_ptr; 458 } 459 } 460 461 462 463 template <class DT> 464 void FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const 465 { 466 assert (df != 0); 467 assert (x != 0); 468 assert (df != x); 469 470 const long * const bit_rev_lut_ptr = get_br_ptr (); 471 long coef_index = 0; 472 do 473 { 474 const long rev_index_0 = bit_rev_lut_ptr [coef_index]; 475 const long rev_index_1 = bit_rev_lut_ptr [coef_index + 1]; 476 const long rev_index_2 = bit_rev_lut_ptr [coef_index + 2]; 477 const long rev_index_3 = bit_rev_lut_ptr [coef_index + 3]; 478 479 DataType * const df2 = df + coef_index; 480 df2 [1] = x [rev_index_0] - x [rev_index_1]; 481 df2 [3] = x [rev_index_2] - x [rev_index_3]; 482 483 const DataType sf_0 = x [rev_index_0] + x [rev_index_1]; 484 const DataType sf_2 = x [rev_index_2] + x [rev_index_3]; 485 486 df2 [0] = sf_0 + sf_2; 487 df2 [2] = sf_0 - sf_2; 488 489 coef_index += 4; 490 } 491 while (coef_index < _length); 492 } 493 494 495 496 template <class DT> 497 void FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const 498 { 499 assert (df != 0); 500 assert (sf != 0); 501 assert (df != sf); 502 503 const DataType sqrt2_2 = DataType (SQRT2 * 0.5); 504 long coef_index = 0; 505 do 506 { 507 DataType v; 508 509 df [coef_index] = sf [coef_index] + sf [coef_index + 4]; 510 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4]; 511 df [coef_index + 2] = sf [coef_index + 2]; 512 df [coef_index + 6] = sf [coef_index + 6]; 513 514 v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2; 515 df [coef_index + 1] = sf [coef_index + 1] + v; 516 df [coef_index + 3] = sf [coef_index + 1] - v; 517 518 v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2; 519 df [coef_index + 5] = v + sf [coef_index + 3]; 520 df [coef_index + 7] = v - sf [coef_index + 3]; 521 522 coef_index += 8; 523 } 524 while (coef_index < _length); 525 } 526 527 528 529 template <class DT> 530 void FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const 531 { 532 assert (df != 0); 533 assert (sf != 0); 534 assert (df != sf); 535 assert (pass >= 3); 536 assert (pass < _nbr_bits); 537 538 if (pass <= TRIGO_BD_LIMIT) 539 { 540 compute_direct_pass_n_lut (df, sf, pass); 541 } 542 else 543 { 544 compute_direct_pass_n_osc (df, sf, pass); 545 } 546 } 547 548 549 550 template <class DT> 551 void FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const 552 { 553 assert (df != 0); 554 assert (sf != 0); 555 assert (df != sf); 556 assert (pass >= 3); 557 assert (pass < _nbr_bits); 558 559 const long nbr_coef = 1 << pass; 560 const long h_nbr_coef = nbr_coef >> 1; 561 const long d_nbr_coef = nbr_coef << 1; 562 long coef_index = 0; 563 const DataType * const cos_ptr = get_trigo_ptr (pass); 564 do 565 { 566 const DataType * const sf1r = sf + coef_index; 567 const DataType * const sf2r = sf1r + nbr_coef; 568 DataType * const dfr = df + coef_index; 569 DataType * const dfi = dfr + nbr_coef; 570 571 // Extreme coefficients are always real 572 dfr [0] = sf1r [0] + sf2r [0]; 573 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] = 574 dfr [h_nbr_coef] = sf1r [h_nbr_coef]; 575 dfi [h_nbr_coef] = sf2r [h_nbr_coef]; 576 577 // Others are conjugate complex numbers 578 const DataType * const sf1i = sf1r + h_nbr_coef; 579 const DataType * const sf2i = sf1i + nbr_coef; 580 for (long i = 1; i < h_nbr_coef; ++ i) 581 { 582 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef); 583 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef); 584 DataType v; 585 586 v = sf2r [i] * c - sf2i [i] * s; 587 dfr [i] = sf1r [i] + v; 588 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] = 589 590 v = sf2r [i] * s + sf2i [i] * c; 591 dfi [i] = v + sf1i [i]; 592 dfi [nbr_coef - i] = v - sf1i [i]; 593 } 594 595 coef_index += d_nbr_coef; 596 } 597 while (coef_index < _length); 598 } 599 600 601 602 template <class DT> 603 void FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const 604 { 605 assert (df != 0); 606 assert (sf != 0); 607 assert (df != sf); 608 assert (pass > TRIGO_BD_LIMIT); 609 assert (pass < _nbr_bits); 610 611 const long nbr_coef = 1 << pass; 612 const long h_nbr_coef = nbr_coef >> 1; 613 const long d_nbr_coef = nbr_coef << 1; 614 long coef_index = 0; 615 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)]; 616 do 617 { 618 const DataType * const sf1r = sf + coef_index; 619 const DataType * const sf2r = sf1r + nbr_coef; 620 DataType * const dfr = df + coef_index; 621 DataType * const dfi = dfr + nbr_coef; 622 623 osc.clear_buffers (); 624 625 // Extreme coefficients are always real 626 dfr [0] = sf1r [0] + sf2r [0]; 627 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] = 628 dfr [h_nbr_coef] = sf1r [h_nbr_coef]; 629 dfi [h_nbr_coef] = sf2r [h_nbr_coef]; 630 631 // Others are conjugate complex numbers 632 const DataType * const sf1i = sf1r + h_nbr_coef; 633 const DataType * const sf2i = sf1i + nbr_coef; 634 for (long i = 1; i < h_nbr_coef; ++ i) 635 { 636 osc.step (); 637 const DataType c = osc.get_cos (); 638 const DataType s = osc.get_sin (); 639 DataType v; 640 641 v = sf2r [i] * c - sf2i [i] * s; 642 dfr [i] = sf1r [i] + v; 643 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] = 644 645 v = sf2r [i] * s + sf2i [i] * c; 646 dfi [i] = v + sf1i [i]; 647 dfi [nbr_coef - i] = v - sf1i [i]; 648 } 649 650 coef_index += d_nbr_coef; 651 } 652 while (coef_index < _length); 653 } 654 655 656 657 // Transform in several pass 658 template <class DT> 659 void FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const 660 { 661 assert (f != 0); 662 assert (f != use_buffer ()); 663 assert (x != 0); 664 assert (x != use_buffer ()); 665 assert (x != f); 666 667 DataType * sf = const_cast <DataType *> (f); 668 DataType * df; 669 DataType * df_temp; 670 671 if (_nbr_bits & 1) 672 { 673 df = use_buffer (); 674 df_temp = x; 675 } 676 else 677 { 678 df = x; 679 df_temp = use_buffer (); 680 } 681 682 for (int pass = _nbr_bits - 1; pass >= 3; -- pass) 683 { 684 compute_inverse_pass_n (df, sf, pass); 685 686 if (pass < _nbr_bits - 1) 687 { 688 DataType * const temp_ptr = df; 689 df = sf; 690 sf = temp_ptr; 691 } 692 else 693 { 694 sf = df; 695 df = df_temp; 696 } 697 } 698 699 compute_inverse_pass_3 (df, sf); 700 compute_inverse_pass_1_2 (x, df); 701 } 702 703 704 705 template <class DT> 706 void FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const 707 { 708 assert (df != 0); 709 assert (sf != 0); 710 assert (df != sf); 711 assert (pass >= 3); 712 assert (pass < _nbr_bits); 713 714 if (pass <= TRIGO_BD_LIMIT) 715 { 716 compute_inverse_pass_n_lut (df, sf, pass); 717 } 718 else 719 { 720 compute_inverse_pass_n_osc (df, sf, pass); 721 } 722 } 723 724 725 726 template <class DT> 727 void FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const 728 { 729 assert (df != 0); 730 assert (sf != 0); 731 assert (df != sf); 732 assert (pass >= 3); 733 assert (pass < _nbr_bits); 734 735 const long nbr_coef = 1 << pass; 736 const long h_nbr_coef = nbr_coef >> 1; 737 const long d_nbr_coef = nbr_coef << 1; 738 long coef_index = 0; 739 const DataType * const cos_ptr = get_trigo_ptr (pass); 740 do 741 { 742 const DataType * const sfr = sf + coef_index; 743 const DataType * const sfi = sfr + nbr_coef; 744 DataType * const df1r = df + coef_index; 745 DataType * const df2r = df1r + nbr_coef; 746 747 // Extreme coefficients are always real 748 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef] 749 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef] 750 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2; 751 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2; 752 753 // Others are conjugate complex numbers 754 DataType * const df1i = df1r + h_nbr_coef; 755 DataType * const df2i = df1i + nbr_coef; 756 for (long i = 1; i < h_nbr_coef; ++ i) 757 { 758 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i] 759 df1i [i] = sfi [i] - sfi [nbr_coef - i]; 760 761 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef); 762 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef); 763 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i] 764 const DataType vi = sfi [i] + sfi [nbr_coef - i]; 765 766 df2r [i] = vr * c + vi * s; 767 df2i [i] = vi * c - vr * s; 768 } 769 770 coef_index += d_nbr_coef; 771 } 772 while (coef_index < _length); 773 } 774 775 776 777 template <class DT> 778 void FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const 779 { 780 assert (df != 0); 781 assert (sf != 0); 782 assert (df != sf); 783 assert (pass > TRIGO_BD_LIMIT); 784 assert (pass < _nbr_bits); 785 786 const long nbr_coef = 1 << pass; 787 const long h_nbr_coef = nbr_coef >> 1; 788 const long d_nbr_coef = nbr_coef << 1; 789 long coef_index = 0; 790 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)]; 791 do 792 { 793 const DataType * const sfr = sf + coef_index; 794 const DataType * const sfi = sfr + nbr_coef; 795 DataType * const df1r = df + coef_index; 796 DataType * const df2r = df1r + nbr_coef; 797 798 osc.clear_buffers (); 799 800 // Extreme coefficients are always real 801 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef] 802 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef] 803 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2; 804 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2; 805 806 // Others are conjugate complex numbers 807 DataType * const df1i = df1r + h_nbr_coef; 808 DataType * const df2i = df1i + nbr_coef; 809 for (long i = 1; i < h_nbr_coef; ++ i) 810 { 811 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i] 812 df1i [i] = sfi [i] - sfi [nbr_coef - i]; 813 814 osc.step (); 815 const DataType c = osc.get_cos (); 816 const DataType s = osc.get_sin (); 817 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i] 818 const DataType vi = sfi [i] + sfi [nbr_coef - i]; 819 820 df2r [i] = vr * c + vi * s; 821 df2i [i] = vi * c - vr * s; 822 } 823 824 coef_index += d_nbr_coef; 825 } 826 while (coef_index < _length); 827 } 828 829 830 831 template <class DT> 832 void FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const 833 { 834 assert (df != 0); 835 assert (sf != 0); 836 assert (df != sf); 837 838 const DataType sqrt2_2 = DataType (SQRT2 * 0.5); 839 long coef_index = 0; 840 do 841 { 842 df [coef_index] = sf [coef_index] + sf [coef_index + 4]; 843 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4]; 844 df [coef_index + 2] = sf [coef_index + 2] * 2; 845 df [coef_index + 6] = sf [coef_index + 6] * 2; 846 847 df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3]; 848 df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7]; 849 850 const DataType vr = sf [coef_index + 1] - sf [coef_index + 3]; 851 const DataType vi = sf [coef_index + 5] + sf [coef_index + 7]; 852 853 df [coef_index + 5] = (vr + vi) * sqrt2_2; 854 df [coef_index + 7] = (vi - vr) * sqrt2_2; 855 856 coef_index += 8; 857 } 858 while (coef_index < _length); 859 } 860 861 862 863 template <class DT> 864 void FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const 865 { 866 assert (x != 0); 867 assert (sf != 0); 868 assert (x != sf); 869 870 const long * bit_rev_lut_ptr = get_br_ptr (); 871 const DataType * sf2 = sf; 872 long coef_index = 0; 873 do 874 { 875 { 876 const DataType b_0 = sf2 [0] + sf2 [2]; 877 const DataType b_2 = sf2 [0] - sf2 [2]; 878 const DataType b_1 = sf2 [1] * 2; 879 const DataType b_3 = sf2 [3] * 2; 880 881 x [bit_rev_lut_ptr [0]] = b_0 + b_1; 882 x [bit_rev_lut_ptr [1]] = b_0 - b_1; 883 x [bit_rev_lut_ptr [2]] = b_2 + b_3; 884 x [bit_rev_lut_ptr [3]] = b_2 - b_3; 885 } 886 { 887 const DataType b_0 = sf2 [4] + sf2 [6]; 888 const DataType b_2 = sf2 [4] - sf2 [6]; 889 const DataType b_1 = sf2 [5] * 2; 890 const DataType b_3 = sf2 [7] * 2; 891 892 x [bit_rev_lut_ptr [4]] = b_0 + b_1; 893 x [bit_rev_lut_ptr [5]] = b_0 - b_1; 894 x [bit_rev_lut_ptr [6]] = b_2 + b_3; 895 x [bit_rev_lut_ptr [7]] = b_2 - b_3; 896 } 897 898 sf2 += 8; 899 coef_index += 8; 900 bit_rev_lut_ptr += 8; 901 } 902 while (coef_index < _length); 903 } 904 905 906 907 } // namespace ffft 908 909 910 911 #endif // ffft_FFTReal_CODEHEADER_INCLUDED 912 913 #undef ffft_FFTReal_CURRENT_CODEHEADER 914 915 916 917 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/