1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
/* ScummVM - Graphic Adventure Engine
*
* ScummVM is the legal property of its developers, whose names
* are too numerous to list here. Please refer to the COPYRIGHT
* file distributed with this source distribution.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*/
// Based on eos' (I)RDFT code which is in turn
// Based upon the (I)RDFT code in FFmpeg
// Copyright (c) 2009 Alex Converse <alex dot converse at gmail dot com>
#include "math/rdft.h"
#include "math/fft.h"
#include "math/utils.h"
namespace Math {
RDFT::RDFT(int bits, TransformType trans) : _bits(bits), _sin(1 << bits), _cos(1 << bits), _fft(nullptr) {
assert((_bits >= 4) && (_bits <= 16));
_inverse = trans == IDFT_C2R || trans == DFT_C2R;
_signConvention = trans == IDFT_R2C || trans == DFT_C2R ? 1 : -1;
_fft = new FFT(bits - 1, trans == IDFT_C2R || trans == IDFT_R2C);
int n = 1 << bits;
_tSin = _sin.getTable() + (trans == DFT_R2C || trans == DFT_C2R) * (n >> 2);
_tCos = _cos.getTable();
}
RDFT::~RDFT() {
delete _fft;
}
void RDFT::calc(float *data) {
const int n = 1 << _bits;
const float k1 = 0.5f;
const float k2 = 0.5f - _inverse;
if (!_inverse) {
_fft->permute((Complex *)data);
_fft->calc ((Complex *)data);
}
Complex ev, od;
/* i=0 is a special case because of packing, the DC term is real, so we
are going to throw the N/2 term (also real) in with it. */
ev.re = data[0];
data[0] = ev.re + data[1];
data[1] = ev.re - data[1];
int i;
for (i = 1; i < (n >> 2); i++) {
int i1 = 2 * i;
int i2 = n - i1;
/* Separate even and odd FFTs */
ev.re = k1 * (data[i1 ] + data[i2 ]);
od.im = -k2 * (data[i1 ] - data[i2 ]);
ev.im = k1 * (data[i1 + 1] - data[i2 + 1]);
od.re = k2 * (data[i1 + 1] + data[i2 + 1]);
/* Apply twiddle factors to the odd FFT and add to the even FFT */
data[i1 ] = ev.re + od.re * _tCos[i] - od.im * _tSin[i];
data[i1 + 1] = ev.im + od.im * _tCos[i] + od.re * _tSin[i];
data[i2 ] = ev.re - od.re * _tCos[i] + od.im * _tSin[i];
data[i2 + 1] = -ev.im + od.im * _tCos[i] + od.re * _tSin[i];
}
data[2 * i + 1] = _signConvention * data[2 * i + 1];
if (_inverse) {
data[0] *= k1;
data[1] *= k1;
_fft->permute((Complex *)data);
_fft->calc ((Complex *)data);
}
}
} // End of namespace Math
|
