ResidualVM logo ResidualVM website - Forums - Contact us BuildBot - Doxygen - Wiki curved edge

fft.cpp

Go to the documentation of this file.
00001 /* ScummVM - Graphic Adventure Engine
00002  *
00003  * ScummVM is the legal property of its developers, whose names
00004  * are too numerous to list here. Please refer to the COPYRIGHT
00005  * file distributed with this source distribution.
00006  *
00007  * This program is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU General Public License
00009  * as published by the Free Software Foundation; either version 2
00010  * of the License, or (at your option) any later version.
00011  *
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00020  *
00021  */
00022 
00023 // Based on eos' (I)FFT code which is in turn
00024 // Based upon the (I)FFT code in FFmpeg
00025 // Copyright (c) 2008 Loren Merritt
00026 // Copyright (c) 2002 Fabrice Bellard
00027 // Partly based on libdjbfft by D. J. Bernstein
00028 
00029 #include "common/cosinetables.h"
00030 #include "common/fft.h"
00031 #include "common/util.h"
00032 #include "common/textconsole.h"
00033 
00034 namespace Common {
00035 
00036 FFT::FFT(int bits, int inverse) : _bits(bits), _inverse(inverse) {
00037     assert((_bits >= 2) && (_bits <= 16));
00038 
00039     int n = 1 << bits;
00040     int nPoints;
00041 
00042     _tmpBuf = new Complex[n];
00043     _expTab = new Complex[n / 2];
00044     _revTab = new uint16[n];
00045 
00046     _splitRadix = 1;
00047 
00048     for (int i = 0; i < n; i++)
00049         _revTab[-splitRadixPermutation(i, n, _inverse) & (n - 1)] = i;
00050 
00051     for (int i = 0; i < ARRAYSIZE(_cosTables); i++) {
00052         if (i + 4 <= _bits) {
00053             nPoints = 1 << (i + 4);
00054             _cosTables[i] = new Common::CosineTable(nPoints);
00055         }
00056         else
00057             _cosTables[i] = nullptr;
00058     }
00059 }
00060 
00061 FFT::~FFT() {
00062     for (int i = 0; i < ARRAYSIZE(_cosTables); i++) {
00063         delete _cosTables[i];
00064     }
00065 
00066     delete[] _revTab;
00067     delete[] _expTab;
00068     delete[] _tmpBuf;
00069 }
00070 
00071 const uint16 *FFT::getRevTab() const {
00072     return _revTab;
00073 }
00074 
00075 void FFT::permute(Complex *z) {
00076     int np = 1 << _bits;
00077 
00078     if (_tmpBuf) {
00079         for (int j = 0; j < np; j++)
00080             _tmpBuf[_revTab[j]] = z[j];
00081 
00082         memcpy(z, _tmpBuf, np * sizeof(Complex));
00083 
00084         return;
00085     }
00086 
00087     // Reverse
00088     for (int j = 0; j < np; j++) {
00089         int k = _revTab[j];
00090 
00091         if (k < j)
00092             SWAP(z[k], z[j]);
00093     }
00094 }
00095 
00096 int FFT::splitRadixPermutation(int i, int n, int inverse) {
00097     if (n <= 2)
00098         return i & 1;
00099 
00100     int m = n >> 1;
00101 
00102     if (!(i & m))
00103         return splitRadixPermutation(i, m, inverse) * 2;
00104 
00105     m >>= 1;
00106 
00107     if (inverse == !(i & m))
00108         return splitRadixPermutation(i, m, inverse) * 4 + 1;
00109 
00110     return splitRadixPermutation(i, m, inverse) * 4 - 1;
00111 }
00112 
00113 #define sqrthalf (float)M_SQRT1_2
00114 
00115 #define BF(x, y, a, b) { \
00116     x = a - b; \
00117     y = a + b; \
00118 }
00119 
00120 #define BUTTERFLIES(a0, a1, a2, a3) { \
00121     BF(t3, t5, t5, t1); \
00122     BF(a2.re, a0.re, a0.re, t5); \
00123     BF(a3.im, a1.im, a1.im, t3); \
00124     BF(t4, t6, t2, t6); \
00125     BF(a3.re, a1.re, a1.re, t4); \
00126     BF(a2.im, a0.im, a0.im, t6); \
00127 }
00128 
00129 // force loading all the inputs before storing any.
00130 // this is slightly slower for small data, but avoids store->load aliasing
00131 // for addresses separated by large powers of 2.
00132 #define BUTTERFLIES_BIG(a0, a1, a2, a3) { \
00133     float r0 = a0.re, i0 = a0.im, r1 = a1.re, i1 = a1.im; \
00134     BF(t3, t5, t5, t1); \
00135     BF(a2.re, a0.re, r0, t5); \
00136     BF(a3.im, a1.im, i1, t3); \
00137     BF(t4, t6, t2, t6); \
00138     BF(a3.re, a1.re, r1, t4); \
00139     BF(a2.im, a0.im, i0, t6); \
00140 }
00141 
00142 #define TRANSFORM(a0, a1, a2, a3, wre, wim) { \
00143     t1 = a2.re * wre + a2.im * wim; \
00144     t2 = a2.im * wre - a2.re * wim; \
00145     t5 = a3.re * wre - a3.im * wim; \
00146     t6 = a3.im * wre + a3.re * wim; \
00147     BUTTERFLIES(a0, a1, a2, a3) \
00148 }
00149 
00150 #define TRANSFORM_ZERO(a0, a1, a2, a3) { \
00151     t1 = a2.re; \
00152     t2 = a2.im; \
00153     t5 = a3.re; \
00154     t6 = a3.im; \
00155     BUTTERFLIES(a0, a1, a2, a3) \
00156 }
00157 
00158 /* z[0...8n-1], w[1...2n-1] */
00159 #define PASS(name) \
00160 static void name(Complex *z, const float *wre, unsigned int n) { \
00161     float t1, t2, t3, t4, t5, t6; \
00162     int o1 = 2 * n; \
00163     int o2 = 4 * n; \
00164     int o3 = 6 * n; \
00165     const float *wim = wre + o1; \
00166     n--; \
00167     \
00168     TRANSFORM_ZERO(z[0], z[o1], z[o2], z[o3]); \
00169     TRANSFORM(z[1], z[o1 + 1], z[o2 + 1], z[o3 + 1], wre[1], wim[-1]); \
00170     do { \
00171         z += 2; \
00172         wre += 2; \
00173         wim -= 2; \
00174         TRANSFORM(z[0], z[o1], z[o2], z[o3], wre[0], wim[0]);\
00175         TRANSFORM(z[1], z[o1 + 1], z[o2 + 1], z[o3 + 1], wre[1], wim[-1]);\
00176     } while(--n);\
00177 }
00178 
00179 PASS(pass)
00180 #undef BUTTERFLIES
00181 #define BUTTERFLIES BUTTERFLIES_BIG
00182 PASS(pass_big)
00183 
00184 void FFT::fft4(Complex *z) {
00185     float t1, t2, t3, t4, t5, t6, t7, t8;
00186 
00187     BF(t3, t1, z[0].re, z[1].re);
00188     BF(t8, t6, z[3].re, z[2].re);
00189     BF(z[2].re, z[0].re, t1, t6);
00190     BF(t4, t2, z[0].im, z[1].im);
00191     BF(t7, t5, z[2].im, z[3].im);
00192     BF(z[3].im, z[1].im, t4, t8);
00193     BF(z[3].re, z[1].re, t3, t7);
00194     BF(z[2].im, z[0].im, t2, t5);
00195 }
00196 
00197 void FFT::fft8(Complex *z) {
00198     float t1, t2, t3, t4, t5, t6, t7, t8;
00199 
00200     fft4(z);
00201 
00202     BF(t1, z[5].re, z[4].re, -z[5].re);
00203     BF(t2, z[5].im, z[4].im, -z[5].im);
00204     BF(t3, z[7].re, z[6].re, -z[7].re);
00205     BF(t4, z[7].im, z[6].im, -z[7].im);
00206     BF(t8, t1, t3, t1);
00207     BF(t7, t2, t2, t4);
00208     BF(z[4].re, z[0].re, z[0].re, t1);
00209     BF(z[4].im, z[0].im, z[0].im, t2);
00210     BF(z[6].re, z[2].re, z[2].re, t7);
00211     BF(z[6].im, z[2].im, z[2].im, t8);
00212 
00213     TRANSFORM(z[1], z[3], z[5], z[7], sqrthalf, sqrthalf);
00214 }
00215 
00216 void FFT::fft16(Complex *z) {
00217     float t1, t2, t3, t4, t5, t6;
00218 
00219     fft8(z);
00220     fft4(z + 8);
00221     fft4(z + 12);
00222 
00223     assert(_cosTables[0]);
00224     const float * const cosTable = _cosTables[0]->getTable();
00225 
00226     TRANSFORM_ZERO(z[0], z[4], z[8], z[12]);
00227     TRANSFORM(z[2], z[6], z[10], z[14], sqrthalf, sqrthalf);
00228     TRANSFORM(z[1], z[5], z[9], z[13], cosTable[1],cosTable[3]);
00229     TRANSFORM(z[3], z[7], z[11], z[15], cosTable[3], cosTable[1]);
00230 }
00231 
00232 void FFT::fft(int n, int logn, Complex *z) {
00233     switch (logn) {
00234     case 2:
00235         fft4(z);
00236         break;
00237     case 3:
00238         fft8(z);
00239         break;
00240     case 4:
00241         fft16(z);
00242         break;
00243     default:
00244         fft((n / 2), logn - 1, z);
00245         fft((n / 4), logn - 2, z + (n / 4) * 2);
00246         fft((n / 4), logn - 2, z + (n / 4) * 3);
00247         assert(_cosTables[logn - 4]);
00248         if (n > 1024)
00249             pass_big(z, _cosTables[logn - 4]->getTable(), (n / 4) / 2);
00250         else
00251             pass(z, _cosTables[logn - 4]->getTable(), (n / 4) / 2);
00252     }
00253 }
00254 
00255 void FFT::calc(Complex *z) {
00256     fft(1 << _bits, _bits, z);
00257 }
00258 
00259 } // End of namespace Common


Generated on Sat Nov 30 2019 05:00:24 for ResidualVM by doxygen 1.7.1
curved edge   curved edge