1 #ifndef INCLUDE_KISSFFT_H 2 #define INCLUDE_KISSFFT_H 40 template<
typename T_scalar,
typename T_complex>
46 T_scalar phinc = (inverse ? 2 : -2) * acos((T_scalar)-1) / nfft;
47 for(
int i = 0;
i < nfft; ++
i)
48 dst[
i] =
exp(std::complex<T_scalar>(0,
i * phinc));
51 void prepare(std::vector<std::complex<T_scalar> >& dst,
int nfft,
bool inverse, std::vector<int>& stageRadix, std::vector<int>& stageRemainder)
78 stageRadix.push_back(p);
79 stageRemainder.push_back(n);
92 template<
typename T_Scalar,
typename T_Complex,
typename T_traits = kissfft_utils::traits<T_Scalar, T_Complex> >
103 kissfft(
int nfft,
bool inverse,
const traits_type & traits = traits_type()) :
104 _nfft(nfft), _inverse(inverse), _traits(traits)
106 _traits.prepare(
_twiddles, _nfft, _inverse, _stageRadix, _stageRemainder);
109 void configure(
int nfft,
bool inverse,
const traits_type & traits = traits_type())
113 _stageRemainder.clear();
118 _traits.prepare(
_twiddles, _nfft, _inverse, _stageRadix, _stageRemainder);
123 kf_work(0, dst, src, 1, 1);
127 void kf_work(
int stage, cpx_type* Fout,
const cpx_type* f,
size_t fstride,
size_t in_stride)
129 int p = _stageRadix[stage];
130 int m = _stageRemainder[stage];
131 cpx_type * Fout_beg = Fout;
132 cpx_type * Fout_end = Fout + p * m;
137 f += fstride * in_stride;
138 }
while(++Fout != Fout_end);
145 kf_work(stage + 1, Fout, f, fstride * p, in_stride);
146 f += fstride * in_stride;
147 }
while((Fout += m) != Fout_end);
155 kf_bfly2(Fout, fstride, m);
158 kf_bfly3(Fout, fstride, m);
161 kf_bfly4(Fout, fstride, m);
164 kf_bfly5(Fout, fstride, m);
167 kf_bfly_generic(Fout, fstride, m, p);
173 void C_ADD(cpx_type& c,
const cpx_type& a,
const cpx_type& b)
177 void C_MUL(cpx_type& c,
const cpx_type& a,
const cpx_type& b)
180 c =
cpx_type(a.real() * b.real() - a.imag() * b.imag(), a.real() * b.imag() + a.imag() * b.real());
182 void C_SUB(cpx_type& c,
const cpx_type& a,
const cpx_type& b)
193 scalar_type
S_MUL(
const scalar_type& a,
const scalar_type& b)
206 void kf_bfly2(cpx_type* Fout,
const size_t fstride,
int m)
208 for(
int k = 0; k < m; ++k) {
211 C_MUL(t, Fout[m + k], _traits.twiddle(k * fstride));
212 Fout[m + k] = Fout[k] - t;
217 void kf_bfly4(cpx_type* Fout,
const size_t fstride,
const size_t m)
220 int negative_if_inverse = _inverse * -2 + 1;
221 for(
size_t k = 0; k < m; ++k) {
223 C_MUL(scratch[0], Fout[k + m], _traits.twiddle(k * fstride));
224 C_MUL(scratch[1], Fout[k + 2 * m], _traits.twiddle(k * fstride * 2));
225 C_MUL(scratch[2], Fout[k + 3 * m], _traits.twiddle(k * fstride * 3));
226 scratch[5] = Fout[k] - scratch[1];
228 Fout[k] += scratch[1];
229 scratch[3] = scratch[0] + scratch[2];
230 scratch[4] = scratch[0] - scratch[2];
231 scratch[4] =
cpx_type(scratch[4].imag() * negative_if_inverse, -scratch[4].real() * negative_if_inverse);
233 Fout[k + 2 * m] = Fout[k] - scratch[3];
234 Fout[k] += scratch[3];
235 Fout[k + m] = scratch[5] + scratch[4];
236 Fout[k + 3 * m] = scratch[5] - scratch[4];
240 void kf_bfly3(cpx_type* Fout,
const size_t fstride,
const size_t m)
243 const size_t m2 = 2 * m;
253 C_FIXDIV(Fout[m], 3);
254 C_FIXDIV(Fout[m2], 3);
256 C_MUL(scratch[1], Fout[m], *tw1);
257 C_MUL(scratch[2], Fout[m2], *tw2);
259 C_ADD(scratch[3], scratch[1], scratch[2]);
260 C_SUB(scratch[0], scratch[1], scratch[2]);
264 Fout[m] =
cpx_type(Fout->real() - HALF_OF(scratch[3].real()), Fout->imag() - HALF_OF(scratch[3].imag()));
266 C_MULBYSCALAR(scratch[0], epi3.imag());
268 C_ADDTO(*Fout, scratch[3]);
270 Fout[m2] =
cpx_type(Fout[m].real() + scratch[0].imag(), Fout[m].imag() - scratch[0].real());
272 C_ADDTO(Fout[m],
cpx_type(-scratch[0].imag(), scratch[0].real()));
277 void kf_bfly5(cpx_type* Fout,
const size_t fstride,
const size_t m)
285 cpx_type scratch[13];
289 ya = twiddles[fstride * m];
290 yb = twiddles[fstride * 2 * m];
294 Fout2 = Fout0 + 2 * m;
295 Fout3 = Fout0 + 3 * m;
296 Fout4 = Fout0 + 4 * m;
299 for(u = 0; u < m; ++u) {
307 C_MUL(scratch[1], *Fout1, tw[u * fstride]);
308 C_MUL(scratch[2], *Fout2, tw[2 * u * fstride]);
309 C_MUL(scratch[3], *Fout3, tw[3 * u * fstride]);
310 C_MUL(scratch[4], *Fout4, tw[4 * u * fstride]);
312 C_ADD(scratch[7], scratch[1], scratch[4]);
313 C_SUB(scratch[10], scratch[1], scratch[4]);
314 C_ADD(scratch[8], scratch[2], scratch[3]);
315 C_SUB(scratch[9], scratch[2], scratch[3]);
317 C_ADDTO(*Fout0, scratch[7]);
318 C_ADDTO(*Fout0, scratch[8]);
320 scratch[5] = scratch[0] +
cpx_type(S_MUL(scratch[7].real(), ya.real()) + S_MUL(scratch[8].real(), yb.real()), S_MUL(scratch[7].imag(), ya.real())
321 + S_MUL(scratch[8].imag(), yb.real()));
323 scratch[6] =
cpx_type(S_MUL(scratch[10].imag(), ya.imag()) + S_MUL(scratch[9].imag(), yb.imag()), -S_MUL(scratch[10].real(), ya.imag()) - S_MUL(
324 scratch[9].real(), yb.imag()));
326 C_SUB(*Fout1, scratch[5], scratch[6]);
327 C_ADD(*Fout4, scratch[5], scratch[6]);
329 scratch[11] = scratch[0] +
cpx_type(S_MUL(scratch[7].real(), yb.real()) + S_MUL(scratch[8].real(), ya.real()), S_MUL(scratch[7].imag(), yb.real())
330 + S_MUL(scratch[8].imag(), ya.real()));
332 scratch[12] =
cpx_type(-S_MUL(scratch[10].imag(), yb.imag()) + S_MUL(scratch[9].imag(), ya.imag()), S_MUL(scratch[10].real(), yb.imag()) - S_MUL(
333 scratch[9].real(), ya.imag()));
335 C_ADD(*Fout2, scratch[11], scratch[12]);
336 C_SUB(*Fout3, scratch[11], scratch[12]);
356 cpx_type* scratchbuf =
new cpx_type[p];
358 for(u = 0; u < m; ++u) {
360 for(q1 = 0; q1 < p; ++q1) {
361 scratchbuf[q1] = Fout[k];
362 C_FIXDIV(scratchbuf[q1], p);
367 for(q1 = 0; q1 < p; ++q1) {
369 Fout[k] = scratchbuf[0];
370 for(q = 1; q < p; ++q) {
371 twidx += fstride * k;
374 C_MUL(t, scratchbuf[q], twiddles[twidx]);
void C_SUB(cpx_type &c, const cpx_type &a, const cpx_type &b)
traits_type::cpx_type cpx_type
void fill_twiddles(std::complex< T_scalar > *dst, int nfft, bool inverse)
void kf_bfly3(cpx_type *Fout, const size_t fstride, const size_t m)
scalar_type S_MUL(const scalar_type &a, const scalar_type &b)
scalar_type HALF_OF(const scalar_type &a)
std::vector< int > _stageRadix
void kf_bfly_generic(cpx_type *Fout, const size_t fstride, int m, int p)
Fixed< IntType, IntBits > exp(Fixed< IntType, IntBits > const &x)
kissfft(int nfft, bool inverse, const traits_type &traits=traits_type())
std::vector< cpx_type > _twiddles
void C_ADD(cpx_type &c, const cpx_type &a, const cpx_type &b)
void prepare(std::vector< std::complex< T_scalar > > &dst, int nfft, bool inverse, std::vector< int > &stageRadix, std::vector< int > &stageRemainder)
void C_MUL(cpx_type &c, const cpx_type &a, const cpx_type &b)
void configure(int nfft, bool inverse, const traits_type &traits=traits_type())
traits_type::scalar_type scalar_type
std::vector< int > _stageRemainder
void C_FIXDIV(cpx_type &, int)
void kf_bfly5(cpx_type *Fout, const size_t fstride, const size_t m)
const cpx_type twiddle(int i)
void kf_bfly4(cpx_type *Fout, const size_t fstride, const size_t m)
void transform(const cpx_type *src, cpx_type *dst)
void C_MULBYSCALAR(cpx_type &c, const scalar_type &a)
std::vector< cpx_type > _twiddles
void kf_bfly2(cpx_type *Fout, const size_t fstride, int m)
void kf_work(int stage, cpx_type *Fout, const cpx_type *f, size_t fstride, size_t in_stride)
void C_ADDTO(cpx_type &c, const cpx_type &a)