diff --git a/fft.hh b/fft.hh index 4a2cf01..65db4c3 100644 --- a/fft.hh +++ b/fft.hh @@ -1125,5 +1125,40 @@ public: } }; +template +class RealToHalfComplexTransform +{ + static_assert(BINS%2==0, "BINS must be even"); + static const int N = BINS / 2; + TYPE factors[N]; + TYPE A[N], B[N]; +public: + typedef typename TYPE::value_type value_type; + RealToHalfComplexTransform() + { + for (int n = 0; n < N; ++n) + factors[n] = TYPE(UnitCircle::cos(n, N), -UnitCircle::sin(n, N)); + for (int n = 0; n < N; ++n) { + TYPE sincos( + UnitCircle::sin(n, BINS), + UnitCircle::cos(n, BINS) + ); + A[n] = value_type(0.5) * (TYPE(1) - sincos); + B[n] = value_type(0.5) * (TYPE(1) + sincos); + } + } + inline void operator ()(TYPE *out, const value_type *in) + { + FFT::Dit::dit(out, reinterpret_cast(in), factors); + out[N] = value_type(0.5) * (out[0].real() - out[0].imag()); + out[0] = value_type(0.5) * (out[0].real() + out[0].imag()); + for (int i = 1; i <= N/2; ++i) { + TYPE tmp = out[i]*A[i] + conj(out[N-i])*B[i]; + out[N-i] = out[N-i]*A[N-i] + conj(out[i])*B[N-i]; + out[i] = tmp; + } + } +}; + }