Commit 2da563b7 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

cleanup fft.

parent 699d33ae
......@@ -435,40 +435,32 @@ inverse_of_matrix(double **a, double **b, int n)
}
*/
void
fft(double *real, double *imag, int n, int sign)
fft(double *restrict real, double *restrict imag, int n, int sign)
{
/* n must be a power of 2 */
/* sign should be 1 (FT) or -1 (reverse FT) */
int i, j, j1, j2;
// n must be a power of 2
// sign should be 1 (FT) or -1 (reverse FT)
int j, j1, j2;
int bit;
int step;
double temp_r, temp_i, norm;
double w_r, w_i, ww_r, ww_i;
double temp_r, temp_i;
/* Bit reversal part */
for (i = j = 0; i < n; i++) /* The bit pattern of i and j are reverse */
// Bit reversal part
for (int i = j = 0; i < n; i++) // The bit pattern of i and j are reverse
{
if (i > j)
{
// swap real part
std::swap(real[i], real[j]);
// swap imaginary part
std::swap(imag[i], imag[j]);
}
if (i > j) std::swap(real[i], real[j]);
if (i > j) std::swap(imag[i], imag[j]);
for (bit = n >> 1; j & bit; bit >>= 1) j ^= bit;
j |= bit;
}
/* Danielson-Lanczos Part */
for (step = 1; step < n; step <<= 1)
// Danielson-Lanczos Part
for (int step = 1; step < n; step <<= 1)
{
w_r = std::cos(M_PI / step);
w_i = std::sin(M_PI / step) * sign;
ww_r = 1;
ww_i = 0;
for (i = 0; i < step; i++)
const double w_r = std::cos(M_PI / step);
const double w_i = std::sin(M_PI / step) * sign;
double ww_r = 1;
double ww_i = 0;
for (int i = 0; i < step; i++)
{
for (j1 = i, j2 = i + step; j2 < n; j1 += 2 * step, j2 += 2 * step)
{
......@@ -485,22 +477,17 @@ fft(double *real, double *imag, int n, int sign)
}
}
norm = 1. / std::sqrt(n);
for (i = 0; i < n; i++)
{
real[i] *= norm;
imag[i] *= norm;
}
const double norm = 1. / std::sqrt(n);
for (int i = 0; i < n; i++) real[i] *= norm;
for (int i = 0; i < n; i++) imag[i] *= norm;
}
void
ft(double *real, double *imag, int n, int sign)
{
/* sign should be 1 (FT) or -1 (reverse FT) */
int j, k;
// sign should be 1 (FT) or -1 (reverse FT)
static double *work_r = 0, *work_i = 0;
double sum_r, sum_i, norm;
double w_r, w_i, ww_r, ww_i, temp_r;
double temp_r;
if (!work_r)
{
......@@ -513,15 +500,15 @@ ft(double *real, double *imag, int n, int sign)
/* free_at_exit (work_i); */
}
for (k = 0; k < n; k++)
for (int k = 0; k < n; k++)
{
w_r = std::cos(2 * M_PI * k / n);
w_i = std::sin(2 * M_PI * k / n) * sign;
ww_r = 1;
ww_i = 0;
sum_r = 0;
sum_i = 0;
for (j = 0; j < n; j++)
const double w_r = std::cos(2 * M_PI * k / n);
const double w_i = std::sin(2 * M_PI * k / n) * sign;
double ww_r = 1;
double ww_i = 0;
double sum_r = 0;
double sum_i = 0;
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
sum_i += real[j] * ww_i + imag[j] * ww_r;
......@@ -533,32 +520,27 @@ ft(double *real, double *imag, int n, int sign)
work_i[k] = sum_i;
}
norm = 1. / std::sqrt(n);
for (k = 0; k < n; k++)
{
real[k] = work_r[k] * norm;
imag[k] = work_i[k] * norm;
}
const double norm = 1. / std::sqrt(n);
for (int k = 0; k < n; k++) real[k] = work_r[k] * norm;
for (int k = 0; k < n; k++) imag[k] = work_i[k] * norm;
}
/* reentrant version of ft */
// reentrant version of ft
void
ft_r(double *restrict real, double *restrict imag, int n, int sign, double *restrict work_r, double *restrict work_i)
{
/* sign should be 1 (FT) or -1 (reverse FT) */
int j, k;
double sum_r, sum_i, norm;
double w_r, w_i, ww_r, ww_i, temp_r;
// sign should be 1 (FT) or -1 (reverse FT)
double temp_r;
for (k = 0; k < n; k++)
for (int k = 0; k < n; k++)
{
w_r = std::cos(2 * M_PI * k / n);
w_i = std::sin(2 * M_PI * k / n) * sign;
ww_r = 1;
ww_i = 0;
sum_r = 0;
sum_i = 0;
for (j = 0; j < n; j++)
const double w_r = std::cos(2 * M_PI * k / n);
const double w_i = std::sin(2 * M_PI * k / n) * sign;
double ww_r = 1;
double ww_i = 0;
double sum_r = 0;
double sum_i = 0;
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
sum_i += real[j] * ww_i + imag[j] * ww_r;
......@@ -570,12 +552,9 @@ ft_r(double *restrict real, double *restrict imag, int n, int sign, double *rest
work_i[k] = sum_i;
}
norm = 1. / std::sqrt(n);
for (k = 0; k < n; k++)
{
real[k] = work_r[k] * norm;
imag[k] = work_i[k] * norm;
}
const double norm = 1. / std::sqrt(n);
for (int k = 0; k < n; k++) real[k] = work_r[k] * norm;
for (int k = 0; k < n; k++) imag[k] = work_i[k] * norm;
}
double
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment