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
| template<class G>
Vector<G> &complex_radix2_fft_in_place(Vector<G> &X)
{
typedef Vector<G> array_type;
typedef typename array_type::size_type size_type;
typedef PROMOTE2(complex<float>,typename array_type::value_type) value_type;
typedef PROMOTE2(complex<float>,typename value_type_traits<value_type>::value_type) complex_type;
typedef typename complex_type::value_type real_type;
size_type n = X.size();
assert(log2n(n)>=0);
bitreverse_order(X);
for (int dual=1; dual<n; dual*=2)
{
real_type theta = -PI/dual;
real_type s = sin(theta);
real_type s2 = real_type(2)*sqr(sin(real_type(0.5)*theta));
complex_type w(1,0);
value_type wd;
for (size_t a=0; a<dual; ++a, w+=s*complex_type(-w.imag(),w.real())-s2*w)
for (size_t b=0; b<n; b+=2*dual)
{
value_type &Ti=X[b+a];
value_type &Tj=X[b+a+dual];
wd=w*Tj;
Tj=Ti-wd;
Ti+=wd;
}
}
return X;
} |
Partager