96
Jc[3] = reverse256[Ic[0]];
J >>= (32 - T);
if (I < J)
{
S = x[I];
x[I] = x[J];
x[J] = S;
}
}
//main loop
for(N = 2, Nd2 = 1, Skew = Nmax >> 1, Step= 1; N <= Nmax; Nd2 = N, N += N,
Skew >>= 1, Step++)
{
for(Warray = Wstore, k = 0; k < Nd2; k++, Warray += Skew)
{
for(m = k; m < Nmax; m += N)
{
Temp = *Warray;
if (complement)
Temp.im= -Temp.im;
mpNd2= m + Nd2;
Temp *= x[mpNd2];
x[mpNd2] = x[m];
x[mpNd2] -= Temp;
x[m] += Temp;
}
}
}
}
97
/*
x: x - array of items
N: N - number of items in array
complement: false - normal (direct) transformation, true - reverse transformation
*/
void universal_fft(ShortComplex *x, int N, bool complement)
{
ShortComplex *x_;
ShortComplex *w;
ShortComplex *Wstore;
int T;
T= (int)floor(log((double)N) / log(2.0) + 0.5);
if (1 << T == N)
{
Wstore= createWstore(N);
fft_step(x, T, complement, Wstore);
if (complement)
{
for(int i= 0; i < N; i++)
x[i]/= N;
}
delete [] Wstore;
return;
}
//find N', T
int N2= N+N;
int N_;
long double arg;
for(N_= 1, T= 0; N_ < N2; N_+= N_, T++)
{
}
98
//find --2pi/N/2 = pi/N
long double piN= M_PI / N;
if (complement)
piN= -piN;
//find x_[n] = x[n]*e^--2*j*pi*n*n/N/2 = x[n]*e^j*piN*n*n
x_= new ShortComplex[N_];
Complex v;
int n;
for(n= 0; n < N; ++n)
{
arg= piN*n*n;
v.re= cosl(arg);
v.im= sinl(arg);
complex_mul(x_ + n, x + n, &v);
}
for(; n < N_; ++n)
x_[n].re= x_[n].im= 0;
//find w[n] = e^-j*2*pi*(2*N-2-n)^2/N/2= e^-j*piN*(2*N-2-n)^2
w= new ShortComplex[N_];
int N22= 2*N - 2;
for(n= 0; n < N_; ++n, --N22)
{
arg= -piN*N22*N22;
w[n].re= (double)cos(arg);
w[n].im= (double)sin(arg);
}
//FFT1
Wstore= createWstore(N_);
fft_step(x_, T, false, Wstore);
//FFT2
fft_step(w, T, false, Wstore);
99
//svertka
for(n= 0; n < N_; ++n)
x_[n]*= w[n];
//FFT3 (complement)
fft_step(x_, T, true, Wstore);
//find X[n] = X_[n]*e^--j*2*pi*n*n/N/2 = X_[n]*e^j*piN*n*n
for(n= 0, N22= 2*N - 2; n < N; ++n, --N22)
{
arg= piN*n*n;
v.re= cosl(arg);
v.im= sinl(arg);
v/= N_;
if (complement)
v/= N;
complex_mul(x + n, x_ + N22, &v);
}
delete [] x_;
delete [] w;
delete [] Wstore;}