| 1 | /***************************************************************************** |
| 2 | |
| 3 | FFTRealPassInverse.hpp |
| 4 | Copyright (c) 2005 Laurent de Soras |
| 5 | |
| 6 | --- Legal stuff --- |
| 7 | |
| 8 | This library is free software; you can redistribute it and/or |
| 9 | modify it under the terms of the GNU Lesser General Public |
| 10 | License as published by the Free Software Foundation; either |
| 11 | version 2.1 of the License, or (at your option) any later version. |
| 12 | |
| 13 | This library is distributed in the hope that it will be useful, |
| 14 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 16 | Lesser General Public License for more details. |
| 17 | |
| 18 | You should have received a copy of the GNU Lesser General Public |
| 19 | License along with this library; if not, write to the Free Software |
| 20 | Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
| 21 | |
| 22 | *Tab=3***********************************************************************/ |
| 23 | |
| 24 | |
| 25 | |
| 26 | #if defined (FFTRealPassInverse_CURRENT_CODEHEADER) |
| 27 | #error Recursive inclusion of FFTRealPassInverse code header. |
| 28 | #endif |
| 29 | #define |
| 30 | |
| 31 | #if ! defined (FFTRealPassInverse_CODEHEADER_INCLUDED) |
| 32 | #define |
| 33 | |
| 34 | |
| 35 | |
| 36 | /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
| 37 | |
| 38 | #include "FFTRealUseTrigo.h" |
| 39 | |
| 40 | |
| 41 | |
| 42 | /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
| 43 | |
| 44 | |
| 45 | |
| 46 | template <int PASS> |
| 47 | 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 []) |
| 48 | { |
| 49 | process_internal ( |
| 50 | len, |
| 51 | dest_ptr, |
| 52 | src_ptr: f_ptr, |
| 53 | cos_ptr, |
| 54 | cos_len, |
| 55 | br_ptr, |
| 56 | osc_list |
| 57 | ); |
| 58 | FFTRealPassInverse <PASS - 1>::process_rec ( |
| 59 | len, |
| 60 | src_ptr, |
| 61 | dest_ptr, |
| 62 | cos_ptr, |
| 63 | cos_len, |
| 64 | br_ptr, |
| 65 | osc_list |
| 66 | ); |
| 67 | } |
| 68 | |
| 69 | |
| 70 | |
| 71 | template <int PASS> |
| 72 | 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 []) |
| 73 | { |
| 74 | process_internal ( |
| 75 | len, |
| 76 | dest_ptr, |
| 77 | src_ptr, |
| 78 | cos_ptr, |
| 79 | cos_len, |
| 80 | br_ptr, |
| 81 | osc_list |
| 82 | ); |
| 83 | FFTRealPassInverse <PASS - 1>::process_rec ( |
| 84 | len, |
| 85 | src_ptr, |
| 86 | dest_ptr, |
| 87 | cos_ptr, |
| 88 | cos_len, |
| 89 | br_ptr, |
| 90 | osc_list |
| 91 | ); |
| 92 | } |
| 93 | |
| 94 | template <> |
| 95 | 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 []) |
| 96 | { |
| 97 | // Stops recursion |
| 98 | } |
| 99 | |
| 100 | |
| 101 | |
| 102 | template <int PASS> |
| 103 | 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 []) |
| 104 | { |
| 105 | const long dist = 1L << (PASS - 1); |
| 106 | const long c1_r = 0; |
| 107 | const long c1_i = dist; |
| 108 | const long c2_r = dist * 2; |
| 109 | const long c2_i = dist * 3; |
| 110 | const long cend = dist * 4; |
| 111 | const long table_step = cos_len >> (PASS - 1); |
| 112 | |
| 113 | enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT }; |
| 114 | enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 }; |
| 115 | |
| 116 | long coef_index = 0; |
| 117 | do |
| 118 | { |
| 119 | const DataType * const sf = src_ptr + coef_index; |
| 120 | DataType * const df = dest_ptr + coef_index; |
| 121 | |
| 122 | // Extreme coefficients are always real |
| 123 | df [c1_r] = sf [c1_r] + sf [c2_r]; |
| 124 | df [c2_r] = sf [c1_r] - sf [c2_r]; |
| 125 | df [c1_i] = sf [c1_i] * 2; |
| 126 | df [c2_i] = sf [c2_i] * 2; |
| 127 | |
| 128 | FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]); |
| 129 | |
| 130 | // Others are conjugate complex numbers |
| 131 | for (long i = 1; i < dist; ++ i) |
| 132 | { |
| 133 | df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i]; |
| 134 | df [c1_i + i] = sf [c2_r + i] - sf [cend - i]; |
| 135 | |
| 136 | DataType c; |
| 137 | DataType s; |
| 138 | FFTRealUseTrigo <TRIGO_DIRECT>::iterate ( |
| 139 | osc_list [TRIGO_OSC], |
| 140 | c, |
| 141 | s, |
| 142 | cos_ptr, |
| 143 | i * table_step, |
| 144 | (dist - i) * table_step |
| 145 | ); |
| 146 | |
| 147 | const DataType vr = sf [c1_r + i] - sf [c2_r - i]; |
| 148 | const DataType vi = sf [c2_r + i] + sf [cend - i]; |
| 149 | |
| 150 | df [c2_r + i] = vr * c + vi * s; |
| 151 | df [c2_i + i] = vi * c - vr * s; |
| 152 | } |
| 153 | |
| 154 | coef_index += cend; |
| 155 | } |
| 156 | while (coef_index < len); |
| 157 | } |
| 158 | |
| 159 | template <> |
| 160 | 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 []) |
| 161 | { |
| 162 | // Antepenultimate pass |
| 163 | const DataType sqrt2_2 = DataType (SQRT2 * 0.5); |
| 164 | |
| 165 | long coef_index = 0; |
| 166 | do |
| 167 | { |
| 168 | dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4]; |
| 169 | dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4]; |
| 170 | dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2; |
| 171 | dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2; |
| 172 | |
| 173 | dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3]; |
| 174 | dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7]; |
| 175 | |
| 176 | const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3]; |
| 177 | const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7]; |
| 178 | |
| 179 | dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2; |
| 180 | dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2; |
| 181 | |
| 182 | coef_index += 8; |
| 183 | } |
| 184 | while (coef_index < len); |
| 185 | } |
| 186 | |
| 187 | template <> |
| 188 | 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 []) |
| 189 | { |
| 190 | // Penultimate and last pass at once |
| 191 | const long qlen = len >> 2; |
| 192 | |
| 193 | long coef_index = 0; |
| 194 | do |
| 195 | { |
| 196 | const long ri_0 = br_ptr [coef_index >> 2]; |
| 197 | |
| 198 | const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2]; |
| 199 | const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2]; |
| 200 | const DataType b_1 = src_ptr [coef_index + 1] * 2; |
| 201 | const DataType b_3 = src_ptr [coef_index + 3] * 2; |
| 202 | |
| 203 | dest_ptr [ri_0 ] = b_0 + b_1; |
| 204 | dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1; |
| 205 | dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3; |
| 206 | dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3; |
| 207 | |
| 208 | coef_index += 4; |
| 209 | } |
| 210 | while (coef_index < len); |
| 211 | } |
| 212 | |
| 213 | |
| 214 | |
| 215 | /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
| 216 | |
| 217 | |
| 218 | |
| 219 | /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
| 220 | |
| 221 | |
| 222 | |
| 223 | #endif // FFTRealPassInverse_CODEHEADER_INCLUDED |
| 224 | |
| 225 | #undef FFTRealPassInverse_CURRENT_CODEHEADER |
| 226 | |
| 227 | |
| 228 | |
| 229 | /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
| 230 | |