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

mdct.cpp

Go to the documentation of this file.
00001 /* ResidualVM - A 3D game interpreter
00002  *
00003  * ResidualVM is the legal property of its developers, whose names
00004  * are too numerous to list here. Please refer to the AUTHORS
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 xoreos' (I)RDFT code which is in turn
00024 // Based upon the (I)MDCT code in FFmpeg
00025 // Copyright (c) 2002 Fabrice Bellard
00026 
00031 #include "common/math.h"
00032 #include "common/util.h"
00033 #include "common/fft.h"
00034 #include "common/mdct.h"
00035 
00036 namespace Common {
00037 
00038 MDCT::MDCT(int bits, bool inverse, double scale) : _bits(bits), _fft(0) {
00039     _size = 1 << bits;
00040 
00041     _fft = new FFT(_bits - 2, inverse);
00042 
00043     const int size2 = _size >> 1;
00044     const int size4 = _size >> 2;
00045 
00046     _tCos = new float[size2];
00047     _tSin = _tCos + size4;
00048 
00049     const double theta = 1.0 / 8.0 + (scale < 0 ? size4 : 0);
00050 
00051     scale = sqrt(ABS(scale));
00052     for (int i = 0; i < size4; i++) {
00053         const double alpha = 2 * M_PI * (i + theta) / _size;
00054 
00055         _tCos[i] = -cos(alpha) * scale;
00056         _tSin[i] = -sin(alpha) * scale;
00057     }
00058 }
00059 
00060 MDCT::~MDCT() {
00061     delete[] _tCos;
00062     delete _fft;
00063 }
00064 
00065 #define CMUL(dre, dim, are, aim, bre, bim) do { \
00066         (dre) = (are) * (bre) - (aim) * (bim);  \
00067         (dim) = (are) * (bim) + (aim) * (bre);  \
00068     } while (0)
00069 
00070 void MDCT::calcMDCT(float *output, const float *input) {
00071     Complex *x = (Complex *) output;
00072 
00073     const int size2 = _size >> 1;
00074     const int size4 = _size >> 2;
00075     const int size8 = _size >> 3;
00076     const int size3 = _size * 3;
00077 
00078     const uint16 *revTab = _fft->getRevTab();
00079 
00080     // Pre rotation
00081     for (int i = 0; i < size8; i++) {
00082         float re = -input[2 * i + size3] - input[size3 - 1 - 2 * i];
00083         float im = -input[2 * i + size4] + input[size4 - 1 - 2 * i];
00084         int    j = revTab[i];
00085 
00086         CMUL(x[j].re, x[j].im, re, im, -_tCos[i], _tSin[i]);
00087 
00088         re =  input[2 * i        ] - input[size2 - 1 - 2 * i];
00089         im = -input[2 * i + size2] - input[_size - 1 - 2 * i];
00090          j = revTab[size8 + i];
00091 
00092         CMUL(x[j].re, x[j].im, re, im, -_tCos[size8 + i], _tSin[size8 + i]);
00093     }
00094 
00095     _fft->calc(x);
00096 
00097     // Post rotation
00098     for (int i = 0; i < size8; i++) {
00099         float r0, i0, r1, i1;
00100 
00101         CMUL(i1, r0, x[size8-i-1].re, x[size8-i-1].im, -_tSin[size8-i-1], -_tCos[size8-i-1]);
00102         CMUL(i0, r1, x[size8+i  ].re, x[size8+i  ].im, -_tSin[size8+i  ], -_tCos[size8+i  ]);
00103 
00104         x[size8 - i - 1].re = r0;
00105         x[size8 - i - 1].im = i0;
00106         x[size8 + i    ].re = r1;
00107         x[size8 + i    ].im = i1;
00108     }
00109 }
00110 
00111 void MDCT::calcIMDCT(float *output, const float *input) {
00112     const int size2 = _size >> 1;
00113     const int size4 = _size >> 2;
00114 
00115     calcHalfIMDCT(output + size4, input);
00116 
00117     for (int k = 0; k < size4; k++) {
00118         output[        k    ] = -output[size2 - k - 1];
00119         output[_size - k - 1] =  output[size2 + k    ];
00120     }
00121 }
00122 
00123 void MDCT::calcHalfIMDCT(float *output, const float *input) {
00124     Complex *z = (Complex *) output;
00125 
00126     const int size2 = _size >> 1;
00127     const int size4 = _size >> 2;
00128     const int size8 = _size >> 3;
00129 
00130     const uint16 *revTab = _fft->getRevTab();
00131 
00132     // Pre rotation
00133     const float *in1 = input;
00134     const float *in2 = input + size2 - 1;
00135     for (int k = 0; k < size4; k++) {
00136         const int j = revTab[k];
00137 
00138         CMUL(z[j].re, z[j].im, *in2, *in1, _tCos[k], _tSin[k]);
00139 
00140         in1 += 2;
00141         in2 -= 2;
00142     }
00143 
00144     _fft->calc(z);
00145 
00146     // Post rotation + reordering
00147     for (int k = 0; k < size8; k++) {
00148         float r0, i0, r1, i1;
00149 
00150         CMUL(r0, i1, z[size8-k-1].im, z[size8-k-1].re, _tSin[size8-k-1], _tCos[size8-k-1]);
00151         CMUL(r1, i0, z[size8+k  ].im, z[size8+k  ].re, _tSin[size8+k  ], _tCos[size8+k  ]);
00152 
00153         z[size8 - k - 1].re = r0;
00154         z[size8 - k - 1].im = i0;
00155         z[size8 + k    ].re = r1;
00156         z[size8 + k    ].im = i1;
00157     }
00158 }
00159 
00160 } // End of namespace Common


Generated on Sat Feb 16 2019 05:00:57 for ResidualVM by doxygen 1.7.1
curved edge   curved edge