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 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
|
template<class G>
void real_radix2_fft_in_place(Vector<G> &X)
{
typedef Vector<G> array_type;
typedef typename array_type::size_type size_type;
typedef PROMOTE2(float,typename array_type::value_type) value_type;
typedef PROMOTE2(complex<float>,typename array_type::value_type) complex_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 (size_t p_1=1, p=2, q=n/2; p_1<n; p_1=p, p*=2, q/=2)
{
complex_type w(1,0);
real_type theta = -2*PI/p;
real_type s = sin (theta);
real_type s2 = 2*sqr(sin(theta/2));
{
value_type t0;
for (size_type b=0; b<q; b++)
{
typename array_type::reference r0 = X[b*p];
typename array_type::reference r1 = X[b*p+p_1];
t0 = r0+r1;
r1 = r0-r1;
r0 = t0;
}
}
complex_value_type z0, z1, z2;
for (size_t a=1; a<p_1/2; ++a)
{
w+=s*complex_type(-w.imag(),w.real())-s2*w;
for (size_t b=0; b<q; ++b)
{
typename array_type::reference r0 = X[b*p+a];
typename array_type::reference r1 = X[b*p+p_1-a];
typename array_type::reference r2 = X[b*p+p_1+a];
typename array_type::reference r3 = X[b*p+p-a];
real(z0)=r2; imag(z0)=r3;
real(z1)=r0; imag(z1)=r1;
z0 *= w;
z2 = z1+z0;
z1 -= z0;
r0 = real(z2);
r3 = imag(z2);
r1 = real(z1);
r2 = -imag(z1);
}
}
if (p_1>1) for (size_t b=0; b<q; b++) X[b*p+p-p_1/2]*=double(-1);
}
}
template<class G>
void real_radix2_ifft_in_place(Vector<G> &X)
{
typedef Vector<G> array_type;
typedef typename array_type::size_type size_type;
typedef PROMOTE2(float,typename array_type::value_type) value_type;
typedef PROMOTE2(complex<float>,typename array_type::value_type) complex_value_type;
typedef PROMOTE2(complex<float>,typename value_type_traits<value_type>::value_type) complex_type;
typedef typename complex_type::value_type real_type;
//typedef Vector<G> array_type;
//typedef typename array_type::size_type size_type;
//typedef PROMOTE2(complex<float>,typename Vector<G>::value_type) value_type;
//typedef typename value_type::value_type real_type;
size_type n = X.size();
assert(log2n(n)>=0);
for (size_t p_1=n/2, p=n, q=1; q<n; p=p_1, p_1/=2, q*=2)
{
complex_type w(1,0);
real_type theta = 2*PI/p;
real_type s = sin (theta);
real_type s2 = 2*sqr(sin(theta/2));
{
value_type t0;
for (size_t b=0; b<q; b++)
{
typename array_type::reference r0 = X[b*p];
typename array_type::reference r1 = X[b*p+p_1];
t0 = r0+r1;
r1 = r0-r1;
r0 = t0;
}
}
complex_value_type z0;
for (size_t a=1; a<p_1/2; ++a)
{
w+=s*complex_type(-w.imag(),w.real())-s2*w;
for (size_t b=0; b<q; ++b)
{
typename array_type::reference r0 = X[b*p+a];
typename array_type::reference r1 = X[b*p+p-a];
typename array_type::reference r2 = X[b*p+p_1-a];
typename array_type::reference r3 = X[b*p+p_1+a];
real(z0)=r0-r2; imag(z0)=r1+r3;
z0 *= w;
r0 += r2;
r2 = r1-r3;
r3 = real(z0);
r1 = imag(z0);
}
}
if (p_1>1) for (size_t b=0; b<q; b++) { X[b*p+p_1/2]*=2; X[b*p+p_1+p_1/2]*=-2; }
}
bitreverse_order(X);
X/=n;
} |
Partager