1 /*
2 * TinyFFT.cpp
3 * -----------
4 * Purpose: A simple FFT implementation for power-of-two FFTs
5 * Notes : This is a C++ adaption of Ryuhei Mori's BSD 2-clause licensed TinyFFT
6 * available from https://github.com/ryuhei-mori/tinyfft
7 * Authors: Ryuhei Mori
8 * OpenMPT Devs
9 * The OpenMPT source code is released under the BSD license. Read LICENSE for more details.
10 */
11
12
13 #include "stdafx.h"
14 #include "TinyFFT.h"
15
16 OPENMPT_NAMESPACE_BEGIN
17
GenerateTwiddleFactors(uint32 i,uint32 b,std::complex<double> z)18 void TinyFFT::GenerateTwiddleFactors(uint32 i, uint32 b, std::complex<double> z)
19 {
20 if(b == 0)
21 w[i] = z;
22 else
23 {
24 GenerateTwiddleFactors(i, b >> 1, z);
25 GenerateTwiddleFactors(i | b, b >> 1, z * w[b]);
26 }
27 }
28
29
TinyFFT(const uint32 fftSize)30 TinyFFT::TinyFFT(const uint32 fftSize)
31 : w(std::size_t(1) << (fftSize - 1))
32 , k(fftSize)
33 {
34 const uint32 m = 1 << k;
35 constexpr double PI2_ = 6.28318530717958647692;
36 const double arg = -PI2_ / m;
37 for(uint32 i = 1, j = m / 4; j; i <<= 1, j >>= 1)
38 {
39 w[i] = std::exp(I * (arg * j));
40 }
41 GenerateTwiddleFactors(0, m / 4, 1);
42 }
43
44
Size() const45 uint32 TinyFFT::Size() const noexcept
46 {
47 return 1 << k;
48 }
49
50
51 // Computes in-place FFT of size 2^k of A, result is in bit-reversed order.
FFT(std::vector<std::complex<double>> & A) const52 void TinyFFT::FFT(std::vector<std::complex<double>> &A) const
53 {
54 MPT_ASSERT(A.size() == (std::size_t(1) << k));
55 const uint32 m = 1 << k;
56 uint32 u = 1;
57 uint32 v = m / 4;
58 if(k & 1)
59 {
60 for(uint32 j = 0; j < m / 2; j++)
61 {
62 auto Ajv = A[j + (m / 2)];
63 A[j + (m / 2)] = A[j] - Ajv;
64 A[j] += Ajv;
65 }
66 u <<= 1;
67 v >>= 1;
68 }
69 for(uint32 i = k & ~1; i > 0; i -= 2)
70 {
71 for(uint32 jh = 0; jh < u; jh++)
72 {
73 auto wj = w[jh << 1];
74 auto wj2 = w[jh];
75 auto wj3 = wj2 * wj;
76 for(uint32 j = jh << i, je = j + v; j < je; j++)
77 {
78 auto tmp0 = A[j];
79 auto tmp1 = wj * A[j + v];
80 auto tmp2 = wj2 * A[j + 2 * v];
81 auto tmp3 = wj3 * A[j + 3 * v];
82
83 auto ttmp0 = tmp0 + tmp2;
84 auto ttmp2 = tmp0 - tmp2;
85 auto ttmp1 = tmp1 + tmp3;
86 auto ttmp3 = -I * (tmp1 - tmp3);
87
88 A[j] = ttmp0 + ttmp1;
89 A[j + v] = ttmp0 - ttmp1;
90 A[j + 2 * v] = ttmp2 + ttmp3;
91 A[j + 3 * v] = ttmp2 - ttmp3;
92 }
93 }
94 u <<= 2;
95 v >>= 2;
96 }
97 }
98
99
100 // Computes in-place IFFT of size 2^k of A, input is expected to be in bit-reversed order.
IFFT(std::vector<std::complex<double>> & A) const101 void TinyFFT::IFFT(std::vector<std::complex<double>> &A) const
102 {
103 MPT_ASSERT(A.size() == (std::size_t(1) << k));
104 const uint32 m = 1 << k;
105 uint32 u = m / 4;
106 uint32 v = 1;
107 for(uint32 i = 2; i <= k; i += 2)
108 {
109 for(uint32 jh = 0; jh < u; jh++)
110 {
111 auto wj = std::conj(w[jh << 1]);
112 auto wj2 = std::conj(w[jh]);
113 auto wj3 = wj2 * wj;
114 for(uint32 j = jh << i, je = j + v; j < je; j++)
115 {
116 auto tmp0 = A[j];
117 auto tmp1 = A[j + v];
118 auto tmp2 = A[j + 2 * v];
119 auto tmp3 = A[j + 3 * v];
120
121 auto ttmp0 = tmp0 + tmp1;
122 auto ttmp1 = tmp0 - tmp1;
123 auto ttmp2 = tmp2 + tmp3;
124 auto ttmp3 = I * (tmp2 - tmp3);
125
126 A[j] = ttmp0 + ttmp2;
127 A[j + v] = wj * (ttmp1 + ttmp3);
128 A[j + 2 * v] = wj2 * (ttmp0 - ttmp2);
129 A[j + 3 * v] = wj3 * (ttmp1 - ttmp3);
130 }
131 }
132 u >>= 2;
133 v <<= 2;
134 }
135 if(k & 1)
136 {
137 for(uint32 j = 0; j < m / 2; j++)
138 {
139 auto Ajv = A[j + (m / 2)];
140 A[j + (m / 2)] = A[j] - Ajv;
141 A[j] += Ajv;
142 }
143 }
144 }
145
146
Normalize(std::vector<std::complex<double>> & data)147 void TinyFFT::Normalize(std::vector<std::complex<double>> &data)
148 {
149 const double s = static_cast<double>(data.size());
150 for(auto &v : data)
151 v /= s;
152 }
153
154 OPENMPT_NAMESPACE_END
155