| 127 | */ |
| 128 | |
| 129 | void FFT(size_t NumSamples, |
| 130 | bool InverseTransform, |
| 131 | const float *RealIn, const float *ImagIn, |
| 132 | float *RealOut, float *ImagOut) |
| 133 | { |
| 134 | double angle_numerator = 2.0 * M_PI; |
| 135 | double tr, ti; /* temp real, temp imaginary */ |
| 136 | |
| 137 | if (!IsPowerOfTwo(NumSamples)) { |
| 138 | wxFprintf(stderr, "%ld is not a power of two\n", NumSamples); |
| 139 | exit(1); |
| 140 | } |
| 141 | |
| 142 | if (!gFFTBitTable) |
| 143 | InitFFT(); |
| 144 | |
| 145 | if (!InverseTransform) |
| 146 | angle_numerator = -angle_numerator; |
| 147 | |
| 148 | /* Number of bits needed to store indices */ |
| 149 | auto NumBits = NumberOfBitsNeeded(NumSamples); |
| 150 | |
| 151 | /* |
| 152 | ** Do simultaneous data copy and bit-reversal ordering into outputs... |
| 153 | */ |
| 154 | |
| 155 | for (size_t i = 0; i < NumSamples; i++) { |
| 156 | auto j = FastReverseBits(i, NumBits); |
| 157 | RealOut[j] = RealIn[i]; |
| 158 | ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i]; |
| 159 | } |
| 160 | |
| 161 | /* |
| 162 | ** Do the FFT itself... |
| 163 | */ |
| 164 | |
| 165 | size_t BlockEnd = 1; |
| 166 | for (size_t BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) { |
| 167 | |
| 168 | double delta_angle = angle_numerator / (double) BlockSize; |
| 169 | |
| 170 | double sm2 = sin(-2 * delta_angle); |
| 171 | double sm1 = sin(-delta_angle); |
| 172 | double cm2 = cos(-2 * delta_angle); |
| 173 | double cm1 = cos(-delta_angle); |
| 174 | double w = 2 * cm1; |
| 175 | double ar0, ar1, ar2, ai0, ai1, ai2; |
| 176 | |
| 177 | for (size_t i = 0; i < NumSamples; i += BlockSize) { |
| 178 | ar2 = cm2; |
| 179 | ar1 = cm1; |
| 180 | |
| 181 | ai2 = sm2; |
| 182 | ai1 = sm1; |
| 183 | |
| 184 | for (size_t j = i, n = 0; n < BlockEnd; j++, n++) { |
| 185 | ar0 = w * ar1 - ar2; |
| 186 | ar2 = ar1; |