FFTRealPassInverse.hpp 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230
  1. /*****************************************************************************
  2. FFTRealPassInverse.hpp
  3. By Laurent de Soras
  4. --- Legal stuff ---
  5. This program is free software. It comes without any warranty, to
  6. the extent permitted by applicable law. You can redistribute it
  7. and/or modify it under the terms of the Do What The Fuck You Want
  8. To Public License, Version 2, as published by Sam Hocevar. See
  9. http://sam.zoy.org/wtfpl/COPYING for more details.
  10. *Tab=3***********************************************************************/
  11. #if defined (ffft_FFTRealPassInverse_CURRENT_CODEHEADER)
  12. #error Recursive inclusion of FFTRealPassInverse code header.
  13. #endif
  14. #define ffft_FFTRealPassInverse_CURRENT_CODEHEADER
  15. #if ! defined (ffft_FFTRealPassInverse_CODEHEADER_INCLUDED)
  16. #define ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
  17. /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
  18. #include "ffft/FFTRealUseTrigo.h"
  19. namespace ffft
  20. {
  21. /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
  22. template <int PASS>
  23. void FFTRealPassInverse <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType f_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  24. {
  25. process_internal (
  26. len,
  27. dest_ptr,
  28. f_ptr,
  29. cos_ptr,
  30. cos_len,
  31. br_ptr,
  32. osc_list
  33. );
  34. FFTRealPassInverse <PASS - 1>::process_rec (
  35. len,
  36. src_ptr,
  37. dest_ptr,
  38. cos_ptr,
  39. cos_len,
  40. br_ptr,
  41. osc_list
  42. );
  43. }
  44. template <int PASS>
  45. void FFTRealPassInverse <PASS>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  46. {
  47. process_internal (
  48. len,
  49. dest_ptr,
  50. src_ptr,
  51. cos_ptr,
  52. cos_len,
  53. br_ptr,
  54. osc_list
  55. );
  56. FFTRealPassInverse <PASS - 1>::process_rec (
  57. len,
  58. src_ptr,
  59. dest_ptr,
  60. cos_ptr,
  61. cos_len,
  62. br_ptr,
  63. osc_list
  64. );
  65. }
  66. template <>
  67. inline void FFTRealPassInverse <0>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  68. {
  69. // Stops recursion
  70. }
  71. template <int PASS>
  72. void FFTRealPassInverse <PASS>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  73. {
  74. const long dist = 1L << (PASS - 1);
  75. const long c1_r = 0;
  76. const long c1_i = dist;
  77. const long c2_r = dist * 2;
  78. const long c2_i = dist * 3;
  79. const long cend = dist * 4;
  80. const long table_step = cos_len >> (PASS - 1);
  81. enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT };
  82. enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 };
  83. long coef_index = 0;
  84. do
  85. {
  86. const DataType * const sf = src_ptr + coef_index;
  87. DataType * const df = dest_ptr + coef_index;
  88. // Extreme coefficients are always real
  89. df [c1_r] = sf [c1_r] + sf [c2_r];
  90. df [c2_r] = sf [c1_r] - sf [c2_r];
  91. df [c1_i] = sf [c1_i] * 2;
  92. df [c2_i] = sf [c2_i] * 2;
  93. FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
  94. // Others are conjugate complex numbers
  95. for (long i = 1; i < dist; ++ i)
  96. {
  97. df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i];
  98. df [c1_i + i] = sf [c2_r + i] - sf [cend - i];
  99. DataType c;
  100. DataType s;
  101. FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
  102. osc_list [TRIGO_OSC],
  103. c,
  104. s,
  105. cos_ptr,
  106. i * table_step,
  107. (dist - i) * table_step
  108. );
  109. const DataType vr = sf [c1_r + i] - sf [c2_r - i];
  110. const DataType vi = sf [c2_r + i] + sf [cend - i];
  111. df [c2_r + i] = vr * c + vi * s;
  112. df [c2_i + i] = vi * c - vr * s;
  113. }
  114. coef_index += cend;
  115. }
  116. while (coef_index < len);
  117. }
  118. template <>
  119. inline void FFTRealPassInverse <2>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  120. {
  121. // Antepenultimate pass
  122. const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
  123. long coef_index = 0;
  124. do
  125. {
  126. dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
  127. dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
  128. dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2;
  129. dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2;
  130. dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3];
  131. dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7];
  132. const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3];
  133. const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7];
  134. dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2;
  135. dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2;
  136. coef_index += 8;
  137. }
  138. while (coef_index < len);
  139. }
  140. template <>
  141. inline void FFTRealPassInverse <1>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
  142. {
  143. // Penultimate and last pass at once
  144. const long qlen = len >> 2;
  145. long coef_index = 0;
  146. do
  147. {
  148. const long ri_0 = br_ptr [coef_index >> 2];
  149. const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2];
  150. const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2];
  151. const DataType b_1 = src_ptr [coef_index + 1] * 2;
  152. const DataType b_3 = src_ptr [coef_index + 3] * 2;
  153. dest_ptr [ri_0 ] = b_0 + b_1;
  154. dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1;
  155. dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3;
  156. dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3;
  157. coef_index += 4;
  158. }
  159. while (coef_index < len);
  160. }
  161. /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
  162. /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
  163. } // namespace ffft
  164. #endif // ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
  165. #undef ffft_FFTRealPassInverse_CURRENT_CODEHEADER
  166. /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/